使用lme4计算负二项GLMM预测置信区间的方法合理性验证

现有代码正确性评估 #

• 你当前的代码计算的是固定效应部分在对数链接尺度上的总体平均置信区间，仅纳入了固定效应估计的不确定性，没有包含随机效应方差、负二项分布离散参数的估计不确定性，最终得到的区间会偏窄，样本量偏小时偏差更明显。如果你需要计算给定随机效应水平的条件预测区间，还需要额外纳入随机效应的预测方差。
• 关于poly()变量的兼容性：你的写法在预测数据集和建模数据集一致时是适用的，model.matrix()会基于输入数据集生成和建模时一致的正交多项式基；如果是对新数据集做预测，不建议手动构造poly()公式，避免多项式基缩放不匹配的问题。

更合理的实现方案 #

方案1：参数bootstrap法（推荐） #

这是lme4框架下精度最高的方案，可以自动纳入所有模型参数的估计不确定性，不需要依赖正态分布假设：

library(lme4)

# 定义bootstrap预测函数，返回链接尺度的总体平均预测值
boot_pred <- function(mod) {
  mm <- model.matrix(~ f1 + f2 + poly(f3,2), mod@frame)
  return(mm %*% fixef(mod))
}

# 运行bootstrap，迭代次数可根据算力调整，推荐≥500次
set.seed(123)
boot_res <- bootMer(m.nb, FUN = boot_pred, nsim = 500)

# 计算95%分位数置信区间
dd$y_est <- apply(boot_res$t, 2, mean, na.rm = TRUE)
dd$CIlower <- apply(boot_res$t, 2, quantile, 0.025, na.rm = TRUE)
dd$CIupper <- apply(boot_res$t, 2, quantile, 0.975, na.rm = TRUE)

# 如需转换为原始计数响应尺度，对所有值做指数转换即可
dd$y_est_resp <- exp(dd$y_est)
dd$CIlower_resp <- exp(dd$CIlower)
dd$CIupper_resp <- exp(dd$CIupper)

方案2：快速正态近似法 #

如果样本量较大、bootstrap运行成本过高，可以用近似法直接计算标准误：

# 直接用predict方法获取标准误，re.form=NA表示计算总体平均预测值
pred_res <- predict(m.nb, newdata = dd, re.form = NA, se.fit = TRUE)
dd$y_est <- pred_res$fit
dd$CIlower <- pred_res$fit - 1.96 * pred_res$se.fit
dd$CIupper <- pred_res$fit + 1.96 * pred_res$se.fit

如果需要计算给定随机效应g的条件预测区间，去掉re.form = NA参数即可，此时计算的区间会自动纳入随机效应的预测方差。

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