tidymodels中logistic_reg()系数估计与实际趋势不符的疑问

问题描述 #

我正在使用tidymodels基于Palmer企鹅数据集开发逻辑回归模型，希望将female作为结果变量的参考类，模型的阳性结果设为female。但模型系数估计显示bill_depth_mm和bill_length_mm的增加会提升企鹅为female的概率，而可视化结果却呈现相反趋势。不过计算特异性（specificity）和敏感性（sensitivity）时，模型将female视为真阳性结果，ROC曲线也要求参考类对应female的预测概率。请问：

1. 为何模型估计结果看似反转？
2. tidymodels中是否有选项可指定或确认系数对应的预测类别？
3. 还是我的模型理解有误，实际上喙长和喙深的增加确实会提升企鹅为雌性的概率？

附代码实现：

penguins <- read_csv("Data/penguins.csv")

penguins <- penguins %>% 
  mutate(sex = factor(sex))

set.seed(57475)
penguins_splits <- initial_split(penguins,
                                 strata = sex)

penguins_training <- training(penguins_splits) 
penguins_testing <- testing(penguins_splits)

penguins_recipe <- recipe(sex ~ ., data = penguins_training) %>% 
  step_naomit(all_predictors()) %>%
  step_rm(rowid, skip = TRUE) %>%
  step_string2factor(all_nominal_predictors()) %>%
  step_mutate(year = factor(year)) %>%
  step_dummy(all_nominal_predictors())


logistic_model <- logistic_reg() %>% 
  set_engine("glm") %>%
  set_mode("classification")

logistic_wf <- workflow() %>% 
  add_model(logistic_model) %>% 
  add_recipe(penguins_recipe)

logistic_fit <- logistic_wf %>%
  fit(data = penguins_training)

logistic_fit %>%
  pull_workflow_fit %>%
  tidy()

ggplot(penguins_training, aes(x = bill_length_mm, colour = sex)) +
  geom_boxplot()

模型系数输出：

> logistic_fit %>%
+   pull_workflow_fit %>%
+   tidy()
# A tibble: 11 x 5
   term               estimate std.error statistic      p.value
   <chr>                 <dbl>     <dbl>     <dbl>        <dbl>
 1 (Intercept)       -76.7      13.7        -5.58  0.0000000237
 2 bill_length_mm      0.831     0.176       4.72  0.00000238  
 3 bill_depth_mm       1.53      0.391       3.91  0.0000913   
 4 flipper_length_mm  -0.0290    0.0611     -0.474 0.635       
 5 body_mass_g         0.00593   0.00132     4.49  0.00000724  
 6 species_Chinstrap  -9.17      2.06       -4.45  0.00000843  
 7 species_Gentoo     -8.82      3.22       -2.74  0.00608     
 8 island_Dream        0.874     0.968       0.903 0.367       
 9 island_Torgersen   -0.237     0.974      -0.243 0.808       
10 year_X2008         -0.239     0.701      -0.341 0.733       
11 year_X2009         -0.762     0.736      -1.04  0.300  

解答 #

1. 模型系数与可视化反转的核心原因 #

你误解了glm逻辑回归默认的预测类别对应关系：

• R中默认将因子按字母顺序排序，sex因子的水平顺序是["female", "male"]，即female为第一水平（参考类），male为第二水平。
• glm的逻辑回归模型默认建模的是第二水平的发生概率，对应的logit方程为：

  log(P(male) / P(female)) = 截距 + bill_length_mm*0.831 + bill_depth_mm*1.53 + ...
  

  系数为正意味着变量增加时，male的概率上升——这和你的箱线图趋势完全一致（雄性企鹅的喙更长、更深），并非模型估计错误。你之前误以为系数对应female的概率，所以产生了“反转”的错觉。

2. tidymodels中指定/确认系数对应的预测类别 #

有两种常用方法来控制模型对应的预测类别：

• 调整因子水平顺序：在数据预处理阶段，将你想要作为模型预测目标（阳性结果）的类别设为因子的第二水平。比如要让模型预测female的概率，可修改数据处理步骤：

  penguins <- penguins %>% 
    mutate(sex = factor(sex, levels = c("male", "female")))
  

  或者在recipe中加入step_relevel：

  penguins_recipe <- recipe(sex ~ ., data = penguins_training) %>% 
    # ...其他步骤...
    step_relevel(sex, levels = c("male", "female"))
  

  此时模型的logit方程会变为log(P(female)/P(male))，系数符号会和你预期的一致。
• 验证当前因子水平：直接查看因子水平确认模型对应的参考类和预测类：

  levels(penguins_training$sex)
  

另外，关于评估指标（敏感性/特异性、ROC曲线）的阳性类别指定，可在yardstick包的函数中通过event_level参数设置，比如：

sensitivity(results, truth = sex, estimate = .pred_class, event_level = "second")

3. 关于变量与性别的关系 #

你的可视化结论是正确的：雄性企鹅的喙长、喙深整体大于雌性。因此这些变量的增加，实际上是提升企鹅为雄性的概率，而非雌性。之前的矛盾完全是因为对模型预测类别的误解。

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