获取lme4中的残差方差-协方差矩阵 [英] Get Residual Variance-Covariance Matrix in lme4

问题描述

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

I am fitting a linear mixed effects model using 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).

If we say generically Y = X * B + Z * d + e is the form of a linear mixed effects model, then I am trying to get Var(Y) = Z * Var(d) * Z^t + Var(e) from the results of the model.

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

Is the following formulation the right way to do this?

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

我认为这是一个简单的问题，但是我在任何地方都找不到执行此操作的代码，所以我问我是否做对了.

I think this is a simple question, but I can't find anywhere code for doing this, so I'm asking if I'm doing it right.

推荐答案

如果您知道getME()，则可以更轻松地完成此操作，这是通用的extract-bits-a- lmer-拟合函数.特别是，您可以提取转置的Z矩阵(getME(.,"Zt"))和转置的Lambda矩阵-Lambda矩阵是条件模型(BLUP)的缩放方差-协方差矩阵的Cholesky因子；用您的符号表示，Var(d)是剩余方差乘以Lambda的叉积.

You can do this a bit more easily if you know about getME(), which is a general purpose extract-bits-of-a-lmer-fit function. In particular, you can extract the transposed Z matrix (getME(.,"Zt")) and the transposed Lambda matrix - the Lambda matrix is the Cholesky factor of the scaled variance-covariance matrix of the conditional models (BLUPs); in your notation, Var(d) is the residual variance times the cross-product of Lambda.

在此处引用的答案很好，但是下面的答案稍微更笼统(适用于任何lmer适合).

The answer cited here is pretty good but the answer below is slightly more general (it should work for any lmer fit).

适合的型号:

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中的残差方差-协方差矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！