lmer 输出的方差分量的标准误差

Question

我需要从

standard error

的输出中提取方差分量

lmer

。

library(lme4)
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)

以下产生方差分量的估计：

s2 <- VarCorr(model)$Subject[1]

它不是方差的标准误差。我想要标准错误。我怎样才能拥有它？

编辑：

也许我无法让您理解“方差分量的标准误差”的含义。所以我正在编辑我的帖子。

Douglas C. Montgomery 所著的实验设计与分析一书中的第 12 章“随机因素实验”在该章的末尾，例 12-2 是由 SAS 完成的。在例12-2中，模型是一个二因素阶乘随机效应模型。输出在表12-17中给出

我正在尝试通过

lmer

来拟合 R 中的模型。

library(lme4)
fit <- lmer(y~(1|operator)+(1|part),data=dat)

用于提取

Estimate

的 R 代码，在表 12-17 中用 4 注释：

est_ope=VarCorr(fit)$operator[1]
est_part = VarCorr(fit)$part[1]
sig = summary(fit)$sigma
est_res = sig^2

现在我想从 lmer 输出中提取

Std Errors

的结果，在表 12-17 中用 5 注释。

非常感谢！

Answer 1

我认为您正在寻找方差估计的 Wald 标准误差。请注意，这些（正如 Doug Bates 经常指出的那样）Wald 标准误差通常是对方差不确定性的“非常差”估计，因为似然分布在方差尺度上通常远非二次......我假设你知道自己在做什么并且对这些数字有一些很好的用途...... 这可以（现在）使用

merDeriv

包来完成。

library(lme4)
library(merDeriv)
m1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
sqrt(diag(vcov(m1, full = TRUE)))
vv <- vcov(m1, full = TRUE)
colnames(vv)
## [1] "(Intercept)"                  "Days"                        
## [3] "cov_Subject.(Intercept)"      "cov_Subject.Days.(Intercept)"
## [5] "cov_Subject.Days"             "residual"

因为这代表完整/组合协方差矩阵，所以前两个索引（行/列）{（截距），天数}代表固定效应截距和斜率。以 

cov_

开头的元素是方差分量，格式为

cov_<grp_variable>.<term>

表示方差（例如

cov_Subject.(Intercept)

是截距中的受试者间方差），

cov_<grp_variable>.<term1>.<term2>

表示协方差（

cov_Subject.Days.(Intercept)

是受试者间协方差）斜率 (

Days

) 和截距之间）。

residual

表示残差方差。如果我们想要方差分量的标准误差，我们取对角线的平方根并仅保留元素 3 到 5：

sqrt(diag(vv)[3:5])
## [1] 288.78602  46.67876  14.78208

或更笼统地说

sqrt(diag(vv)[grepl("^cov", colnames(vv))])

（奇怪的是，只有 

colnames()

有效 -

vv

的行名称对于随机效应组件是空白的）。旧答案

library("lme4")
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy, REML=FALSE)

（目前对于 REML 估计来说做到这一点相当困难......）

提取以标准差和相关性而不是 Cholesky 因子参数化的偏差函数（请注意，这是一个内部函数，因此不能保证它将来会继续以相同的方式工作......）

dd.ML <- lme4:::devfun2(model,useSc=TRUE,signames=FALSE)

提取参数作为原始尺度上的标准差：

 vv <- as.data.frame(VarCorr(model)) ## need ML estimates!
 pars <- vv[,"sdcor"]
 ## will need to be careful about order if using this for
 ## a random-slopes model ...

现在计算二阶导数（Hessian）矩阵：

library("numDeriv")
hh1 <- hessian(dd.ML,pars)
vv2 <- 2*solve(hh1)  ## 2* converts from log-likelihood to deviance scale
sqrt(diag(vv2))  ## get standard errors

这些是标准差的标准误差：将它们加倍以获得方差的标准误差（当您变换一个值时，其标准误差根据变换的导数缩放）。

我认为这应该可以，但你可能需要仔细检查......

Answer 2

library(arm)
se.ranef(model)
#$Subject
#    (Intercept)
#308    9.475668
#309    9.475668
#310    9.475668
#330    9.475668
#331    9.475668
#332    9.475668
#333    9.475668
#334    9.475668
#335    9.475668
#337    9.475668
#349    9.475668
#350    9.475668
#351    9.475668
#352    9.475668
#369    9.475668
#370    9.475668
#371    9.475668
#372    9.475668

这实际上是随机效应的条件方差-协方差矩阵的平方根：

sqrt(attr(ranef(model, condVar = TRUE)$Subject, "postVar"))

Answer 3