如何对多重插补数据集的GAM拟合结果进行可视化？

云上老周

2026-6-23

多重插补后GAM模型的合并绘图解决方案

我正在处理含缺失值的数据集，用R的mice包做多重插补，之后用mgcv包在插补后的数据集上拟合广义可加模型（GAM）。我想绘制所有插补数据集上的模型整体拟合结果，但不知道怎么合并各插补数据集的绘图结果。

现有尝试代码

单独绘制各插补模型的拟合图

library(mgcv)
library(mice)
library(ggplot2)

# 设置随机种子保证可复现
set.seed(123)

# 生成模拟数据集
n <- 200
data <- data.frame(
  outcome = rnorm(n),
  exposure = rnorm(n),
  covariate1 = sample(0:1, n, replace = TRUE),
  covariate2 = rnorm(n, 40, 2),
  covariate3 = rnorm(n, 30, 1),
  covariate4 = sample(0:1, n, replace = TRUE),
  covariate5 = sample(0:1, n, replace = TRUE)
)
# 人为添加缺失值
data[sample(1:n, 50), "exposure"] <- NA
data[sample(1:n, 50), "covariate1"] <- NA

# 执行多重插补
imp <- mice(data, m = 5, method = 'pmm', seed = 123)
# 在插补数据集上拟合GAM模型
model <- with(imp, gam(outcome ~ s(exposure) + covariate1 + covariate2 + covariate3 + covariate4 + covariate5, method = 'REML'))

# 单独绘制每个插补模型的拟合图
fitted_models <- model$analyses
par(mfrow = c(2, 3))  # 调整布局以容纳所有图
for (i in 1:length(fitted_models)) {
  plot(fitted_models[[i]], shade = TRUE, rug = TRUE, residuals = TRUE, pch = 1, cex = 1, main = paste("插补数据集", i))
}

上述代码只能生成每个插补数据集的单独拟合图，无法得到合并后的整体结果。

尝试合并模型但失败的代码

# pooled.fit <- pool(model)
# plot.gam(pooled.fit)

pool()函数返回的是合并后的模型参数汇总，不是可直接用plot.gam()绘制的GAM对象，因此无法直接绘图。

解决方案：合并插补模型的平滑项拟合结果并绘图

要得到所有插补数据集上的整体拟合结果，需要提取每个插补模型中平滑项的拟合值、置信区间，然后计算它们的均值和合并后的置信区间，最后用ggplot2绘制。

完整代码

library(mgcv)
library(mice)
library(ggplot2)
library(dplyr)

# 设置随机种子
set.seed(123)

# 生成模拟数据集（同之前）
n <- 200
data <- data.frame(
  outcome = rnorm(n),
  exposure = rnorm(n),
  covariate1 = sample(0:1, n, replace = TRUE),
  covariate2 = rnorm(n, 40, 2),
  covariate3 = rnorm(n, 30, 1),
  covariate4 = sample(0:1, n, replace = TRUE),
  covariate5 = sample(0:1, n, replace = TRUE)
)
data[sample(1:n, 50), "exposure"] <- NA
data[sample(1:n, 50), "covariate1"] <- NA

# 多重插补与模型拟合
imp <- mice(data, m = 5, method = 'pmm', seed = 123)
model <- with(imp, gam(outcome ~ s(exposure) + covariate1 + covariate2 + covariate3 + covariate4 + covariate5, method = 'REML'))

# 提取每个插补模型的平滑项结果
smooth_results <- lapply(model$analyses, function(gam_mod) {
  # 提取s(exposure)的平滑项拟合数据（不直接绘图）
  smooth_dat <- plot(gam_mod, select = 1, seWithMean = TRUE, plot = FALSE)
  # 转换为数据框并添加插补编号
  data.frame(
    exposure = smooth_dat$x,
    fit = smooth_dat$fit,
    se = smooth_dat$se,
    lower = smooth_dat$fit - 1.96 * smooth_dat$se,
    upper = smooth_dat$fit + 1.96 * smooth_dat$se,
    imputation = as.character(which(model$analyses == gam_mod))
  )
})

# 合并所有插补数据集的结果
combined_smooth <- bind_rows(smooth_results)

# 按Rubin规则计算合并后的均值和置信区间
pooled_smooth <- combined_smooth %>%
  group_by(exposure) %>%
  summarize(
    mean_fit = mean(fit),
    # 计算合并标准误：综合插补内和插补间变异
    within_var = mean(se^2),
    between_var = var(fit),
    m = length(model$analyses),
    pooled_se = sqrt(within_var + between_var + between_var/m),
    lower = mean_fit - 1.96 * pooled_se,
    upper = mean_fit + 1.96 * pooled_se
  )

# 绘制合并后的拟合图
ggplot(pooled_smooth, aes(x = exposure, y = mean_fit)) +
  # 可选：添加各插补模型的拟合曲线用于对比
  geom_line(data = combined_smooth, aes(group = imputation), color = "gray", alpha = 0.5) +
  # 添加合并后的均值拟合曲线
  geom_line(color = "red", linewidth = 1) +
  # 添加合并后的置信区间
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "red", alpha = 0.2) +
  # 添加原始数据的rug图展示分布
  geom_rug(data = data, aes(x = exposure), sides = "b", alpha = 0.5) +
  labs(title = "多重插补后GAM模型的合并拟合结果",
       x = "暴露变量",
       y = "平滑项拟合值") +
  theme_minimal()

代码说明

提取平滑项数据：用plot(gam_mod, plot = FALSE)获取每个GAM模型中s(exposure)的拟合值、标准误等数据，避免直接生成单独绘图。
合并结果：将所有插补模型的平滑项数据合并到一个数据框中，便于后续计算。
计算合并统计量：采用Rubin规则计算合并后的拟合均值和标准误，同时考虑插补内和插补间的变异，得到更可靠的置信区间。
可视化：用ggplot2绘制合并后的核心拟合曲线与置信区间，可选添加各插补模型的拟合曲线作为对比，同时用rug图展示原始数据的分布情况。

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

如何对多重插补数据集的GAM拟合结果进行可视化？

现有尝试代码 #

单独绘制各插补模型的拟合图 #

尝试合并模型但失败的代码 #

解决方案：合并插补模型的平滑项拟合结果并绘图 #

完整代码 #