使用 R 中的 confint 计算固定效应的 CI

Question

我想执行引导以获得二项式 GLMM 中固定效应的 95% 顺式：

m <- glmer(cbind(df$Valid.detections, df$Missed.detections) ~ distance + 
              Habitat + Replicate + transmitter.depth + receiver.depth + 
              wind.speed + wtc + Transmitter + (1 | Unit) + 
              (1 | SUR.ID) + distance:Transmitter + 
              distance:Habitat + distance:transmitter.depth + distance:receiver.depth + 
              distance:wind.speed, data = df, family = binomial(link=logit),control=glmerControl(calc.derivs=F))

我发现confint()函数可以实现这个，所以我指定了函数：

confint(m, method = "boot", boot.type = "basic", seed = 123, nsim = 1000)

在我决定终止之前，该函数已经运行了 8 个多小时。终止后，我收到以下警告消息 (x10)：

Warning messages:
1: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations

我的问题是： 1) 我需要担心这些警告信息吗？如果是这样，我该如何处理它们？，2）因为 8 小时后它仍在运行，我不知道执行此功能需要多长时间。因此，在执行此功能时最好有某种进度条。我读到 confint() 可以从 bootMer 获取参数，所以我包含了参数 .progress = "txt"，结果是：

confint(m, method = "boot", boot.type = "basic", seed = 123, nsim = 1000, .progress = "txt")

但它似乎不起作用。或者，如果有更好的方法来实现相同的目标，我愿意接受建议。

谢谢你的帮助

Answer 1

下面是一个例子：

library("lme4")
(t1 <- system.time(
    gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                 data = cbpp, family = binomial)))
##    user  system elapsed 
##   0.188   0.000   0.186

nranpars <- length(getME(gm1,"theta"))
nfixpars <- length(fixef(gm1))

(t2 <- system.time(c1 <- confint(gm1,method="boot", nsim=1000,
                  parm=(nranpars+1):(nranpars+nfixpars),
                  .progress="txt")))

##    user  system elapsed 
## 221.958   0.164 222.187

请注意，此进度条只有 80 个字符长，因此它仅在每 1000/80=12 次引导迭代后递增。如果您的原始模型需要一个小时才能适应，那么您不应该期望在 12 小时后看到任何进度条活动......

(t3 <- system.time(c2 <- confint(gm1,
                  parm=(nranpars+1):(nranpars+nfixpars))))

##    user  system elapsed 
##   5.212   0.012   5.236 

1000 次 bootstrap 代表可能有点过头了——如果你的模型拟合很慢，你可能可以从 200 次 bootstrap 代表中获得合理的CI。

我也试过这个optimizer="nloptwrap"，希望它能加快速度。确实如此，尽管有一个缺点（见下文）。

(t4 <- system.time(
    gm1B <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                 data = cbpp, family = binomial, 
                 control=glmerControl(optimizer="nloptwrap"))))
##   user  system elapsed 
##  0.064   0.008   0.075 

(t5 <- system.time(c3 <- confint(gm1B,method="boot", nsim=1000,
                  parm=(nranpars+1):(nranpars+nfixpars),
                  .progress="txt",PBargs=list(style=3))))
##
##   user  system elapsed 
## 65.264   2.160  67.504

这要快得多，但会发出有关收敛的警告（在本例中为 37）。根据 ，以all.equal()这种方式计算的置信区间只有大约 2% 的差异。（包装本身还有一些皱纹需要整理……）

加速这一过程的最佳选择是并行化——不幸的是，这样你就失去了使用进度条的能力......

(t6 <- system.time(c4 <- confint(gm1,method="boot", nsim=1000,
                  parm=(nranpars+1):(nranpars+nfixpars),
                  parallel="multicore", ncpus=4)))

## 
##     user  system elapsed 
##  310.355   0.916 116.917

这需要更多的用户时间（它计算所有内核上使用的时间），但经过的时间减少了一半。（用 4 个内核做得更好会更好，但两倍的速度仍然很好。这些是虚拟 Linux 机器上的虚拟内核，真实（非虚拟）内核可能具有更好的性能。）

与nloptwrap和多核组合的I可以得到时间降低至91秒（用户）/ 36秒（消逝）。