如何在gratia包compare_smooths绘图中添加原始数据并转换y轴为概率

1. 将y轴转换为logistic概率（0-1区间） #

compare_smooths()返回的平滑项估计处于线性预测器尺度（即logit尺度），要转换为直观的0-1概率尺度，需要结合模型截距与logistic逆链接函数处理：

简洁实现方式 #

直接在draw()中传递fun和constant参数：fun指定逆链接函数plogis()，constant传入每个模型的截距值，gratia会自动为对应模型的平滑项加上截距后完成概率转换：

draw(comp, 
     fun = plogis, 
     constant = list(m1 = coef(m1)[1], m2 = coef(m2)[1])) +
  labs(y = "Response Probability")

也可以用inv_link(m1)替代plogis()，效果完全一致，更贴合模型的链接函数定义。

2. 添加原始数据到对比图中 #

先整理两个模型的原始数据集，添加模型标识后，用ggplot的几何图层叠加到平滑曲线之上：

# 整理原始数据，添加模型区分列
raw_data <- dplyr::bind_rows(
  df1 %>% dplyr::mutate(model = "m1"),
  df2 %>% dplyr::mutate(model = "m2")
)

然后在绘图时叠加抖动散点（避免0/1数据点完全重叠）：

draw(comp, 
     fun = plogis, 
     constant = list(m1 = coef(m1)[1], m2 = coef(m2)[1])) +
  ggplot2::geom_jitter(data = raw_data, 
              aes(x = .data[[unique(comp$variable)]], y = response, color = model),
              alpha = 0.3, width = 0.2, height = 0.02) +
  labs(y = "Response Probability")

• .data[[unique(comp$variable)]]：动态匹配当前子图的x变量（如p1或p2），无需针对不同平滑项单独修改代码
• geom_jitter：通过轻微抖动避免0/1数据点重叠，alpha设置透明度防止遮挡平滑曲线

完整可运行代码 #

library(gratia)
library(dplyr)
library(ggplot2)

# simulate data
set.seed(9)
df1 <- data.frame(response=c(rep(0,500), rep(1,500)),
                  p1 = c(rnorm(500,100,5),rnorm(500,110,3)),
                  p2 = rnorm(1000,10,2))


df2 <- data.frame(response=sample(c(0,1),1000, replace = T),
                  p1 = rnorm(1000,80,10),
                  p2 = rnorm(1000,15,1))

# create models
m1 <- gam(response ~ s(p1,bs="cr") + s(p2, bs="cr"),
          data = df1, family=binomial, method = "REML", select=T)
m2 <- gam(response ~ s(p1,bs="cr") + s(p2, bs="cr"),
          data = df2, family=binomial, method = "REML", select=T)

# compare smooths
comp <- compare_smooths(m1,m2)

# 整理原始数据
raw_data <- bind_rows(
  df1 %>% mutate(model = "m1"),
  df2 %>% mutate(model = "m2")
)

# 绘制带概率y轴和原始数据的对比图
draw(comp, 
     fun = inv_link(m1), 
     constant = list(m1 = coef(m1)[1], m2 = coef(m2)[1])) +
  geom_jitter(data = raw_data, 
              aes(x = .data[[unique(comp$variable)]], y = response, color = model),
              alpha = 0.3, width = 0.2, height = 0.02) +
  labs(y = "Response Probability") +
  theme_bw()

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