如何在ggplot中绘制多个GAMM平滑曲线？含fs项与全局项合并

一、当前因子平滑项（s(CYR.std, fSite, bs = "fs")）的绘制是否准确？ #

当前的绘制方式并不准确。你用fitted_values()得到的是整个模型的完整拟合值，包含了s(sal)、季节效应、全局s(CYR.std)等所有项的贡献，并非单独的站点偏差因子平滑项。

fs类型的平滑项本质是每个站点相对于全局CYR.std趋势的特异性偏差，如果要单独展示这部分偏差，应该直接提取该平滑项的估计值，而非全模型拟合值：

# 提取因子平滑项的偏差估计（响应尺度）
fs_smooth <- smooth_estimates(m, smooth = "s(CYR.std,fSite)", scale = "response")

# 绘制仅站点偏差的曲线
ggplot(fs_smooth, aes(x = CYR.std, y = est, color = fSite)) +
  geom_line(alpha = 0.6) +
  labs(title = "站点特异性时间偏差（相对于全局趋势）")

如果你的目标是展示包含全局效应的每个站点完整拟合曲线，那当前fitted_values()的做法是可行的——此时曲线已经叠加了全局趋势和其他协变量的固定效应（默认会把其他协变量设为均值/参考水平）。

二、如何将全局项与因子平滑项合并到同一图中？ #

分两种场景处理：

场景1：展示全局时间趋势 + 各站点的偏差曲线 #

先提取全局s(CYR.std)的估计值，再和站点偏差曲线叠加，能清晰对比全局趋势和站点特异性差异：

# 提取全局CYR.std平滑项的估计（含置信区间）
global_smooth <- smooth_estimates(m, smooth = "s(CYR.std)", scale = "response")

# 合并绘图
ggplot() +
  # 绘制所有站点的偏差曲线（降低透明度避免杂乱）
  geom_line(data = fs_smooth, aes(x = CYR.std, y = est, color = fSite), alpha = 0.3) +
  # 绘制全局趋势曲线（加粗突出）
  geom_line(data = global_smooth, aes(x = CYR.std, y = est), color = "black", linewidth = 1.2) +
  # 全局趋势的置信区间
  geom_ribbon(data = global_smooth, aes(x = CYR.std, ymin = lower, ymax = upper), alpha = 0.2, fill = "black") +
  labs(x = "标准化年份", y = "响应值", title = "全局时间趋势与站点偏差对比") +
  theme_minimal()

场景2：展示全局时间趋势 + 各站点完整拟合曲线 #

如果想展示每个站点的完整拟合（已包含全局趋势），同时叠加单独的全局曲线，直接合并两个数据集即可：

# 给全局拟合数据添加分组标识，方便区分
fv2$fSite <- "全局效应"

# 合并站点拟合和全局拟合的数据
combined_fv <- rbind(fv, fv2)

# 绘图
ggplot(combined_fv, aes(x = CYR.std, y = fitted, color = fSite)) +
  # 绘制所有站点的完整拟合曲线（降低透明度）
  geom_line(data = subset(combined_fv, fSite != "全局效应"), alpha = 0.4) +
  # 绘制全局效应曲线（加粗突出）
  geom_line(data = subset(combined_fv, fSite == "全局效应"), linewidth = 1.2) +
  # 添加全局效应的置信区间
  geom_ribbon(data = fv2, aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "black", color = NA) +
  labs(x = "标准化年份", y = "预测计数", title = "全局效应与各站点完整拟合对比") +
  theme_minimal() +
  # 调整图例，让全局效应的线条不透明
  guides(color = guide_legend(override.aes = list(alpha = 1)))

补充提示 #

• gratia包的smooth_estimates()可以精准提取单个平滑项的估计，避免其他协变量的干扰，是单独展示特定效应的最佳方式。
• 如果模型包含by交互项（比如你代码中的季节交互），提取平滑项时可以通过terms参数指定，确保拿到的是你需要的部分。

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