如何在GAM模型中获取不同X1取值的预测值并绘制对应优势比图

步骤1：构造标准化预测数据集 #

要单独提取x1的效应，需要将模型中其他协变量（这里是x0）固定为恒定值（通常用样本均值/中位数，消除其他变量的影响），同时生成覆盖x1实际取值范围的均匀序列：

# 加载所需包
library(mgcv)
library(ggplot2)

# 你的原始建模代码
dat <- mgcv::gamSim(1,n=400,dist="binary",scale=2)
adjY <- gam(y~ s(x0)+s(x1),family=binomial(link="logit"),data=dat)

# 构造预测数据集
new_dat <- data.frame(
  x1 = seq(min(dat$x1), max(dat$x1), length.out = 200), # 覆盖x1全部取值的200个点，绘图更平滑
  x0 = mean(dat$x0) # x0固定为样本均值，控制其影响
)

步骤2：计算预测值及优势比 #

在logit链接尺度下做预测，两个不同x1取值对应的预测值差值的指数就是二者的优势比，我们以x1=0为参照组：

# 预测对数优势（link尺度）及标准误
pred_link <- predict(adjY, newdata = new_dat, type = "link", se.fit = TRUE)

# 找到x1=0对应的预测值作为参照
ref_row <- which.min(abs(new_dat$x1 - 0))
ref_logit <- pred_link$fit[ref_row]

# 计算各x1对应的优势比及95%置信区间
new_dat$or <- exp(pred_link$fit - ref_logit)
new_dat$or_lower <- exp(pred_link$fit - 1.96*pred_link$se.fit - ref_logit)
new_dat$or_upper <- exp(pred_link$fit + 1.96*pred_link$se.fit - ref_logit)

步骤3：绘制X1对应的优势比曲线 #

用ggplot2绘图的示例代码：

ggplot(new_dat, aes(x = x1, y = or)) +
  geom_line(linewidth = 1, color = "steelblue") +
  geom_ribbon(aes(ymin = or_lower, ymax = or_upper), alpha = 0.2, fill = "steelblue") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red") + # OR=1的参考线，代表无效应
  labs(x = "X1取值", y = "优势比(OR, 参照组: X1=0)", title = "X1对应优势比变化曲线") +
  theme_bw()

注意：如果x1的实际取值范围不包含0，建议选择x1的中位数或者观测范围内的典型值作为参照组，避免外推带来的误差。

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