lme4 计算协方差的置信区间
lme4 calculate confidence intervals of covariances
请参阅 Ben Bolker 16/05/2016 的回答以获得适当的解决方案。下面的OP。
我正在用 lme4 拟合多个多级模型。我想报告随机效应的方差和协方差，并自动执行此过程。
我知道我可以用 as.data.frame(VarCorr(mymodel)) 得到方差，我知道我可以用 confint(mymodel) 得到置信区间。显然，我可以 merge/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"、"Time"、"I(Time^2)" 和“(Intercept)”变量，然后按两列进行匹配，但我没有到达任何地方。问题是每次模型更改时，您可能会有不同数量的方差和协方差，以及不同的三角矩阵。
感谢您的帮助，
k.
请注意 lme4 摘要中随机效应的标准差 不是 方差的标准差！只是方差的平方根！
如果您需要随机效应方差的置信区间，则需要 profile() 似然。参见 ?lme4::profile。
已解决，感谢贡献。我会更新初始post。可以使用 Snijders & Bosker 提供的数据集测试结果 here.
导入
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
现在我们有一个方差-协方差 table，它使用非标准化随机效应及其上下置信边界。我相信有更好的方法可以做到这一点，但这是一个开始...
k.
你的计算似乎给出了合理的答案，但它没有意义（对我来说；我随时准备成为 corrected/enlightened ...）。假设cov = corr*var1*var2。假设 ci(.) 是数量的（下限或上限）置信限度。 ci(cov) = ci(corr)*ci(var1)*ci(var2) 绝不是真的（它 是 有趣的是你得到了合理的答案；我认为这最有可能发生在数量近似不相关的时候......）如果你有每个分量的方差和它们之间的协方差（我不是指随机效应方差和协方差本身，而是 他们的 采样 variances/covariances）你可以传播他们 大约 使用 delta 方法，但这些很难获得（参见 here）。
据我所知，"right" 方法是在方差-协方差尺度而不是标准偏差-相关尺度上进行似然曲线计算。这在以前是不可能的，但现在是（开发版本在 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/variance) 个元素的平方会怎样？
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
此处的 "gold standard" 似乎介于我的配置文件结果和您的结果之间：此处的协方差下限为 -0.067 与 -0.0755（配置文件）或 -0.0588。
请参阅 Ben Bolker 16/05/2016 的回答以获得适当的解决方案。下面的OP。
我正在用 lme4 拟合多个多级模型。我想报告随机效应的方差和协方差，并自动执行此过程。
我知道我可以用 as.data.frame(VarCorr(mymodel)) 得到方差，我知道我可以用 confint(mymodel) 得到置信区间。显然，我可以 merge/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"、"Time"、"I(Time^2)" 和“(Intercept)”变量，然后按两列进行匹配，但我没有到达任何地方。问题是每次模型更改时，您可能会有不同数量的方差和协方差，以及不同的三角矩阵。
感谢您的帮助，
k.
请注意 lme4 摘要中随机效应的标准差 不是 方差的标准差！只是方差的平方根！
如果您需要随机效应方差的置信区间，则需要 profile() 似然。参见 ?lme4::profile。
已解决，感谢贡献。我会更新初始post。可以使用 Snijders & Bosker 提供的数据集测试结果 here.
导入
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
现在我们有一个方差-协方差 table，它使用非标准化随机效应及其上下置信边界。我相信有更好的方法可以做到这一点，但这是一个开始...
k.
你的计算似乎给出了合理的答案，但它没有意义（对我来说；我随时准备成为 corrected/enlightened ...）。假设cov = corr*var1*var2。假设 ci(.) 是数量的（下限或上限）置信限度。 ci(cov) = ci(corr)*ci(var1)*ci(var2) 绝不是真的（它 是 有趣的是你得到了合理的答案；我认为这最有可能发生在数量近似不相关的时候......）如果你有每个分量的方差和它们之间的协方差（我不是指随机效应方差和协方差本身，而是 他们的 采样 variances/covariances）你可以传播他们 大约 使用 delta 方法，但这些很难获得（参见 here）。
据我所知，"right" 方法是在方差-协方差尺度而不是标准偏差-相关尺度上进行似然曲线计算。这在以前是不可能的，但现在是（开发版本在 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/variance) 个元素的平方会怎样？
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
此处的 "gold standard" 似乎介于我的配置文件结果和您的结果之间：此处的协方差下限为 -0.067 与 -0.0755（配置文件）或 -0.0588。