lme4计算协方差的置信区间

Question

请参阅Ben Bolker 16/05/2016的答案,了解相应的解决方案.OP下面.

---

我正在使用lme4安装几个多级模型.我想报告随机效应的方差和协方差,并自动化这个过程.

我知道我可以得到差异as.data.frame(VarCorr(mymodel)),我知道我可以获得置信区间confint(mymodel).显然,我可以合并/绑定两个表,并通过简单地将输出confint()放在适当的行和列上来将置信区间放在方差周围,但是如果不是手动的话,我还无法找到一种令人信服的方法来计算协方差.

说结果confint是:

conf <- NULL
a <- c(6.2,-0.4,2.2,1.5,-0.4,-0.5,2.8,-0.9,1.3,3.9)
b <- c(6.8,-0.2,2.5,2.5,0.1,0.2,4.8,-0.7,2.3,5)
conf <- data.frame(a,b,row.names = c("sd_(Intercept)|ID","cor_Time.(Intercept)|ID","sd_Time|ID","sd_(Intercept)|Group","cor_Time.(Intercept)|Group","cor_I(Time^2).(Intercept)|Group","sd_Time|Group","cor_I(Time^2).Time|Group","sd_I(Time^2)|Group","sigma"))
colnames(conf) <- c("2.5%","97.5%")
conf

如何自动执行各种乘法以获得协方差

cov.time.intercept <- conf[1,2]*conf[1,1]*conf[1,3]

？

我试过分割标准偏差和相关性,创建"ID","时间","我(时间^ 2)"和"(拦截)"变量然后匹配两列,但我没有得到任何结果.问题是,每次模型更改时,您可能会有不同数量的方差和协方差,以及不同的三角矩阵.

感谢您的任何帮助,

ķ.

Answer 1

你的计算似乎给出了看似合理的答案，但它没有意义 （对我来说；我随时准备被纠正/启发......）。认为cov = corr*var1*var2。假设这ci(.)是某个数量的（下或上）置信限。这绝不是真的ci(cov) = ci(corr)*ci(var1)*ci(var2)（有趣的是你得到了合理的答案；我认为当数量近似不相关时最有可能发生这种情况......）如果你有每个成分的方差和它们之间的协方差（我并不意味着随机效应方差和协方差本身，而是它们的采样方差/协方差）您可以使用增量方法近似传播它们，但这些很难获得（请参见此处）。

据我所知，“正确”的方法是在方差-协方差尺度而不是标准差-相关尺度上进行似然剖面计算。这在以前是不可能的，但现在可以了（Github 上有开发版本）。

安装最新版本：

library(remotes) ## for install_github (or library(devtools))
install_github("lme4/lme4")

预赛：

chap12 <- foreign::read.dta(file = "ch12.dta")
library(lme4)
snijders <- lmer(prox_pup ~ 1 + prox_sel + (1 + occ|teacher),
                 data = chap12)

as.data.frame(VarCorr(snijders))
##        grp        var1 var2        vcov      sdcor
## 1  teacher (Intercept) <NA>  0.15617962  0.3951957
## 2  teacher         occ <NA>  0.01205317  0.1097869
## 3  teacher (Intercept)  occ -0.03883458 -0.8950676
## 4 Residual        <NA> <NA>  0.04979762  0.2231538

比较结果时我们必须小心一点，因为profile.merMod，我们很快就会使用它，自动（并且默默地！）将拟合从默认 REML 转换为最大似然拟合（因为基于 REML 的配置文件可能在统计上存在风险）；然而，这看起来并没有太大的区别。

s2 <- refitML(snijders)
as.data.frame(VarCorr(s2))
##        grp        var1 var2        vcov      sdcor
## 1  teacher (Intercept) <NA>  0.15426049  0.3927601
## 2  teacher         occ <NA>  0.01202631  0.1096645
## 3  teacher (Intercept)  occ -0.03884427 -0.9018483
## 4 Residual        <NA> <NA>  0.04955549  0.2226106

p.sd <- profile(s2,which="theta_",
              signames=FALSE)
p.vcov <- profile(s2,which="theta_",prof.scale="varcov",
              signames=FALSE)

我们收到一些关于非单调配置文件的警告......

confint(p.vcov)
##                                    2.5 %     97.5 %
## var_(Intercept)|teacher      0.08888931  0.26131067
## cov_occ.(Intercept)|teacher -0.07553263 -0.01589043
## var_occ|teacher              0.00000000  0.02783863
## sigma                        0.03463184  0.07258777

如果我们检查相关（标准差/方差）元素的平方会怎样？

confint(p.sd)[c(1,3,4),]^2
##                              2.5 %     97.5 %
## sd_(Intercept)|teacher 0.089089363 0.26130970
## sd_occ|teacher         0.002467408 0.02779329
## sigma                  0.034631759 0.07263869

除了方差的下限外，它们匹配得很好occ；它们也符合您上面的结果。然而，协方差结果（这是我声称很难的结果）对我来说是 (-0.0755,-0.0159)，对你来说是 (-0.0588,-0.0148)，大约有 20% 的差异。这可能不是什么大问题，具体取决于您想要做什么。

我们也尝试一下暴力破解：

sumfun <- function(x) {
    vv <- as.data.frame(VarCorr(x),order="lower.tri")[,"vcov"]
    ## cheating a bit here, using internal lme4 naming functions ...
    return(setNames(vv,
       c(lme4:::tnames(x,old=FALSE,prefix=c("var","cov")),
         "sigmasq")))
}

cc <- confint(s2,method="boot",nsim=1000,FUN=sumfun,seed=101,
        .progress="txt", PBargs=list(style=3))
## .progress/PBargs just cosmetic ...

##                                    2.5 %      97.5 %
## var_(Intercept)|teacher      0.079429623  0.24053633
## cov_occ.(Intercept)|teacher -0.067063911 -0.01479572
## var_occ|teacher              0.002733402  0.02378310
## sigmasq                      0.031952508  0.06736664

这里的“黄金标准”似乎介于我的配置文件结果和您的结果之间：协方差的下限是 -0.067，与 -0.0755（配置文件）或 -0.0588 相比。