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

【问题标题】：Standard Error of variance component from the output of lmerlmer 输出的方差分量的标准误差
【发布时间】：2015-07-29 08:11:55
【问题描述】：

我需要从lmer 的输出中提取方差分量的standard error。

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

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

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

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

编辑：

也许我无法让您理解我所说的“方差分量的标准误差”是什么意思。所以我正在编辑我的帖子。

在 Douglas C. Montgomery 的书 Design and Analysis of Experiments 的第 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 注释。

非常感谢！

【问题讨论】：

你能发布这个例子的数据吗？这也有助于了解您将使用标准误差来做什么（正如我在下面指出的那样，方差估计的标准误差是不可靠的不确定性指标 - 配置文件间隔会更好）
为什么你想要一个非对称分布的参数的标准误差。您应该重新构建您提出的问题。不要模仿 SAS 的错误。如果要进行低测试，请使用 anove 功能。如果您需要 CI，请使用配置文件或引导 CI。 lmer 不提供您要求的数字是有原因的。虽然本告诉你如何得到它。不要让 SAS 或 Stata 报告影响您的事实。

标签： r lme4

【解决方案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"

如果我们想要方差分量的标准误，我们取对角线的平方根，只保留最后三个元素：

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

旧答案

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

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

我认为应该这样做，但您可能需要仔细检查...

【讨论】：

能否请您给我一份 REML 估算的参考资料。表中的结果是针对方差分量的 REML 估计的标准误差。也许这就是为什么我没有得到答案。我已经正确指定了模型，并通过VarCorr(model) 获得了方差分量的 REML 估计值，它与表中的结果相匹配。如果您向我介绍如何计算 lmer(Reaction ~ Days + (Days||Subject), sleepstudy) 标准差的标准误差，那将非常有帮助
在 lmer 中是否可以计算 REML 估计值的标准差的标准误？
如果我写REML=TRUE 并以同样的方式继续，它不会工作吗？
我想可能不会。我可能会在github.com/lme4/lme4/issues 上打开一个问题......有几种方法可以做到这一点，但没有非常简单的解决方案。

【解决方案2】：

我不太确定“方差分量的标准误差”是什么意思。我最好的猜测（根据您的代码）是您想要随机效应的标准误差。您可以使用 package arm 来获得它：

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"))

【讨论】：

【解决方案3】：

mn2=lmer(pun~ pre + (pre|pro), REML = TRUE, data = pro)
summary(mn2)
coe2=coef(mn2)
coe2

# Matriz de varianza-covarianza (covarianza)
as.data.frame(VarCorr(mn2))

# Extraer coeficientes fijos
fixef(mn2)

# Extraer desvios de a - alfa y b - beta
re=as.data.frame(ranef(mn2))

【讨论】：

向随机效应选项 (ranef) 询问 as.data.frame 转换为长格式 condsd 对应条件标准差