如何用logistic regression绘制概率图？复刻目标图及计算育幼概率困惑

一、需求与数据集 #

需要复刻的目标图展示了不同P_Treatment组下，胚胎存在与否（二分类变量Embryo_Presence_vs._Absence）随连续变量Final_Size变化的logistic回归关系。现有数据集如下：

print.data.frame(x[c("P_Treatment", "Embryo_Presence_vs._Absence", "Final_Size")])

P_Treatment Embryo_Presence_vs._Absence Final_Size
1   Intermediate                           0       3.80
2   Intermediate                           0       3.76
3   Intermediate                           0       3.70
4   Intermediate                           0       3.73
5   Intermediate                           0       3.10
6   Intermediate                           0       3.77
7   Intermediate                           0       3.59
8   Intermediate                           0       3.73
9   Intermediate                           0       3.36
10  Intermediate                           0       3.81
11  Intermediate                           0       3.72
12  Intermediate                           0       3.92
13  Intermediate                           0       3.56
14  Intermediate                           0       3.78
15  Intermediate                           0       3.63
16  Intermediate                           0       3.27
17  Intermediate                           0       3.60
18  Intermediate                           0       3.74
19  Intermediate                           0       3.60
20  Intermediate                           0       3.65
21  Intermediate                           0       3.59
22  Intermediate                           0       3.75
23  Intermediate                           0       3.78
24  Intermediate                           0       3.55
25  Intermediate                           0       3.65
26  Intermediate                           0       3.65
27           Low                           0       3.91
28           Low                           0       3.72
29           Low                           0       3.77
30           Low                           0       3.73
31           Low                           0       3.65
32           Low                           0       3.57
33           Low                           0       3.82
34           Low                           0       3.76
35           Low                           0       3.88
36           Low                           0       4.10
37           Low                           0       3.71
38           Low                           0       3.57
39           Low                           0       3.77
40           Low                           0       3.60
41           Low                           0       3.70
42           Low                           0       3.55
43           Low                           0       4.00
44           Low                           0       3.56
45           Low                           0       3.71
46           Low                           0       3.61
47           Low                           0       3.63
48           Low                           0       3.72
49           Low                           0       3.80
50           Low                           0       3.86
51           Low                           0       3.08
52           Low                           0       3.81
53           Low                           0       3.73
54           Low                           0       3.84
55           Low                           0       3.76
56           Low                           0       3.66
57           Low                           0       3.70
58           Low                           0       3.71
59           Low                           0       3.60
60           Low                           0       3.75
61           Low                           0       3.74
62  Intermediate                           1       3.89
63           Low                           1       3.87
64  Intermediate                           1       3.99
65  Intermediate                           1       3.93
66  Intermediate                           1       3.65
67  Intermediate                           1       3.64
68  Intermediate                           1       3.81
69           Low                           1       3.67
70           Low                           1       3.52
71  Intermediate                           1       3.70
72           Low                           1       3.83
73           Low                           1       3.65
74           Low                           1       3.94
75  Intermediate                           1       3.75
76           Low                           1       3.71
77  Intermediate                           1       3.56
78  Intermediate                           1       3.89
79           Low                           1       3.72
80           Low                           1       3.25
81  Intermediate                           1       3.79
82  Intermediate                           1       3.60
83  Intermediate                           1       3.88
84  Intermediate                           1       3.75
85  Intermediate                           1       3.75
86           Low                           1       3.58
87  Intermediate                           1       3.75
88  Intermediate                           1       3.65
89           Low                           1       3.60
90  Intermediate                           1       3.68
91           Low                           1       3.65
92  Intermediate                           1       3.88
93  Intermediate                           1       3.77
94  Intermediate                           1       3.63
95           Low                           1       3.68
96  Intermediate                           1       3.83
97  Intermediate                           1       3.85
98           Low                           1       3.81
99  Intermediate                           1       3.56
100 Intermediate                           1       3.70
101          Low                           1       3.74
102          Low                           1       3.65
103 Intermediate                           1       3.69
104 Intermediate                           1       3.89
105 Intermediate                           1       3.70
106 Intermediate                           1       3.70
107 Intermediate                           1       3.76
108 Intermediate                           1       3.66
109 Intermediate                           1       3.78
110 Intermediate                           1       4.22
111 Intermediate                           1       3.68
112          Low                           1       3.82
113          Low                           1       3.82
114          Low                           1       3.76
115          Low                           1       3.91
116 Intermediate                           1       3.56
117          Low                           1       3.66
118 Intermediate                           1       3.65
119          Low                           1       3.61
120 Intermediate                           1       3.65
121          Low                           1       4.16
122 Intermediate                           1       3.74
123 Intermediate                           1       3.60
124          Low                           1       3.50
125          Low                           1       3.76
126          Low                           1       3.85
127          Low                           1       3.83
128          Low                           1       3.60

