R中GLMM获取随机截距水平的边际效应置信区间方法问询

现象原因说明 #

1. 预测区间过宽：merTools的predictInterval默认同时纳入了固定效应估计不确定性、随机效应估计不确定性、以及结局层面的残差变异（二分类结局对应伯努利分布的随机误差），相当于给虚拟患者的个体预测区间，而非站点效应本身的区间，所以远宽于预期的仅反映站点绩效差异的区间。
2. 多次运行结果波动：predictInterval返回的点估计是采样模拟结果的均值，而非原始模型固定效应+随机效应BLUP的固定组合，随机采样过程会导致每次结果存在微小差异，设置随机种子即可解决复现问题。

可行方案 #

方案1：快速正态近似法（秒出结果，适合临床发表场景） #

忽略方差参数的估计不确定性（样本量3000情况下偏差极小），直接提取随机效应的条件标准误计算OR的置信区间：

# 提取站点水平随机效应及条件方差
re <- ranef(mod, condVar = TRUE)
site_re <- re[["AmbulanceZone:AmbulanceStation"]]
# 提取每个站点随机效应的条件标准误
site_re$se <- attr(re[["AmbulanceZone:AmbulanceStation"]], "postVar")[1, 1, ]

# 计算站点相对于平均水平的OR及95%置信区间
site_re$or <- exp(site_re[["(Intercept)"]])
site_re$or_lower <- exp(site_re[["(Intercept)"]] - 1.96 * site_re$se)
site_re$or_upper <- exp(site_re[["(Intercept)"]] + 1.96 * site_re$se)

如果需要虚拟患者的预测概率区间，直接把固定效应的线性部分加上随机效应再做logit转换即可，点估计完全固定无波动。

方案2：参数模拟法（更准确，速度远快于全模型bootstrap） #

仅对模型系数的联合分布采样，无需每次重拟合模型，兼顾固定效应和随机效应的不确定性，同时保证结果可复现：

library(arm)
# 设置随机种子保证结果固定
set.seed(123)
# 采样1000次模型系数（含固定、随机效应）
sims <- sim(mod, n.sims = 1000)

# 手动计算虚拟患者的固定效应线性部分
fix_val <- fixef(mod)[["(Intercept)"]] + fixef(mod)[["genderMale"]] + 
  fixef(mod)[["logresponsetime"]] * 2.0794 + fixef(mod)[["rhythmShockable"]]
# 固定点估计（无波动）
site_re$pred_prob <- plogis(fix_val + site_re[["(Intercept)"]])

# 计算模拟后的95%预测区间
site_sim <- t(apply(sims@ranef[["AmbulanceZone:AmbulanceStation"]], 2, function(x) {
  fix_sim <- fixef(sims)[,"(Intercept)"] + fixef(sims)[,"genderMale"] + 
    fixef(sims)[,"logresponsetime"] * 2.0794 + fixef(sims)[,"rhythmShockable"]
  plogis(fix_sim + x)
}))
site_re$pred_lower <- apply(site_sim, 1, quantile, 0.025)
site_re$pred_upper <- apply(site_sim, 1, quantile, 0.975)

如果不需要纳入残差变异，关掉对应参数后区间宽度完全符合站点效应的预期范围。

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