如何在R中绘制鸟类性别判别分析的1维概率轮廓线？

你要的不是2D等高线，而是逻辑回归模型下的概率等值直线——Green & Theobald 1989里的分割线本质是线性模型的概率阈值线，直接用模型系数推导方程即可，不用geom_contour。

核心步骤 #

1. 拟合二元逻辑回归模型（性别为因变量，形态特征为自变量）
2. 推导对应概率的直线方程
3. 在ggplot散点图上添加这些直线

具体代码示例 #

假设你的Dunlin数据框名为dunlin，包含：

• sex：性别因子（"Male"/"Female"）
• bill_length：喙长（连续变量）
• wing_length：翅长（连续变量）

1. 拟合模型

# 将性别转为二元变量（1=雄性，0=雌性）
dunlin$male <- as.numeric(dunlin$sex == "Male")
# 拟合逻辑回归
glm_fit <- glm(male ~ bill_length + wing_length, data = dunlin, family = binomial(link = "logit"))

2. 计算概率阈值对应的直线

逻辑回归的概率公式为：
$$P(\text{Male}=1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 \times \text{bill_length} + \beta_2 \times \text{wing_length})}}$$

整理成直线方程（以bill_length为x轴，wing_length为y轴）：
$$\text{wing_length} = \frac{\text{logit}(p) - \beta_0 - \beta_1 \times \text{bill_length}}{\beta_2}$$

其中logit(p) = log(p/(1-p))，对应p=0.01、0.5、0.99的logit值分别为≈-4.595、0、4.595。

# 提取模型系数
coefs <- coef(glm_fit)
# 生成x轴的序列（覆盖数据中喙长的范围）
bill_seq <- seq(min(dunlin$bill_length), max(dunlin$bill_length), length.out = 100)

# 计算各概率对应的翅长
# P(Male)=0.01
logit_01 <- log(0.01 / 0.99)
wing_01 <- (logit_01 - coefs[1] - coefs[2] * bill_seq) / coefs[3]
# P(Male)=0.5（中性分割线）
logit_05 <- log(0.5 / 0.5) # 等于0
wing_05 <- (logit_05 - coefs[1] - coefs[2] * bill_seq) / coefs[3]
# P(Male)=0.99
logit_99 <- log(0.99 / 0.01)
wing_99 <- (logit_99 - coefs[1] - coefs[2] * bill_seq) / coefs[3]

3. 绘图（散点+分割线）

library(ggplot2)
ggplot(dunlin, aes(x = bill_length, y = wing_length, color = sex)) +
  geom_point(alpha = 0.7) +
  # 添加三条分割线
  geom_line(aes(x = bill_seq, y = wing_01, color = NULL), linetype = "dashed", color = "red") +
  geom_line(aes(x = bill_seq, y = wing_05, color = NULL), linetype = "solid", color = "black", size = 1) +
  geom_line(aes(x = bill_seq, y = wing_99, color = NULL), linetype = "dashed", color = "blue") +
  # 添加标注
  annotate("text", x = max(bill_seq), y = wing_01[length(wing_01)],
           label = "P(Male) = 0.01", hjust = 1, vjust = 1, color = "red") +
  annotate("text", x = max(bill_seq), y = wing_05[length(wing_05)],
           label = "P(Male) = 0.5", hjust = 1, vjust = 1, color = "black") +
  annotate("text", x = max(bill_seq), y = wing_99[length(wing_99)],
           label = "P(Male) = 0.99", hjust = 1, vjust = 1, color = "blue") +
  labs(x = "喙长(mm)", y = "翅长(mm)", color = "性别") +
  theme_bw()

关键说明 #

• 为什么不用geom_contour？因为geom_contour是针对连续的2D预测曲面（比如用predict生成所有(x,y)组合的概率后画等高线），但线性逻辑回归的概率等值线本身就是直线，直接推导方程更简洁，也完全匹配Green & Theobald的图。
• 如果你的模型是线性判别分析（LDA），推导逻辑类似，只是系数来自LDA模型而非逻辑回归。
• 如果需要交换x/y轴，只需调整直线方程的整理方式即可。

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