逻辑回归

分类问题建模

统计建模篇

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

问题

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

library(tidyverse)
gredata <- read_csv("demo_data/gredata.csv")
gredata

我们能用学生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)

得到gre的系数是

coef(model_logit)[2]

怎么理解这个0.003582呢？

模型的输出

为了更好地理解模型的输出，这里用三种不同的度量方式(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$

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

The odds scale

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

exp(0.003582)

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)
    )

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；

::: {.rmdnote} 提请注意的是，以上三种度量的方式中：

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

:::

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

预测

预测与拟合

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

gredata %>%
  mutate(
    pred = predict(model_logit),
    fitted = fitted(model_logit)
    )

线性模型中，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) )
    )

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

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

gredata %>%
  mutate(
    pred = predict(model_logit, type = "response"),
    fitted = fitted(model_logit)
    )

预测

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

newdata <- tibble(
  gre = c(550, 660, 700, 780)
)
newdata

前面讲到predict()中type参数有若干选项type = c("link", "response", "terms"), 可参考相关章节

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]
    
  )

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) )
  )

type = "terms"，返回一个矩阵，具体计算为：模型的设计矩阵中心化以后，与模型返回的系数相乘而得到的新矩阵^[https://stackoverflow.com/questions/37963904/what-does-predict-glm-type-terms-actually-do]。

predict(model_logit, gredata, type = 'terms') %>% 
  as_tibble() %>% 
  head()

我们复盘一次

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()

# remove the objects
# rm(list=ls())
rm(gredata, logit_log_odds, logit_odds, logit_probs, model_logit, newdata, X)

pacman::p_unload(pacman::p_loaded(), character.only = TRUE)