第 58 章 logistic回归模型

本章讲广义线性模型中的logistic回归模型。

58.1 问题

假定这里有一组数据,包含学生GRE成绩和被录取的状态(admit = 1,就是被录取;admit = 0,没有被录取)

library(tidyverse)
gredata <- read_csv("demo_data/gredata.csv")
gredata
## # A tibble: 400 × 4
##   admit   gre   gpa  rank
##   <dbl> <dbl> <dbl> <dbl>
## 1     0   380  3.61     3
## 2     1   660  3.67     3
## 3     1   800  4        1
## 4     1   640  3.19     4
## 5     0   520  2.93     4
## 6     1   760  3        2
## # … with 394 more rows

我们能用学生GRE的成绩预测录取状态吗?回答这个问题需要用到logistic回归模型,

$$ \[\begin{align*} \text{admit}_{i} &\sim \mathrm{Binomial}(1, p_{i}) \\ \text{logit}(p_{i}) &= \log\Big(\frac{p_{i}}{1 - p_{i}}\Big) = \alpha + \beta \cdot \text{gre}_{i} \\ \text{equivalent to,} \quad p_{i} &= \frac{1}{1 + \exp[- (\alpha + \beta \cdot \text{gre}_{i})]} \\ & = \frac{\exp(\alpha + \beta \cdot \text{gre}_{i})}{1 + \exp (\alpha + \beta \cdot\text{gre}_{i})} \\ \end{align*}\] $$

这里 \(p_i\) 就是被录取的概率。预测因子 \(gre\) 的线性组合模拟的是 \(log\Big(\frac{p_{i}}{1 - p_{i}}\Big)\) ,即对数比率(log-odds).

按照上面表达式,用glm函数写代码,

model_logit <- glm(admit ~ gre,
  data = gredata,
  family = binomial(link = "logit")
)

summary(model_logit)
## 
## Call:
## glm(formula = admit ~ gre, family = binomial(link = "logit"), 
##     data = gredata)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.162  -0.905  -0.755   1.349   1.988  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.901344   0.606038   -4.79  1.7e-06 ***
## gre          0.003582   0.000986    3.63  0.00028 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 486.06  on 398  degrees of freedom
## AIC: 490.1
## 
## Number of Fisher Scoring iterations: 4

得到gre的系数是

coef(model_logit)[2]
##      gre 
## 0.003582

怎么理解这个0.003582呢?

58.2 模型的输出

为了更好地理解模型的输出,这里用三种不同的度量方式(scales)来计算系数。

假定 \(p\) 为录取的概率

num scale formula
1 The log-odds scale \(log \frac{p}{1 - p}\)
2 The odds scale \(\frac{p}{1 -p}\)
3 The probability scale \(p\)

58.2.1 The log-odds scale

模型给出系数(0.003582)事实上就是log-odds度量方式的结果,具体来说, 这里系数(0.003582)代表着: GRE考试成绩每增加1个单位,那么log-odds(录取概率)就会增加0.003582. (注意,不是录取概率增加0.003582,而是log-odds(录取概率)增加0.003582)

为了更清楚的理解,我们把log-odds的结果与gre成绩画出来看看

logit_log_odds <- broom::augment_columns(
  model_logit,
  data = gredata,
  type.predict = c("link")
) %>%
  rename(log_odds = .fitted) 
library(latex2exp)
logit_log_odds %>% 
    ggplot(aes(x = gre, y = log_odds)) +
    geom_path(color = "#771C6D", size = 2) +
    labs(title = "Log odds", 
        subtitle = "This is linear!",
        x = NULL,
        y = TeX("$log \\frac{p}{1 - p}$")) +
    theme_minimal() +
    theme(
      plot.title = element_text(face = "bold"),
      axis.title.y = element_text(angle = 90)
          )

由图看到,GRE成绩 与log-odds(录取概率)这个值的关系是线性的,斜率就是模型给出的系数0.003582

58.2.2 The odds scale

第二种odds scale度量方式可能要比第一种要好理解点,我们先求系数的指数,exp(0.003582) = 1.003588

exp(0.003582)
## [1] 1.004

