在 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()，您可以更轻松地执行此操作，这是一个通用的 extract-bits-of-a-lmer-fit 函数。特别是，您可以提取转置的 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)