在 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()，您可以更轻松地执行此操作，这是一个通用的提取位的 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)

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

Get Residual Variance-Covariance Matrix in lme4

r

lme4