GLMER混合逻辑回归模型预测值的置信区间求解问题

针对带有随机效应的混合广义逻辑回归（glmer）模型，普通glm的置信区间计算方法不适用，以下是几种可靠的解决方案，同时解决你遇到的字符型随机因子失效、preds$fit缺失的问题：

方法1：用emmeans计算边际置信区间 #

适合按分组（如火灾强度）计算总体平均预测值的置信区间，自动处理字符型随机因子，支持logit到概率尺度的转换。

# 加载依赖包
library(lme4)
library(emmeans)
library(dplyr)

# 模拟匹配你场景的可复现数据
set.seed(123)
n <- 500
dat <- data.frame(
  plot = sample(paste0("plot_", 1:20), n, replace = TRUE), # 字符型随机因子
  fire_intensity = sample(c("低", "中", "高"), n, replace = TRUE),
  fire_height = rnorm(n, 10, 3),
  tree_diam = rnorm(n, 20, 5),
  mortality = rbinom(n, 1, plogis(-5 + 0.2*fire_height + 0.8*as.numeric(fire_intensity) + 0.1*tree_diam + rnorm(20, 0.5)[match(dat$plot, unique(dat$plot))]))
)

# 拟合glmer模型
model <- glmer(mortality ~ fire_height + fire_intensity + tree_diam + (1|plot), 
               data = dat, family = binomial)

# 按火灾强度分组，固定其他变量为均值计算预测值
emm <- emmeans(model, ~ fire_intensity, 
               at = list(fire_height = mean(dat$fire_height), 
                         tree_diam = mean(dat$tree_diam)))

# 转换为概率尺度的预测值及95%置信区间
emm_probs <- as.data.frame(regrid(emm, transform = "response"))
print(emm_probs)

输出结果包含emmean（预测值）、lower.CL和upper.CL（置信区间上下限），直接对应死亡率的概率估计。

方法2：用merTools::predictInterval计算置信/预测区间 #

通过模拟方法生成区间，可灵活选择是否包含随机效应（置信区间对应固定效应不确定性，预测区间包含随机效应波动）。

library(merTools)

# 构造预测数据集（按火灾强度分组，固定协变量为均值）
pred_dat <- expand.grid(
  fire_intensity = c("低", "中", "高"),
  fire_height = mean(dat$fire_height),
  tree_diam = mean(dat$tree_diam),
  plot = unique(dat$plot)[1] # 仅需提供任意一个随机因子水平，不影响置信区间计算
)

# 计算95%置信区间（不包含随机效应）
preds <- predictInterval(model, newdata = pred_dat, 
                         level = 0.95, type = "prob", 
                         include.random = FALSE, n.sims = 1000)

# 整理结果
pred_result <- cbind(pred_dat[, c("fire_intensity")], preds)
print(pred_result)

结果中的fit是预测值，lwr/upr是置信区间上下限；若需包含随机效应的预测区间，设置include.random = TRUE即可。

方法3：自助法（Bootstrapping）自定义计算 #

适合复杂场景，通过重抽样聚类单元（如地块）生成预测值分布，计算分位数作为置信区间。

library(purrr)

# 定义自助抽样函数
boot_pred <- function(dat, formula) {
  # 按地块重抽样（混合模型需聚类抽样）
  sampled_plots <- sample(unique(dat$plot), replace = TRUE)
  boot_dat <- dat %>% filter(plot %in% sampled_plots) %>% 
    group_by(plot) %>% sample_frac(replace = TRUE) %>% ungroup()
  # 尝试拟合模型（避免抽样失败）
  boot_model <- tryCatch(glmer(formula, data = boot_dat, family = binomial),
                         error = function(e) NULL)
  if(is.null(boot_model)) return(NA)
  # 生成预测
  pred_dat <- expand.grid(
    fire_intensity = c("低", "中", "高"),
    fire_height = mean(dat$fire_height),
    tree_diam = mean(dat$tree_diam),
    plot = unique(dat$plot)[1]
  )
  predict(boot_model, newdata = pred_dat, type = "response")
}

# 运行1000次自助抽样
set.seed(123)
boot_results <- replicate(1000, boot_pred(dat, mortality ~ fire_height + fire_intensity + tree_diam + (1|plot)))

# 计算95%置信区间
boot_ci <- apply(boot_results, 1, function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE))
# 整理结果
boot_result_df <- data.frame(
  fire_intensity = c("低", "中", "高"),
  fit = colMeans(boot_results, na.rm = TRUE),
  lwr = boot_ci[1,],
  upr = boot_ci[2,]
)
print(boot_result_df)

问题溯源与解决 #

• preds$fit缺失：通常是误用了glm的predict函数或未指定type="response"，以上方法均明确返回带fit（或对应列）的结果。
• REquantile失效：该函数对随机因子类型要求严格，将字符型转为因子型（dat$plot <- as.factor(dat$plot)）可尝试修复，但上述方法无需转换即可直接支持字符型随机因子。

内容的提问来源于stack exchange，提问作者Camille Revertégat