如何展示Bootstrapping GLM结果中各因子水平的置信区间？

问题背景 #

对Gamma分布的广义线性混合模型（GLMM）执行Bootstrapping后，仅能得到参考水平的置信区间（CI），需要获取diet（二元因子）和Reproductive_output（三水平因子）所有水平（含交互组合）的置信区间，且结果布局类似ANOVA表。

模型及Bootstrap代码如下：

library("glmmTMB")
mod <- glmmTMB(body_phosphorus_content ~ diet*Reproductive_output, data = data, family = "Gamma")
summary(mod)
car::Anova(mod, type = c("III"))

# 残差检验
library("DHARMa")
sim <- simulateResiduals(mod)
plot(sim)
plotResiduals(sim, data$Reproductive_output)
plotResiduals(sim, data$diet)

# Bootstrap
library("car")
betahat.bootG <- Boot(mod, R=1000) 
summary(betahat.bootG)
confint(betahat.bootG)

解决方案 #

默认Boot()函数输出的是模型系数（以参考水平为基准的对比值），要获取各因子水平的绝对置信区间，需自定义统计量函数，提取对应水平的预测值或边际均值，再执行Bootstrap。

方法1：获取所有因子组合的预测值置信区间 #

适用于查看交互作用下每个因子组合的置信区间：

# 1. 创建包含所有因子水平组合的新数据框
new_data <- expand.grid(
  diet = unique(data$diet),
  Reproductive_output = unique(data$Reproductive_output)
)

# 2. 定义统计量函数：返回各组合的响应尺度预测值（Gamma模型默认log链接，需转换回原始尺度）
get_predictions <- function(model) {
  predict(model, newdata = new_data, type = "response")
}

# 3. 运行Bootstrap
boot_pred <- Boot(mod, R = 1000, statistic = get_predictions)

# 4. 计算并整理置信区间
pred_ci <- confint(boot_pred)
result_df <- cbind(new_data, pred_ci)
colnames(result_df) <- c("diet", "Reproductive_output", "95%_lower_CI", "95%_upper_CI")

# 查看结果
print(result_df)

方法2：获取因子主效应的边际均值置信区间 #

更接近ANOVA表的布局，展示每个因子单独的主效应置信区间（平均另一个因子的效应）：

library(emmeans)

# 1. 定义统计量函数：提取diet和Reproductive_output的边际均值（响应尺度）
get_emmeans <- function(model) {
  # 获取diet的边际均值
  diet_emm <- emmeans(model, ~ diet, type = "response")
  # 获取Reproductive_output的边际均值
  repro_emm <- emmeans(model, ~ Reproductive_output, type = "response")
  # 返回合并的均值向量
  c(as.data.frame(diet_emm)$emmean, as.data.frame(repro_emm)$emmean)
}

# 2. 运行Bootstrap
boot_emmeans <- Boot(mod, R = 1000, statistic = get_emmeans)

# 3. 计算并整理置信区间
emm_ci <- confint(boot_emmeans)
# 为CI命名对应效应
names(emm_ci) <- c(
  paste0("diet_", unique(data$diet)),
  paste0("Reproductive_output_", unique(data$Reproductive_output))
)

# 转换为易读的数据框
emm_result_df <- data.frame(
  Effect = names(emm_ci),
  95%_lower_CI = emm_ci[,1],
  95%_upper_CI = emm_ci[,2]
)

print(emm_result_df)

说明 #

• 两种方法均基于Bootstrap样本计算置信区间，默认采用百分位数法，可通过confint()的method参数调整（如method="bca"）。
• 若模型无交互项，两种方法的结果会更简洁，主效应边际均值直接对应各因子水平的估计值。

内容的提问来源于stack exchange，提问作者bribina