如何为含多平滑项的二项式广义可加模型（GAM）自定义绘制逻辑概率图

嘿，我完全懂你的需求啦——你想用ggplot手动绘制多平滑项二项式GAM的逻辑概率图，要和base R里plot(fit, trans = plogis, shift = coef(fit)[1])的效果一模一样，对吧？别担心，咱们一步步来实现：

第一步：先确认模型拟合 #

首先咱们先把你给出的多平滑项GAM模型跑起来：

# 加载所需库
library(mgcv)
library(gamair)
library(tidyverse)

# 加载数据集
data("wesdr")

# 拟合含多个平滑项的二项式GAM
fit <- gam(
  ret ~ s(dur) + s(bmi),
  method = "REML",
  family = binomial,
  data = wesdr
)

第二步：处理平滑项估计，转换到概率尺度 #

这里的核心是：我们需要把平滑项的部分效应加上模型截距，再用plogis()函数转换为概率尺度——这和base R里plot.gam的shift+trans参数效果完全一致。因为mgcv的平滑项是中心化的，当其他协变量取均值时，它们的平滑项贡献为0，所以总线性预测就是截距加上目标平滑项的估计，转成概率就是我们要的结果。

# 获取平滑项估计并添加置信区间
sm <- smooth_estimates(fit) %>%
  add_confint()

# 提取模型的截距项
intercept <- coef(fit)[1]

# 针对目标平滑项（比如s(dur)）处理：转换到概率尺度
sm_dur <- sm %>%
  filter(smooth == "s(dur)") %>%
  mutate(
    # 计算线性预测尺度的总效应：截距 + 平滑项估计
    lin_pred_est = intercept + est,
    lin_pred_lower = intercept + lower_ci,
    lin_pred_upper = intercept + upper_ci,
    # 用plogis转换为概率
    prob_est = plogis(lin_pred_est),
    prob_lower = plogis(lin_pred_lower),
    prob_upper = plogis(lin_pred_upper)
  )

# 如果要处理另一个平滑项s(bmi)，只需要把filter里的条件改成smooth == "s(bmi)"即可

如果你想把部分残差也转换到概率尺度（可选），可以这么做：

# 给原数据添加部分残差
wesdr <- wesdr %>%
  add_partial_residuals(fit)

# 将s(dur)的部分残差转换为概率尺度
wesdr <- wesdr %>%
  mutate(pr_prob_dur = plogis(pr_s(dur) + intercept))

第三步：绘制自定义概率图 #

现在用ggplot把结果画出来，完全还原base R的效果：

p <- ggplot(sm_dur, aes(x = dur)) +
  # 添加x轴的rug图，展示原数据的dur分布
  geom_rug(data = wesdr, sides = "b", length = grid::unit(0.02, "npc")) +
  # 绘制概率的置信区间
  geom_ribbon(aes(ymin = prob_lower, ymax = prob_upper), alpha = 0.2) +
  # 绘制概率的拟合线
  geom_line(aes(y = prob_est), lwd = 1.2) +
  # 可选：添加转换后的部分残差点（如果需要的话）
  # geom_point(data = wesdr, aes(y = pr_prob_dur), alpha = 0.3, size = 0.8) +
  labs(
    y = "预测概率（ret=1）",
    title = "平滑项s(dur)对应的ret概率效应（其他协变量取均值）",
    x = "糖尿病病程（dur）"
  ) +
  theme_minimal()

# 展示绘图结果
p

额外说明 #

• 这个方法适用于任意数量的平滑项：只要你想绘制某个平滑项的概率效应，只需要过滤出对应的平滑项，重复上面的转换步骤即可。
• 这里的“其他协变量取均值”是mgcv默认的行为，和plot.gam的逻辑一致——因为平滑项是中心化的，当其他协变量处于均值水平时，它们的平滑项贡献为0，所以我们只需要加截距就能得到正确的概率估计。

备注：内容来源于stack exchange，提问作者Shawn Hemelstrand