1.003588的含义是: GRE考试成绩每增加1个单位,那么odds(录取概率)就会增大1.003588倍;若增加2个单位,那么odds(录取概率)就会增大(1.003588 * 1.003588)倍,也就说是个乘法关系。

有时候,大家喜欢用增长百分比表述,那么就是

(exp(0.003582) - 1) x 100% = (1.003588 - 1) x 100% = 0.36%

即,GRE考试成绩每增加1个点,那么odds(录取概率)就会增长百分之0.36.

同样,我们把odds(录取概率)的结果与GRE成绩画出来看看

logit_odds <- broom::augment_columns(
  model_logit,
  data = gredata,
  type.predict = c("link")
) %>%
  rename(log_odds = .fitted) %>%
  mutate(odds_ratio = exp(log_odds))
logit_odds %>% 
    ggplot(aes(x = gre, y = odds_ratio)) +
    geom_line(color = "#FB9E07", size = 2) +
    labs(title = "Odds", 
        subtitle = "This is curvy, but it's a mathy transformation of a linear value",
        x = NULL,
        y = TeX("$\\frac{p}{1 - p}$")) +
    theme_minimal() +
    theme(
      plot.title = element_text(face = "bold"),
      axis.title.y = element_text(angle = 90)
    )

58.2.3 The probability scale.

第三种度量方式是概率度量(probability scale),因为模型假定的是,GRE的分数与log-odds(录取概率)呈线性关系,那么很显然GRE的分数与录取概率就不可能是线性关系了,而是呈非线性关系。我们先看下非线性关系长什么样。

logit_probs <- broom::augment_columns(
  model_logit,
  data = gredata,
  type.predict = c("response")
) %>% 
  rename(pred_prob = .fitted)
logit_probs %>% 
    ggplot(aes(x = gre, y = pred_prob)) +
    #geom_point(aes(x = gre, y = admit)) +
    geom_line(color = "#CF4446", size = 2) +
    labs(title = "Predicted probabilities", 
        sutitle = "Plug values of X into ",
        x = "X (value of explanatory variable)",
        y = TeX("\\hat{P(Y)} ")) +
    theme_minimal() +
    theme(plot.title = element_text(face = "bold"))

可以看到,GRE分数对录取的概率的影响是正的且非线性的,具体来说,

  • GRE分数200分左右,录取概率约0.1;
  • GRE分数500分左右,录取概率约0.25;
  • GRE分数800分左右,录取概率接近0.5;

提请注意的是,以上三种度量的方式中:

  • log_odds scale,预测因子与log_odds()的关系是一个固定值,具有可加性
  • odds scale, 预测因子与odds()的关系是一个固定值,具有乘法性
  • probability,预测因子与probability的关系不再是一个固定值了

用哪种度量方式来理解模型的输出,取决不同的场景。第一种方式容易计算但理解较为困难,第三种方式最容易理解,但不再是线性关系了。

58.3 预测

58.3.1 预测与拟合

先认识下两个常用的函数predict()fitted().

gredata %>%
  mutate(
    pred = predict(model_logit),
    fitted = fitted(model_logit)
    )
## # A tibble: 400 × 6
##   admit   gre   gpa  rank    pred fitted
##   <dbl> <dbl> <dbl> <dbl>   <dbl>  <dbl>
## 1     0   380  3.61     3 -1.54    0.177
## 2     1   660  3.67     3 -0.537   0.369
## 3     1   800  4        1 -0.0356  0.491
## 4     1   640  3.19     4 -0.609   0.352
## 5     0   520  2.93     4 -1.04    0.261
## 6     1   760  3        2 -0.179   0.455
## # … with 394 more rows

线性模型中,predict()fitted()这种写法的返回结果是一样的。但在glm模型中,两者的结果是不同的。predict()返回的是log_odds(录取概率)度量的结果;而fitted()返回的是录取概率。如果想保持一致,需要对predict()返回结果做反向的log_odds计算

\[p = \exp(\alpha) / (1 + \exp(\alpha) )\]

具体如下

gredata %>%
  mutate(
    pred = predict(model_logit),
    fitted = fitted(model_logit),
    pred2 = exp(pred) / (1 + exp(pred) )
    )
