如何获取R中survey包svyglm复杂调查模型预测值的95%置信区间？

prediction包的prediction_summary()函数对survey包生成的复杂调查模型（如svyglm对象）的置信区间计算支持有限，因此直接调用无法返回预期的95%CI。可以通过以下两种方法解决：

方法一：使用survey包原生的predict.svyglm函数 #

survey包自带的predict()方法专门针对复杂调查模型设计，能正确计算预测值的标准误，再手动推导置信区间：

library(survey)
data(api)

# 构建二分类结果变量
apistrat$api00_bin <- ifelse(apistrat$api00 < mean(apistrat$api00), 1, 0)

# 构建调查设计
dstrat <- svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, fpc=~fpc)

# 拟合svyglm模型
api00bin.reg <- svyglm(api00_bin~enroll+avg.ed+yr.rnd+sch.wide, design=dstrat, family = poisson(link = "log"))

# 生成预测数据：固定其他变量为均值，仅改变sch.wide
new_data <- expand.grid(
  enroll = mean(apistrat$enroll),
  avg.ed = mean(apistrat$avg.ed),
  yr.rnd = unique(apistrat$yr.rnd)[1], # 取第一个水平或均值，可按需调整
  sch.wide = c("No", "Yes"),
  stype = unique(apistrat$stype)[1] # 调查设计需保留分层变量，取默认水平
)

# 获取预测值和标准误
preds <- predict(api00bin.reg, newdata = new_data, type = "response", se.fit = TRUE)

# 计算95%置信区间（使用t分布临界值，自由度取模型残差自由度）
crit_val <- qt(0.975, df = api00bin.reg$df.residual)
lower_ci <- preds$fit - crit_val * preds$se.fit
upper_ci <- preds$fit + crit_val * preds$se.fit

# 整理结果
result <- data.frame(
  sch.wide = new_data$sch.wide,
  predicted_value = preds$fit,
  lower_95ci = lower_ci,
  upper_95ci = upper_ci
)

print(result)

方法二：基于prediction包结果手动计算置信区间 #

如果坚持使用prediction_summary()，可以提取其返回的标准误，再手动计算置信区间：

library(prediction)

# 获取预测摘要
pred_sum <- prediction::prediction_summary(api00bin.reg, at = list(sch.wide = c("No", "Yes")), type = "response", level = 0.95, calculate_se = TRUE)

# 提取标准误并计算95%CI
crit_val <- qt(0.975, df = api00bin.reg$df.residual)
pred_sum$lower_95ci <- pred_sum$prediction - crit_val * pred_sum$se
pred_sum$upper_95ci <- pred_sum$prediction + crit_val * pred_sum$se

# 查看包含CI的结果
print(pred_sum)

注意事项 #

• 对于复杂调查模型，优先使用survey包原生方法，因为它会考虑调查设计的权重、分层等因素，计算的标准误更准确。
• 若模型是二分类结果，通常更适合使用family = binomial(link = "logit")而非泊松族，可根据实际需求调整模型族。

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