在mgcv中生成含随机效应的各水平平滑曲线用于绘图

问题背景 #

我们使用mgcv包构建了如下GAM模型，包含一个针对x1的平滑项和针对x2的随机截距项：

gam(y ~ s(x1) + s(x2, bs="re"), method = 'REML', 
    family = 'gaussian', data = my.data) 

通过gratia::smooth_estimates()可以获取x1的均匀间隔平滑预测值，以及x2各水平的随机效应，结果存储在同一数据框中。需要找到简洁方法，生成x2每个水平下x1的完整平滑估计曲线（即把每个x1的平滑值与对应x2的随机效应相加，得到x2各水平对应的可堆叠绘制的曲线）。

简洁实现方法 #

方法一：基于smooth_estimates拆分合并 #

利用dplyr和tidyr的工具拆分平滑项与随机效应，再交叉组合计算总估计值：

library(dplyr)
library(tidyr)

# 获取平滑估计结果
smooth_res <- smooth_estimates(your_gam_model)

# 拆分x1的平滑值和x2的随机效应
x1_smooth_vals <- smooth_res %>%
  filter(smooth == "s(x1)") %>%
  select(x1, smooth_est = est)

x2_re_vals <- smooth_res %>%
  filter(smooth == "s(x2)") %>%
  select(x2, re_est = est)

# 交叉组合并计算每个x2水平下的总平滑估计
grouped_smooth_curves <- x1_smooth_vals %>%
  crossing(x2_re_vals) %>%
  mutate(total_estimate = smooth_est + re_est)

方法二：直接用fitted_values生成全组合预测 #

更高效的方式是用gratia::fitted_values()直接生成包含所有x1和x2组合的拟合值，无需手动拆分相加：

# 生成包含x1全范围和x2所有水平的新数据集
new_pred_data <- expand.grid(
  x1 = seq(min(my.data$x1), max(my.data$x1), length.out = 100), # 自定义x1的间隔数
  x2 = unique(my.data$x2)
)

# 获取所有组合的拟合值（包含平滑项+随机效应）
grouped_fitted_vals <- fitted_values(your_gam_model, data = new_pred_data)

两种方法得到的结果都可以直接用于可视化，比如用ggplot2绘制堆叠曲线：

library(ggplot2)
ggplot(grouped_fitted_vals, aes(x = x1, y = fitted, color = x2)) +
  geom_line(linewidth = 0.8)

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