## # A tibble: 400 × 7
##   admit   gre   gpa  rank    pred fitted pred2
##   <dbl> <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>
## 1     0   380  3.61     3 -1.54    0.177 0.177
## 2     1   660  3.67     3 -0.537   0.369 0.369
## 3     1   800  4        1 -0.0356  0.491 0.491
## 4     1   640  3.19     4 -0.609   0.352 0.352
## 5     0   520  2.93     4 -1.04    0.261 0.261
## 6     1   760  3        2 -0.179   0.455 0.455
## # … with 394 more rows

我这样折腾无非是想让大家知道,在glm中predict()fit()是不同的。

如果想让predict()也返回录取概率,也可以不用那么麻烦,事实上predicttype = "response") 选项已经为我们准备好了。

gredata %>%
  mutate(
    pred = predict(model_logit, type = "response"),
    fitted = fitted(model_logit)
    )
## # A tibble: 400 × 6
##   admit   gre   gpa  rank  pred fitted
##   <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     0   380  3.61     3 0.177  0.177
## 2     1   660  3.67     3 0.369  0.369
## 3     1   800  4        1 0.491  0.491
## 4     1   640  3.19     4 0.352  0.352
## 5     0   520  2.93     4 0.261  0.261
## 6     1   760  3        2 0.455  0.455
## # … with 394 more rows

58.3.2 预测

有时候,我们需要对给定的GRE分数,用建立的模型预测被录取的概率

newdata <- tibble(
  gre = c(550, 660, 700, 780)
)
newdata
## # A tibble: 4 × 1
##     gre
##   <dbl>
## 1   550
## 2   660
## 3   700
## 4   780

前面讲到predict()type参数有若干选项type = c("link", "response", "terms"), 可参考第 60

  • type = "link",预测的是log_odds,实际上就是coef(model_logit)[1] + gre * coef(model_logit)[2]
newdata %>% 
  mutate(
    pred_log_odds = predict(model_logit, newdata = newdata, type = "link"),
    #
    pred_log_odds2 = coef(model_logit)[1] + gre * coef(model_logit)[2]
    
  )
## # A tibble: 4 × 3
##     gre pred_log_odds pred_log_odds2
##   <dbl>         <dbl>          <dbl>
## 1   550        -0.931         -0.931
## 2   660        -0.537         -0.537
## 3   700        -0.394         -0.394
## 4   780        -0.107         -0.107
  • type = "response"",预测的是probabilities, 实际上就是exp(pred_log_odds) / (1 + exp(pred_log_odds) )
newdata %>% 
  mutate(
    pred_log_odds = predict(model_logit, newdata = newdata, type = "link")
    ) %>% 
  mutate(
    pred_prob = predict(model_logit, newdata = newdata, type = "response"),
    #
    pred_prob2 = exp(pred_log_odds) / (1 + exp(pred_log_odds) )
  )
## # A tibble: 4 × 4
##     gre pred_log_odds pred_prob pred_prob2
##   <dbl>         <dbl>     <dbl>      <dbl>
## 1   550        -0.931     0.283      0.283
## 2   660        -0.537     0.369      0.369
## 3   700        -0.394     0.403      0.403
## 4   780        -0.107     0.473      0.473
  • type = "terms",返回一个矩阵,具体计算为:模型的设计矩阵中心化以后,与模型返回的系数相乘而得到的新矩阵10
predict(model_logit, gredata, type = 'terms') %>% 
  as_tibble() %>% 
  head()
## # A tibble: 6 × 1
##      gre
##    <dbl>
## 1 -0.744
## 2  0.259
## 3  0.761
## 4  0.187
## 5 -0.243
## 6  0.617

我们复盘一次

X <- model.matrix(admit ~ gre, data = gredata)

X %>% 
  as.data.frame() %>% 
  mutate(
    across(everything(), ~ .x - mean(.x))
  ) %>% 
  transmute(
    term = coef(model_logit)[2] * gre
    ) %>% 
  head()
##      term
## 1 -0.7440
## 2  0.2590
## 3  0.7605
## 4  0.1873
## 5 -0.2425
## 6  0.6172