带交互项的 lmer 引导答案

【问题标题】：bootstrapping for lmer with interaction term带交互项的 lmer 引导
【发布时间】：2017-03-28 22:17:54
【问题描述】：

我在 R 中使用 lme4 运行混合模型：

full_mod3=lmer(logcptplus1 ~ logdepth*logcobb + (1|fyear) + (1 |flocation),
data=cpt, REML=TRUE)

总结：

Formula: logcptplus1 ~ logdepth * logcobb + (1 | fyear) + (1 | flocation)
   Data: cpt

REML criterion at convergence: 577.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7797 -0.5431  0.0248  0.6562  2.1733 

 Random effects:
 Groups    Name        Variance Std.Dev.
 fyear     (Intercept) 0.2254   0.4748  
 flocation (Intercept) 0.1557   0.3946  
 Residual              0.9663   0.9830  
Number of obs: 193, groups:  fyear, 16; flocation, 16

Fixed effects:
                  Estimate Std. Error t value
(Intercept)        4.3949     1.2319   3.568
logdepth           0.2681     0.4293   0.625
logcobb           -0.7189     0.5955  -1.207
logdepth:logcobb   0.3791     0.2071   1.831

我使用 R 中的 effects 包和函数来计算模型输出的 95% 置信区间。我已经使用effects 包计算并提取了 95% CI 和标准误差，这样我就可以通过将次要预测变量 (logdepth) 保持在中位数来检查重要的预测变量和响应变量之间的关系(2.5)在数据集中：

gm=4.3949 + 0.2681*depth_median + -0.7189*logcobb_range + 0.3791*
(depth_median*logcobb_range)

ef2=effect("logdepth*logcobb",full_mod3,
     xlevels=list(logcobb=seq(log(0.03268),log(0.37980),,200)))

我尝试使用来自here 的代码引导 95% 的 CI。但是，我只需要计算中值深度 (2.5) 的 95% CI。有没有办法在 confint() 代码中指定，以便我可以计算可视化引导结果所需的 CI，如上图所示？

confint(full_mod3,method="boot",nsim=200,boot.type="perc")

【问题讨论】：

标签： r lme4 statistics-bootstrap

【解决方案1】：

您可以通过指定自定义函数来做到这一点：

library(lme4)
?confint.merMod

FUN：引导函数；如果为“NULL”，将使用返回固定效应参数以及标准偏差/相关尺度上的随机效应参数的内部函数。有关详细信息，请参阅“bootMer”。

因此，FUN 可以是一个预测函数 (?predict.merMod)，它使用一个 newdata 参数来改变和修复适当的预测变量。

一个内置数据的例子（不像你的那样有趣，因为有一个连续的预测变量，但我认为它应该足够清楚地说明这种方法）：

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
pframe <- data.frame(Days=seq(0,20,by=0.5))
## predicted values at population level (re.form=NA)
pfun <- function(fit) {
    predict(fit,newdata=pframe,re.form=NA)
}
set.seed(101)
cc <- confint(fm1,method="boot",FUN=pfun)

图片：

par(las=1,bty="l")
matplot(pframe$Days,cc,lty=2,col=1,type="l",
     xlab="Days",ylab="Reaction")

【讨论】：