如何在同一张图中绘制多个GAM模型的平滑曲线

要解决将两个bam模型的平滑曲线合并到同一张ggplot2图中，同时修正y轴刻度的问题，可以按以下步骤操作：

1. 加载必要的R包 #

确保已安装所需包，然后加载：

library(mgcv)
library(gratia)
library(tidyverse)

2. 提取两个模型的平滑曲线预测数据 #

使用gratia::smooth_estimates()分别提取两个模型中s(time)的平滑估计值，添加标识列区分模型：

# 提取TOTAL_S模型的平滑数据
smooth_s <- smooth_estimates(bam_s, smooth = "s(time)") %>%
  mutate(model = "TOTAL_S")

# 提取TOTAL_P模型的平滑数据
smooth_p <- smooth_estimates(bam_p, smooth = "s(time)") %>%
  mutate(model = "TOTAL_P")

# 合并两个数据集
combined_smooths <- bind_rows(smooth_s, smooth_p)

3. 用ggplot2绘制合并后的平滑曲线 #

绘制时指定x轴为时间，y轴为平滑估计值，用颜色区分模型，添加置信区间，并将y轴设置为整数刻度（匹配药物数量的整数属性）：

ggplot(combined_smooths, aes(x = time, y = estimate, color = model)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = model), alpha = 0.2, color = NA) +
  scale_y_continuous(breaks = seq(floor(min(combined_smooths$lower)), ceiling(max(combined_smooths$upper)), by = 1)) +
  labs(x = "事件前年份", y = "药物数量", color = "药物类型", fill = "药物类型") +
  theme_minimal()

关键说明 #

• smooth_estimates()会自动计算平滑曲线的95%置信区间，若需调整置信水平，可添加level = 0.90这类参数
• y轴刻度通过scale_y_continuous()结合seq()生成整数刻度，解决了刻度不符合药物数量属性的问题

完整可运行代码 #

# 加载包
library(mgcv)
library(gratia)
library(tidyverse)

# 生成模拟数据
dset <- data.frame(studynr=rep(c(1:20), each =10),
                   time = rep(c(-9,-8,-7,-6,-5,-4,-3,-2,1,0), times=20), 
                   TOTAL_P = sample(1:10, 200, replace=T),
                   TOTAL_S = sample(1:10, 200, replace=T))

# 构建bam模型
bam_s <- bam(TOTAL_S ~ s(time)
             + s(studynr, bs = 're') 
             + s(studynr, time, bs = 're'),
             data = dset, method = "REML")

bam_p <- bam(TOTAL_P ~ s(time)
             + s(studynr, bs = 're')
             + s(studynr, time, bs = 're'),
             data = dset, method = "REML")

# 提取并合并平滑数据
smooth_s <- smooth_estimates(bam_s, smooth = "s(time)") %>% mutate(model = "TOTAL_S")
smooth_p <- smooth_estimates(bam_p, smooth = "s(time)") %>% mutate(model = "TOTAL_P")
combined_smooths <- bind_rows(smooth_s, smooth_p)

# 绘图
ggplot(combined_smooths, aes(x = time, y = estimate, color = model)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = model), alpha = 0.2, color = NA) +
  scale_y_continuous(breaks = seq(floor(min(combined_smooths$lower)), ceiling(max(combined_smooths$upper)), by = 1)) +
  labs(x = "事件前年份", y = "药物数量", color = "药物类型", fill = "药物类型") +
  theme_minimal()

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