如何在同一张图中绘制多个GAM模型的参数化连续变量？

一、同一参数化变量（如两个模型中的z）同图绘制 #

方法1：使用mgcViz提取绘图数据后合并 #

mgcViz的pterm()函数可以提取参数化项的绘图数据，我们可以将两个模型的数据合并后用ggplot统一绘制：

library(mgcViz)
library(ggplot2)
library(dplyr)

# 将模型转换为gamViz对象
gv1 <- getViz(g1)
gv2 <- getViz(g2)

# 提取z变量的偏效应数据，并标记模型来源
df_g1 <- pterm(gv1, "z")$data %>% mutate(model = "g1")
df_g2 <- pterm(gv2, "z")$data %>% mutate(model = "g2")

# 合并数据并绘制对比图
bind_rows(df_g1, df_g2) %>%
  ggplot(aes(x = x, y = y, color = model)) +
  geom_line(linewidth = 1) +
  # 添加95%置信区间（2倍标准误）
  geom_ribbon(aes(ymin = y - 2*se, ymax = y + 2*se, fill = model), 
              alpha = 0.2, color = NA) +
  labs(x = "z", y = "Partial Effect", color = "Model", fill = "Model") +
  theme_minimal()

方法2：手动提取偏效应（无需mgcViz） #

如果不用mgcViz，可以通过predict.gam()的type="terms"参数直接提取参数化项的偏效应，再整理数据绘图：

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

# 生成覆盖两个模型z变量范围的序列
z_range <- range(c(g1$model$z, g2$model$z))
z_seq <- seq(z_range[1], z_range[2], length.out = 100)

# 提取g1中z的偏效应及标准误
pred_g1 <- predict(g1, 
                   newdata = data.frame(z = z_seq, y = mean(g1$model$y)),  # 控制其他变量为均值
                   type = "terms", 
                   se.fit = TRUE)
df_g1 <- data.frame(
  z = z_seq,
  effect = pred_g1$fit[, "z"],
  se = pred_g1$se.fit[, "z"],
  model = "g1"
)

# 提取g2中z的偏效应及标准误
pred_g2 <- predict(g2, 
                   newdata = data.frame(z = z_seq, y1 = mean(g2$model$y1)),
                   type = "terms", 
                   se.fit = TRUE)
df_g2 <- data.frame(
  z = z_seq,
  effect = pred_g2$fit[, "z"],
  se = pred_g2$se.fit[, "z"],
  model = "g2"
)

# 合并绘图
bind_rows(df_g1, df_g2) %>%
  ggplot(aes(x = z, y = effect, color = model)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = effect - 2*se, ymax = effect + 2*se, fill = model),
              alpha = 0.2, color = NA) +
  labs(x = "z", y = "Partial Effect", color = "Model", fill = "Model") +
  theme_minimal()

二、不同参数化变量（如z和z1）同图绘制 #

可以实现，但需要注意变量的尺度和含义是否适合直接对比。如果变量单位或取值范围差异大，建议先标准化变量，或用分面展示：

方案1：标准化变量后同图对比 #

将z和z1转换为z-score（标准化），让偏效应的尺度具有可比性：

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

# 处理g3的z变量
z_mean <- mean(g3$model$z)
z_sd <- sd(g3$model$z)
z_seq_std <- seq(min((g3$model$z - z_mean)/z_sd), max((g3$model$z - z_mean)/z_sd), length.out = 100)
z_seq_original <- z_seq_std * z_sd + z_mean  # 保留原始值用于预测

pred_g3 <- predict(g3, 
                   newdata = data.frame(z = z_seq_original, y = mean(g3$model$y)),
                   type = "terms", 
                   se.fit = TRUE)
df_g3 <- data.frame(
  variable = "z",
  x_std = z_seq_std,
  effect = pred_g3$fit[, "z"],
  se = pred_g3$se.fit[, "z"]
)

# 处理g4的z1变量
z1_mean <- mean(g4$model$z1)
z1_sd <- sd(g4$model$z1)
z1_seq_std <- seq(min((g4$model$z1 - z1_mean)/z1_sd), max((g4$model$z1 - z1_mean)/z1_sd), length.out = 100)
z1_seq_original <- z1_seq_std * z1_sd + z1_mean

pred_g4 <- predict(g4, 
                   newdata = data.frame(z1 = z1_seq_original, y1 = mean(g4$model$y1)),
                   type = "terms", 
                   se.fit = TRUE)
df_g4 <- data.frame(
  variable = "z1",
  x_std = z1_seq_std,
  effect = pred_g4$fit[, "z1"],
  se = pred_g4$se.fit[, "z1"]
)

# 合并绘图
bind_rows(df_g3, df_g4) %>%
  ggplot(aes(x = x_std, y = effect, color = variable)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = effect - 2*se, ymax = effect + 2*se, fill = variable),
              alpha = 0.2, color = NA) +
  labs(x = "Standardized Variable (z-score)", y = "Partial Effect", 
       color = "Variable", fill = "Variable") +
  theme_minimal()

方案2：分面展示（避免尺度混淆） #

如果不想标准化，用分面可以清晰展示各自的偏效应，同时保留原始变量尺度：

# 基于上面的df_g3和df_g4，调整数据格式
df_g3 <- df_g3 %>% mutate(x_original = z_seq_original)
df_g4 <- df_g4 %>% mutate(x_original = z1_seq_original)

bind_rows(df_g3, df_g4) %>%
  ggplot(aes(x = x_original, y = effect)) +
  geom_line(color = "#2c3e50", linewidth = 1) +
  geom_ribbon(aes(ymin = effect - 2*se, ymax = effect + 2*se),
              alpha = 0.2, fill = "#2c3e50") +
  facet_wrap(~variable, scales = "free_x") +  # 每个变量用独立的x轴
  labs(x = "Variable Value", y = "Partial Effect") +
  theme_minimal()

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