如何使用ggplot绘制含月份项GAM模型的平滑函数

你当前已经完成模型拟合与预测计算，不需要依赖geom_smooth()重复拟合，直接基于预测结果整理数据即可绘制和模型完全匹配的平滑曲线，避免geom_smooth()内部参数设置和你现有模型不一致导致的结果偏差。

第一步：先把预测结果整合到原数据框 #

你之前用predict(type='terms')输出的结果已经包含了每个平滑项的独立拟合值、标准误，直接和原数据按列绑定即可，不用额外计算：

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

carbonD <- carbonD %>%
  mutate(
    # 整体拟合值（加截距）
    fit_all = preds$fit[,1] + preds$fit[,2] + coef(mod)[1],
    # 长期趋势项s(timeStep)的拟合值与95%置信区间
    fit_timestep = preds$fit[,1],
    se_timestep = preds$se.fit[,1],
    lwr_timestep = fit_timestep - 1.96*se_timestep,
    upr_timestep = fit_timestep + 1.96*se_timestep,
    # 月度周期项s(month)的拟合值与95%置信区间
    fit_month = preds$fit[,2],
    se_month = preds$se.fit[,2],
    lwr_month = fit_month - 1.96*se_month,
    upr_month = fit_month + 1.96*se_month
  )

第二步：绘制单个平滑项的独立效应 #

1. 长期时间趋势平滑项 #

直接以timeStep为x轴，叠加置信区间带和平滑线即可：

ggplot(carbonD, aes(x = timeStep)) +
  geom_ribbon(aes(ymin = lwr_timestep + coef(mod)[1], 
                  ymax = upr_timestep + coef(mod)[1]),
              alpha = 0.2, fill = "#2c7fb8") +
  geom_line(aes(y = fit_timestep + coef(mod)[1]), 
            color = "#2c7fb8", linewidth = 1) +
  labs(x = "观测时间步", y = "CO2浓度（校正月度效应）",
       title = "GAM长期趋势平滑项") +
  theme_bw()

2. 月度周期平滑项 #

因为你给month设置的是循环样条(bs="cc")，需要补一个端点值让曲线闭合，避免12月和1月之间出现断档：

# 提取月度项的唯一拟合值，补充端点
month_effect <- carbonD %>%
  distinct(month, fit_month, lwr_month, upr_month) %>%
  add_row(month = 13, 
          fit_month = .$fit_month[.$month == 1],
          lwr_month = .$lwr_month[.$month == 1],
          upr_month = .$upr_month[.$month == 1])

ggplot(month_effect, aes(x = month)) +
  geom_ribbon(aes(ymin = lwr_month, ymax = upr_month),
              alpha = 0.2, fill = "#d95f0e") +
  geom_line(aes(y = fit_month), color = "#d95f0e", linewidth = 1) +
  scale_x_continuous(breaks = 1:12, limits = c(1,13)) +
  labs(x = "月份", y = "CO2浓度周期偏离值（校正长期趋势）",
       title = "GAM月度周期平滑项") +
  theme_bw()

可选：绘制模型整体拟合效果 #

如果需要对照原始观测值展示模型拟合情况，直接用整合好的整体拟合值绘图即可：

ggplot(carbonD, aes(x = timeStep)) +
  geom_point(aes(y = co2), alpha = 0.3, size = 1.2) +
  geom_line(aes(y = fit_all), color = "#de2d26", linewidth = 0.8) +
  labs(x = "观测时间步", y = "CO2浓度",
       title = "GAM模型整体拟合效果") +
  theme_bw()

关于geom_smooth()的说明 #

不推荐直接用geom_smooth()绘制你现有模型的平滑项：该函数需要在绘图内部重新调用模型拟合流程，你需要手动重复一遍模型公式、样条类型、k值、分布族等所有参数，一旦漏设或错配参数，绘出的曲线和你实际拟合的mod结果会出现偏差，远不如直接用已经跑完的预测结果绘图准确。

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