如何从 lme4 获取随机效应（BLUP/条件模式）的协方差矩阵答案

【问题标题】：How to get covariance matrix for random effects (BLUPs/conditional modes) from lme4如何从 lme4 获取随机效应（BLUP/条件模式）的协方差矩阵
【发布时间】：2018-11-11 11:11:46
【问题描述】：

所以，我在 R 中拟合了一个带有两个随机截距的线性混合模型：

Y = X beta  + Z b + e_i,

b ~ MVN (0, Sigma); X 和 Z 分别是固定效应和随机效应模型矩阵，beta 和 b 是固定效应参数和随机效应 BLUP/条件模式。

我想了解b 的底层协方差矩阵，这在lme4 包中似乎不是一件小事。您只能得到VarCorr 的方差，而不是实际的相关矩阵。

根据one of the package vignettes（第2页），可以计算出beta的协方差：e_i * lambda * t(lambda)。您可以从lme4 的输出中提取所有这些组件。

我想知道这是否是要走的路？或者您有什么其他建议？

【问题讨论】：

欢迎来到 StackOverflow。请澄清您的数学符号（Xbeta 实际上是 X * beta 吗？也许您还应该说出 beta、b、sigma 是什么；尽管我不是专家，对于某些这些符号可能很明显）。请记住，您写的问题越清晰，就越有可能更快地得到适当的答案。
是的，Xbeta 应该是 X*beta。 Beta是设计矩阵X的固定效应向量，b是随机效应向量，sigma是b的方差协方差矩阵。感谢您的提示，我会记住这一点的。

标签： r matrix covariance lme4

【解决方案1】：

来自?ranef：

如果“condVar”为“TRUE”，则每个数据帧都有一个属性称为“postVar”，它是一个对称的三维数组面孔；每个面都包含一个方差-协方差矩阵分组因子的特定水平。（此属性的名称是历史文物，有时可能会更改为“condVar” 指向未来。）

建立一个例子：

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
rr <- ranef(fm1,condVar=TRUE)

获取截距的b 值之间的方差-协方差矩阵

pv <- attr(rr[[1]],"postVar")
str(pv)
##num [1:2, 1:2, 1:18] 145.71 -21.44 -21.44 5.31 145.71 ...

所以这是一个 2x2x18 数组；每个切片是特定受试者的条件截距和斜率之间的方差-协方差矩阵（根据定义，每个受试者的截距和斜率独立于所有其他受试者的截距和斜率）。

要将其转换为方差-协方差矩阵（请参阅getMethod("image",sig="dgTMatrix") ...）

library(Matrix)
vc <- bdiag(  ## make a block-diagonal matrix
            lapply(
                ## split 3d array into a list of sub-matrices
                split(pv,slice.index(pv,3)),
                   ## ... put them back into 2x2 matrices
                      matrix,2)) 
image(vc,sub="",xlab="",ylab="",useRaster=TRUE)

【讨论】：

非常感谢！我去看看:)
StackOverflow 弃用 using comments to say "thank you";如果此答案有用，您可以投票（如果您有足够的声誉），并且无论如何如果它令人满意地回答了您的问题，我们鼓励您单击复选标记以接受它。
抱歉，可能需要一段时间。但我试图理解你的话each slice is the variance-covariance matrix among the conditional intercept and slope for a particular subject。那么每个切片到底是什么？你是说海报问的b 的协方差吗？
我问的原因是因为我试图使用你的代码来估计b的协方差矩阵，但结果与我们实际b的协方差有很大不同用于生成数据。非常感谢！
使用上面的代码获取协方差矩阵和从lmer 的输出中读取Std.Dev. 有什么区别？根据stats.stackexchange.com/questions/24452/…，它们应该是同一个意思，但在我的实验中它们的值并不相同。