r - 在 lme4 中获取残差-协方差矩阵

Question

我正在使用 lme4 拟合线性混合效应模型:

library(lme4)
data(Orthodont)
dent <- Orthodont
d.test <- lmer(distance ~ age + (1|Subject), data=dent)

如果我们一般说 Y = X * B + Z * d + e是线性混合效应模型的形式，然后我试图得到 Var(Y) = Z * Var(d) * Z^t + Var(e)从模型的结果来看。

以下公式是正确的方法吗？

k <- table(dent$Subject)[1]
vars <- VarCorr(d.test)
v <- as.data.frame(vars)
sigma <- attr(vars, "sc")
s.tech <- diag(v$vcov[1], nrow=k)
icc <- v$vcov[1]/sum(v$vcov)
s.tech[upper.tri(s.tech)] <- icc
s.tech[lower.tri(s.tech)] <- icc
sI <- diag(sigma^2, nrow=length(dent$age))
var.b <- kronecker(diag(1, nrow=length(dent$age)/k), s.tech)
var.y <- sI + var.b

我认为这是一个简单的问题，但我在任何地方都找不到这样做的代码，所以我在问我是否做得对。

Answer 1

如果您了解 getME()，您可以更轻松地完成此操作，这是一个通用的提取位-a- lmer - 拟合功能。特别是，您可以提取转置 Z 矩阵 ( getME(.,"Zt") ) 和转置 Lambda 矩阵 - Lambda 矩阵是条件模型 (BLUP) 的缩放方差-协方差矩阵的 Cholesky 因子；在您的符号中，Var(d)是剩余方差乘以 Lambda 的叉积。

答案引用here非常好，但下面的答案更笼统(它应该适用于任何 lmer 适合)。

适合型号:

library(lme4)
data(Orthodont,package="nlme")
d.test <- lmer(distance ~ age + (1|Subject), data=Orthodont)

提取成分:

var.d <- crossprod(getME(d.test,"Lambdat"))
Zt <- getME(d.test,"Zt")
vr <- sigma(d.test)^2

组合它们:

var.b <- vr*(t(Zt) %*% var.d %*% Zt)
sI <- vr * Diagonal(nrow(Orthodont))
var.y <- var.b + sI

照片:

image(var.y)