R语言GAM模型基于type='terms'预测结果绘制平滑函数

报错核心原因 #

之前的绘图代码报错是因为predict(..., type='terms')返回的preds是包含拟合值矩阵、标准误矩阵的列表，不能直接作为ggplot的y轴映射变量；且返回的平滑项数值是中心化后的相对贡献值，需要和预测用的newdata按行绑定、补充置信区间计算后才能绘图。

实现步骤 #

1. 整理预测结果为绘图数据框 #

首先提取两个平滑项的拟合值、标准误，计算95%置信区间，和原预测数据集合并：

library(ggplot2)
library(dplyr)
library(tidyr)

# 提取预测值、标准误，计算置信区间
plot_data <- newdata %>%
  mutate(
    # 提取两个平滑项的拟合值
    smooth_timestep = preds$fit[, "s(timeStep)"],
    smooth_month = preds$fit[, "s(month)"],
    # 提取标准误
    se_timestep = preds$se.fit[, "s(timeStep)"],
    se_month = preds$se.fit[, "s(month)"],
    # 计算95%置信区间上下限
    timestep_lwr = smooth_timestep - 1.96 * se_timestep,
    timestep_upr = smooth_timestep + 1.96 * se_timestep,
    month_lwr = smooth_month - 1.96 * se_month,
    month_upr = smooth_month + 1.96 * se_month
  )

2. 分开展示两个平滑项（和内置plot(mod)效果一致） #

把宽格式数据转成长格式，用分面分别绘制长期趋势、月份季节趋势：

# 宽转长适配分面
plot_data_long <- plot_data %>%
  pivot_longer(
    cols = c(smooth_timestep, smooth_month),
    names_to = "smooth_term",
    values_to = "fit_val"
  ) %>%
  mutate(
    # 匹配每个平滑项对应的x轴、置信区间
    x_val = ifelse(smooth_term == "smooth_timestep", timeStep, month),
    lwr = ifelse(smooth_term == "smooth_timestep", timestep_lwr, month_lwr),
    upr = ifelse(smooth_term == "smooth_timestep", timestep_upr, month_upr)
  )

# 绘图
ggplot(plot_data_long, aes(x = x_val, y = fit_val)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "grey80", alpha = 0.7) +
  geom_line(color = "red", linewidth = 1) +
  facet_wrap(~smooth_term, scales = "free_x",
             labeller = as_labeller(c(
               smooth_timestep = "s(timeStep) 长期趋势",
               smooth_month = "s(month) 季节趋势"
             ))) +
  labs(x = "变量取值", y = "平滑项中心化贡献值") +
  theme_bw()

运行后效果和内置plot输出一致，带置信区间阴影。

3. 同图展示趋势+季节效应 #

两个平滑项x轴量纲不同，同图展示更合理的方式是叠加得到最终拟合值，和原始观测、长期趋势对比：

# 提取模型截距，terms类型预测值加截距才是响应变量尺度的预测结果
mod_intercept <- coef(mod)[1]
plot_data <- plot_data %>%
  mutate(
    trend_only = mod_intercept + smooth_timestep,
    fit_total = mod_intercept + smooth_timestep + smooth_month
  )

ggplot(plot_data, aes(x = timeStep)) +
  geom_point(aes(y = co2), alpha = 0.3, size = 0.7) +
  geom_line(aes(y = trend_only, color = "长期趋势"), linewidth = 1) +
  geom_line(aes(y = fit_total, color = "叠加季节效应拟合值"), linewidth = 0.8) +
  scale_color_manual(values = c("长期趋势" = "red", "叠加季节效应拟合值" = "blue")) +
  labs(x = "时间步", y = "CO2浓度", color = NULL) +
  theme_bw()

注意事项 #

• predict(..., type='terms')返回的所有平滑项结果均做了中心化处理（均值为0），单独表示该变量对预测结果的相对贡献，不能直接对应原始CO2浓度尺度，需要加模型截距、其他项的贡献才能得到实际预测值。
• 月份是循环平滑项（bs="cc"），单独绘制月份效应时可以把x轴转为1-12的有序因子，更贴合月份的循环属性。

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