R语言绘制多重插补数据下逻辑回归的显著交互项

针对你这个用MICE多重插补后绘制逻辑回归交互效应的需求，我整理了一套实操性强的R代码方案，专门针对你提到的「4分类自变量（CATIV）参考类与第3类和连续调节变量（MOD）的显著交互」场景：

第一步：处理多重插补的模型结果 #

你已经用with(iX, glm(...))拟合了多个插补数据集的逻辑回归，首先可以先汇总模型结果确认交互项的显著性：

# 加载必要工具包
library(mice)
library(ggplot2)

# 汇总多重插补的模型结果
pooled_fit <- pool(fit)
# 查看汇总后的系数，重点确认参考类与CATIV第3类的交互项p值
summary(pooled_fit)

第二步：构建预测数据集 #

我们只需要聚焦交互显著的两个类别（参考类和CATIV第3类），同时生成连续调节变量MOD的合理取值范围（用均值±2倍标准差覆盖绝大多数数据）：

# 获取MOD的描述性统计（忽略缺失值）
mod_mean <- mean(iX$MOD, na.rm = TRUE)
mod_sd <- sd(iX$MOD, na.rm = TRUE)

# 创建预测数据框：仅保留交互显著的两个类别
# 替换<参考类名称>和<CATIV3名称>为你实际的类别标签（比如"Control"和"Group3"）
pred_data <- expand.grid(
  CATIV = factor(c("<参考类名称>", "<CATIV3名称>"), levels = levels(iX$CATIV)),
  MOD = seq(mod_mean - 2*mod_sd, mod_mean + 2*mod_sd, length.out = 100)
)

第三步：获取带置信区间的预测概率 #

这里提供两种严谨的方式，任选其一即可：

方式1：用汇总后的pooled模型预测 #

# 获取预测概率及标准误
preds <- predict(pooled_fit, newdata = pred_data, type = "response", se.fit = TRUE)

# 计算95%置信区间并合并到数据集
pred_data$pred_prob <- preds$fit
pred_data$lower_ci <- preds$fit - 1.96 * preds$se.fit
pred_data$upper_ci <- preds$fit + 1.96 * preds$se.fit

方式2：直接用mids对象预测（保留插补间变异） #

这种方式会对每个插补数据集单独预测，再取均值和插补间标准差，更贴合多重插补的统计逻辑：

# 对所有插补数据集做预测
preds_mids <- predict(fit, newdata = pred_data, type = "response", se.fit = TRUE)

# 计算所有插补结果的均值（最终预测概率）
pred_data$pred_prob <- rowMeans(preds_mids$fit)
# 计算插补间的标准误
pred_data$se <- apply(preds_mids$fit, 1, sd)
# 计算95%置信区间
pred_data$lower_ci <- pred_data$pred_prob - 1.96 * pred_data$se
pred_data$upper_ci <- pred_data$pred_prob + 1.96 * pred_data$se

第四步：绘制交互效应图（ggplot2） #

我们绘制预测概率随MOD变化的折线图，加上置信区间，直观展示两个类别的差异变化：

ggplot(pred_data, aes(x = MOD, y = pred_prob)) +
  # 绘制置信区间阴影
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci, fill = CATIV), alpha = 0.2, color = NA) +
  # 绘制交互折线
  geom_line(aes(color = CATIV), linewidth = 1.2) +
  # 自定义标签和主题
  labs(
    x = "调节变量(MOD)",
    y = "因变量(BINARYDV=1)的预测概率",
    color = "自变量(CATIV)",
    fill = "自变量(CATIV)",
    title = "CATIV参考类与第3类的交互效应",
    subtitle = "基于多重插补逻辑回归模型的预测结果"
  ) +
  # 可选：添加显著性标注（替换x/y坐标为你图中交互明显的位置）
  annotate("text", x = mod_mean + 1.5*mod_sd, y = 0.7, 
           label = "交互效应显著(p<0.05)", color = "darkred", fontface = "bold") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

关键注意事项 #

• 确保pred_data中的CATIV因子水平与原数据完全一致，否则预测会出错
• 如果需要展示所有4个类别的交互，只需把pred_data中的CATIV改成所有4个水平即可，但聚焦显著的两组会让图表更清晰
• 若想展示logit刻度的交互，可把type="response"改成type="link"，但概率刻度更易被非统计背景的读者理解

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