二、当前绘图问题 #

当前使用的ggplot代码未按P_Treatment分组，导致生成的logistic曲线是整体数据的拟合，无法对应目标图的分组展示，代码如下：

ggplot(x, aes(Final_Size, Embryo_Presence_vs._Absence)) +
  geom_jitter(height = 0.05) +
  stat_smooth(method="glm", se= FALSE, method.args = list(family=binomial))

三、复刻目标图的修正代码 #

调整代码，按P_Treatment分组绘制logistic曲线，同时保留置信区间（与目标图一致）：

library(ggplot2)

ggplot(x, aes(x = Final_Size, y = Embryo_Presence_vs._Absence, 
              color = P_Treatment, group = P_Treatment)) +
  geom_jitter(height = 0.05, alpha = 0.6) + # 降低透明度避免点重叠
  stat_smooth(method = "glm", se = TRUE,
              method.args = list(family = binomial),
              linewidth = 1) +
  labs(x = "Final Size", y = "Probability of Embryo Presence",
       color = "P Treatment") +
  theme_classic()

四、育幼概率计算方法 #

1. 基础R：predict()函数直接计算 #

先拟合包含处理组和最终大小的logistic模型，再生成预测概率：

# 拟合主效应模型，若需考虑交互项可改为 ~ P_Treatment * Final_Size
glm_model <- glm(Embryo_Presence_vs._Absence ~ P_Treatment + Final_Size,
                 data = x, family = binomial)

# 构造预测数据集：覆盖不同处理组和Final_Size范围
new_data <- expand.grid(
  P_Treatment = unique(x$P_Treatment),
  Final_Size = seq(min(x$Final_Size), max(x$Final_Size), length.out = 100)
)

# 预测概率（*type = "response"* 确保输出概率值，而非logit转换值）
new_data$pred_prob <- predict(glm_model, newdata = new_data, type = "response")

# 查看前6行预测结果
head(new_data)

2. emmeans包：计算特定条件下的边际概率 #

适合计算如“不同处理组在某一Final_Size值下的平均育幼概率”这类需求：

library(emmeans)

# 拟合含交互项的模型
glm_model2 <- glm(Embryo_Presence_vs._Absence ~ P_Treatment * Final_Size,
                  data = x, family = binomial)

# 指定Final_Size值，计算各处理组的边际概率
emm <- emmeans(glm_model2, ~ P_Treatment | Final_Size, 
               at = list(Final_Size = c(3.5, 4.0)),
               type = "response")

# 查看结果
summary(emm)

3. ggeffects包：快速生成预测数据与可视化 #

该包可直接生成适配ggplot的预测结果，同时输出概率数据：

library(ggeffects)

# 生成预测数据：按Final_Size和P_Treatment分组
gg_pred <- ggpredict(glm_model, terms = c("Final_Size", "P_Treatment"))

# 查看预测概率数据
head(gg_pred)

# 直接绘制带原始点的预测图
plot(gg_pred) + 
  geom_jitter(data = x, aes(x = Final_Size, y = Embryo_Presence_vs._Absence, color = P_Treatment), 
              height = 0.05, alpha = 0.6)

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