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

Question

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

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

以下是差异部分的估计数：

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

这是，而不是，方差的标准误差。我要标准错误。我怎么能得到它？

编辑：

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

在第12章，随机因素的实验，由道格拉斯C.蒙哥马利(DouglasC.Montgomery)著的“实验设计与分析”一书中，在本章末尾，示例12-2是由SAS完成的。在例子12-2中，该模型是一个双因子因子随机效应模型，.The输出如表12-17所示。

​

​

我正试图在R中拟合lmer的模型。

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

R提取Estimate的代码，表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经常指出的) 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

这些是标准差的标准误差:将它们加倍以得到方差的标准误差(当你转换一个值时，它的标准误差尺度根据转换的导数)。

我觉得这应该行，但你可能想再查一遍.

Answer 2

我不太清楚你所说的“方差分量的标准误差”是什么意思。我的最佳猜测(基于您的代码)是，您希望获得随机效果的标准错误。您可以使用package获得这个信息：

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

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