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

Question

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

6

我正在使用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

我认为这是一个简单的问题，但我找不到任何可以完成这个任务的代码，所以我想问一下我的操作是否正确。

- learner

运行您的代码会抛出一个错误，因为您引用了一个未定义的 x 对象（在sI <- diag(sigma^2, nrow=nrow(x))中）。无论如何，您可以使用 vcov 方法获取任何 merMod 对象（即 lmer 生成的对象）的方差-协方差矩阵。在您的情况下：vcov(d.test)。但是这只提供了固定效应的矩阵。 - Oriol Mirosa

@OriolMirosa 感谢您指出错误。但是，vcov/VarCorr 不是我正在寻找的答案 - 它们为固定和随机效应提供方差/协方差。我试图获取数据 Y 的方差/协方差。 - learner

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

如果您了解getME()，这个通用的提取一个lmer拟合函数的位的函数，那么您可以更加容易地完成这个操作。特别是，您可以提取转置Z矩阵 (getME(.,"Zt")) 和转置Lambda矩阵 - Lambda矩阵是条件模型（BLUPs）的缩放方差-协方差矩阵的Cholesky因子；在您的表示中，Var(d)是剩余方差乘以Lambda的交叉积。

这里引用的答案相当不错，但下面的答案稍微更加通用（它适用于任何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)