lme4计算协方差的置信区间

Question

请参考Ben Bolker在2016年5月16日的答案，以获得适当的解决方案。以下是OP。

---

我正在使用lme4拟合几个多层模型。我想报告随机效应的方差和协方差，并自动化这个过程。

我知道可以使用as.data.frame(VarCorr(mymodel))获取方差，也知道可以使用confint(mymodel)获取置信区间。显然，我可以合并/ rbind两个表格，并通过简单地将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”、“时间”、“I(时间^2)”和“(截距)”变量，然后通过两列匹配，但我一无所获。问题是每次模型更改时，方差和协方差的数量可能不同，三角矩阵也可能不同。

谢谢任何帮助，

k.

- r.kaiza

---

你能再清楚地说明一下你想做什么吗？你想要方差和协方差的置信区间吗？还是你想要方差和协方差本身的方差和协方差？正如@Thierry在下面建议的那样，我认为在尝试找到正确的计算框架之前，你需要先解决/澄清一些问题。 - Ben Bolker

嗨，本，感谢您的回复，如果那不清楚，我很抱歉。实际上，我想将置信区间表示为方差和协方差，而不是作为lme4 :: confint()默认情况下的标准偏差和相关性。 - r.kaiza

只是一个快速的提示，获取方差-协方差尺度上的随机效应的置信区间确实不容易；我正在对该软件包进行一些修改，以帮助解决这个问题。 - Ben Bolker

嗨，本，我不明白，难道不只是将标准化结果乘以SD再平方/相乘的问题吗？最好的祝愿。 - r.kaiza

不要返回个人资料置信区间（仍在努力）。 - Ben Bolker

Answer 1

问题已解决，感谢您的贡献。我将更新初始帖子。您可以使用Snijders & Bosker提供的数据集进行测试，数据集在此处。

导入方式：

library(foreign)
chap12 <- read.dta(file = "<your path>/ch12.dta")

一个临时的模型：

snijders <- lmer(prox_pup ~ 1 + prox_sel + (1 + occ|teacher), data = chap12)

引入该函数：

ExtractVarCovCI <- function(Model) {

v <- NULL
v <- as.data.frame(VarCorr(Model),order = "lower.tri") #Extract variances and covariances

conf <- confint(Model, parm  ="theta_", oldNames = F) #extract CIs

v.conf <- cbind(v,conf) #bind confidence intervals

covs <- as.data.frame(v.conf[!is.na(v[,3]),]) #separate variance from covariance components
vars <- as.data.frame(v.conf[is.na(v[,3]),]) #separate variance from covariance components
vars.sq <- vars[,6:7]^2 #calculate square of variance components
colnames(vars.sq) <- sub("[%]", "% sq.", colnames(vars.sq))

vars2 <- cbind(vars,vars.sq) #bind squares of variance components
covs$`2.5 % sq.` <- c(rep(NA,nrow(covs))) #create empty columns for later
covs$`97.5 % sq.` <- c(rep(NA,nrow(covs))) #create empty columns for later

lcovs <- length(row.names(covs)) #now we re-organise the table so that each covariance is below the variance of its variables
k <- NULL
for (i in seq(1:lcovs)) {
  k <- rbind(k,vars2[vars2$grp %in% covs[i,1] & vars2$var1 %in% covs[i,2],],vars2[vars2$grp %in% covs[i,1] & vars2$var1 %in% covs[i,3],],covs[i,])
}

k2 <- rbind(k,vars2["sigma",]) #bind the level-1 residuals at the end

k2.covrow <- grep("^cor",rownames(k2)) # isolate covariance row position
k2[k2.covrow,8] <- k2[k2.covrow,6]*k2[k2.covrow-1,6]*k2[k2.covrow-2,6] #calculate covariance 2.5%
k2[k2.covrow,9] <- k2[k2.covrow,7]*k2[k2.covrow-1,7]*k2[k2.covrow-2,7] #calculate covariance 97.5%

p <- NULL
p <- k2[,c(4,8:9)] #retain only the estimates and the confidence intervals
rownames(p) <- sub("^sd","var",rownames(p)) #now it's clear that we have proper variances and covariances
rownames(p) <- sub("^cor","cov",rownames(p)) #now it's clear that we have proper variances and covariances
colnames(p) <- c("Estimate", "2.5%", "97.5%")

return(p)
}

运行函数：

ExtractVarCovCI(snijders)

我的输出结果是：

                               Estimate         2.5%       97.5%
var_(Intercept)|teacher      0.15617962  0.089020350  0.26130969
var_occ|teacher              0.01205317  0.002467408  0.02779329
cov_occ.(Intercept)|teacher -0.03883458 -0.014820577 -0.05887660
sigma                        0.04979762  0.034631759  0.07263837

现在我们有一个方差-协方差表，其中使用非标准化的随机效应及其上下置信区间。我确信有更好的方法来做到这一点，但这是一个开始...

好的。

- r.kaiza

Answer 2

你的计算似乎给出了合理的答案，但对我来说它不合理（我准备接受更正/启示...）。假设cov = corr * var1 * var2。假设ci(.)是一个量的（下限或上限）置信区间。这并不意味着ci(cov) = ci(corr) * ci(var1) * ci(var2)（有趣的是你得到了合理的答案; 我认为这最可能发生在数量大致不相关的情况下...）。如果您有每个组件的方差和它们之间的协方差（我不是指随机效应的方差和协方差本身，而是它们的抽样方差/协方差），则可以使用delta方法近似地传播它们，但这些很难获得（请参见此处）。

“正确”的做法，据我所知，是在方差-协方差尺度上进行似然轮廓计算，而不是在标准差-相关性尺度上进行。这以前是不可能的，但现在可以通过 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

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

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%。根据您的目的，这可能并不重要。

让我们也尝试 brute force 方法：

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。

- Ben Bolker

---

谢谢你的深入和有见地的回答，Ben。为了比较起见，使用引导法（nsim = 50）我得到了(-0.014,-0.042)的结果，但如果你的方法更具理论基础，那么我一定会使用它并在我的OP中加上一个便签。也非常感谢在软件包源代码中实现这一点的人。K. - r.kaiza

Ben，我在更新到lme4 1.1-13之后尝试了你的输入，但是我收到了以下错误：s2 <- refitML(snijders) Error in assign(field, forceCopy(current), envir = vEnv) : could not find function "forceCopy"我是否缺少某个包？谢谢，k。 - r.kaiza

哎呀，这很奇怪。 - Ben Bolker

我可以提供任何信息来帮助您跟踪它卡住的位置吗？ - r.kaiza

sessionInfo() 的结果是什么？ - Ben Bolker

Answer 3

请注意，lme4 摘要中的随机效应标准差不是方差的标准误差！它只是方差的平方根！

如果您需要随机效应方差的置信区间，则需要对似然函数进行profile()。请参见?lme4::profile。

- Thierry

---

嗨，Thierry，感谢您的回复。我知道随机效应的标准差不是方差的标准误差，而是平方根。这就是为什么我正在使用confint()计算置信区间。然而，置信区间报告的是随机效应的平方根。对于方差来说，将它们平方是有意义的，但对于协方差，您必须使用标准差，并且我正在尝试自动化这个过程。 - r.kaiza

这似乎是一条注释，不是原帖作者想要的（现在他们已经澄清了）... - Ben Bolker