diff --git a/_bookdown.yml b/_bookdown.yml index 346407e..70b804e 100644 --- a/_bookdown.yml +++ b/_bookdown.yml @@ -50,7 +50,8 @@ rmd_files: [ "tidystats.Rmd", "tests_as_linear.Rmd", "lmm.Rmd", - "glm.Rmd", + "poisson_regression.Rmd", + "logistic_regression.Rmd", "ordinal.Rmd", "tidymodels.Rmd", diff --git a/demo_data/gredata.csv b/demo_data/gredata.csv new file mode 100644 index 0000000..5f2cf4e --- /dev/null +++ b/demo_data/gredata.csv @@ -0,0 +1,401 @@ +admit,gre,gpa,rank +0,380,3.61,3 +1,660,3.67,3 +1,800,4,1 +1,640,3.19,4 +0,520,2.93,4 +1,760,3,2 +1,560,2.98,1 +0,400,3.08,2 +1,540,3.39,3 +0,700,3.92,2 +0,800,4,4 +0,440,3.22,1 +1,760,4,1 +0,700,3.08,2 +1,700,4,1 +0,480,3.44,3 +0,780,3.87,4 +0,360,2.56,3 +0,800,3.75,2 +1,540,3.81,1 +0,500,3.17,3 +1,660,3.63,2 +0,600,2.82,4 +0,680,3.19,4 +1,760,3.35,2 +1,800,3.66,1 +1,620,3.61,1 +1,520,3.74,4 +1,780,3.22,2 +0,520,3.29,1 +0,540,3.78,4 +0,760,3.35,3 +0,600,3.4,3 +1,800,4,3 +0,360,3.14,1 +0,400,3.05,2 +0,580,3.25,1 +0,520,2.9,3 +1,500,3.13,2 +1,520,2.68,3 +0,560,2.42,2 +1,580,3.32,2 +1,600,3.15,2 +0,500,3.31,3 +0,700,2.94,2 +1,460,3.45,3 +1,580,3.46,2 +0,500,2.97,4 +0,440,2.48,4 +0,400,3.35,3 +0,640,3.86,3 +0,440,3.13,4 +0,740,3.37,4 +1,680,3.27,2 +0,660,3.34,3 +1,740,4,3 +0,560,3.19,3 +0,380,2.94,3 +0,400,3.65,2 +0,600,2.82,4 +1,620,3.18,2 +0,560,3.32,4 +0,640,3.67,3 +1,680,3.85,3 +0,580,4,3 +0,600,3.59,2 +0,740,3.62,4 +0,620,3.3,1 +0,580,3.69,1 +0,800,3.73,1 +0,640,4,3 +0,300,2.92,4 +0,480,3.39,4 +0,580,4,2 +0,720,3.45,4 +0,720,4,3 +0,560,3.36,3 +1,800,4,3 +0,540,3.12,1 +1,620,4,1 +0,700,2.9,4 +0,620,3.07,2 +0,500,2.71,2 +0,380,2.91,4 +1,500,3.6,3 +0,520,2.98,2 +0,600,3.32,2 +0,600,3.48,2 +0,700,3.28,1 +1,660,4,2 +0,700,3.83,2 +1,720,3.64,1 +0,800,3.9,2 +0,580,2.93,2 +1,660,3.44,2 +0,660,3.33,2 +0,640,3.52,4 +0,480,3.57,2 +0,700,2.88,2 +0,400,3.31,3 +0,340,3.15,3 +0,580,3.57,3 +0,380,3.33,4 +0,540,3.94,3 +1,660,3.95,2 +1,740,2.97,2 +1,700,3.56,1 +0,480,3.13,2 +0,400,2.93,3 +0,480,3.45,2 +0,680,3.08,4 +0,420,3.41,4 +0,360,3,3 +0,600,3.22,1 +0,720,3.84,3 +0,620,3.99,3 +1,440,3.45,2 +0,700,3.72,2 +1,800,3.7,1 +0,340,2.92,3 +1,520,3.74,2 +1,480,2.67,2 +0,520,2.85,3 +0,500,2.98,3 +0,720,3.88,3 +0,540,3.38,4 +1,600,3.54,1 +0,740,3.74,4 +0,540,3.19,2 +0,460,3.15,4 +1,620,3.17,2 +0,640,2.79,2 +0,580,3.4,2 +0,500,3.08,3 +0,560,2.95,2 +0,500,3.57,3 +0,560,3.33,4 +0,700,4,3 +0,620,3.4,2 +1,600,3.58,1 +0,640,3.93,2 +1,700,3.52,4 +0,620,3.94,4 +0,580,3.4,3 +0,580,3.4,4 +0,380,3.43,3 +0,480,3.4,2 +0,560,2.71,3 +1,480,2.91,1 +0,740,3.31,1 +1,800,3.74,1 +0,400,3.38,2 +1,640,3.94,2 +0,580,3.46,3 +0,620,3.69,3 +1,580,2.86,4 +0,560,2.52,2 +1,480,3.58,1 +0,660,3.49,2 +0,700,3.82,3 +0,600,3.13,2 +0,640,3.5,2 +1,700,3.56,2 +0,520,2.73,2 +0,580,3.3,2 +0,700,4,1 +0,440,3.24,4 +0,720,3.77,3 +0,500,4,3 +0,600,3.62,3 +0,400,3.51,3 +0,540,2.81,3 +0,680,3.48,3 +1,800,3.43,2 +0,500,3.53,4 +1,620,3.37,2 +0,520,2.62,2 +1,620,3.23,3 +0,620,3.33,3 +0,300,3.01,3 +0,620,3.78,3 +0,500,3.88,4 +0,700,4,2 +1,540,3.84,2 +0,500,2.79,4 +0,800,3.6,2 +0,560,3.61,3 +0,580,2.88,2 +0,560,3.07,2 +0,500,3.35,2 +1,640,2.94,2 +0,800,3.54,3 +0,640,3.76,3 +0,380,3.59,4 +1,600,3.47,2 +0,560,3.59,2 +0,660,3.07,3 +1,400,3.23,4 +0,600,3.63,3 +0,580,3.77,4 +0,800,3.31,3 +1,580,3.2,2 +1,700,4,1 +0,420,3.92,4 +1,600,3.89,1 +1,780,3.8,3 +0,740,3.54,1 +1,640,3.63,1 +0,540,3.16,3 +0,580,3.5,2 +0,740,3.34,4 +0,580,3.02,2 +0,460,2.87,2 +0,640,3.38,3 +1,600,3.56,2 +1,660,2.91,3 +0,340,2.9,1 +1,460,3.64,1 +0,460,2.98,1 +1,560,3.59,2 +0,540,3.28,3 +0,680,3.99,3 +1,480,3.02,1 +0,800,3.47,3 +0,800,2.9,2 +1,720,3.5,3 +0,620,3.58,2 +0,540,3.02,4 +0,480,3.43,2 +1,720,3.42,2 +0,580,3.29,4 +0,600,3.28,3 +0,380,3.38,2 +0,420,2.67,3 +1,800,3.53,1 +0,620,3.05,2 +1,660,3.49,2 +0,480,4,2 +0,500,2.86,4 +0,700,3.45,3 +0,440,2.76,2 +1,520,3.81,1 +1,680,2.96,3 +0,620,3.22,2 +0,540,3.04,1 +0,800,3.91,3 +0,680,3.34,2 +0,440,3.17,2 +0,680,3.64,3 +0,640,3.73,3 +0,660,3.31,4 +0,620,3.21,4 +1,520,4,2 +1,540,3.55,4 +1,740,3.52,4 +0,640,3.35,3 +1,520,3.3,2 +1,620,3.95,3 +0,520,3.51,2 +0,640,3.81,2 +0,680,3.11,2 +0,440,3.15,2 +1,520,3.19,3 +1,620,3.95,3 +1,520,3.9,3 +0,380,3.34,3 +0,560,3.24,4 +1,600,3.64,3 +1,680,3.46,2 +0,500,2.81,3 +1,640,3.95,2 +0,540,3.33,3 +1,680,3.67,2 +0,660,3.32,1 +0,520,3.12,2 +1,600,2.98,2 +0,460,3.77,3 +1,580,3.58,1 +1,680,3,4 +1,660,3.14,2 +0,660,3.94,2 +0,360,3.27,3 +0,660,3.45,4 +0,520,3.1,4 +1,440,3.39,2 +0,600,3.31,4 +1,800,3.22,1 +1,660,3.7,4 +0,800,3.15,4 +0,420,2.26,4 +1,620,3.45,2 +0,800,2.78,2 +0,680,3.7,2 +0,800,3.97,1 +0,480,2.55,1 +0,520,3.25,3 +0,560,3.16,1 +0,460,3.07,2 +0,540,3.5,2 +0,720,3.4,3 +0,640,3.3,2 +1,660,3.6,3 +1,400,3.15,2 +1,680,3.98,2 +0,220,2.83,3 +0,580,3.46,4 +1,540,3.17,1 +0,580,3.51,2 +0,540,3.13,2 +0,440,2.98,3 +0,560,4,3 +0,660,3.67,2 +0,660,3.77,3 +1,520,3.65,4 +0,540,3.46,4 +1,300,2.84,2 +1,340,3,2 +1,780,3.63,4 +1,480,3.71,4 +0,540,3.28,1 +0,460,3.14,3 +0,460,3.58,2 +0,500,3.01,4 +0,420,2.69,2 +0,520,2.7,3 +0,680,3.9,1 +0,680,3.31,2 +1,560,3.48,2 +0,580,3.34,2 +0,500,2.93,4 +0,740,4,3 +0,660,3.59,3 +0,420,2.96,1 +0,560,3.43,3 +1,460,3.64,3 +1,620,3.71,1 +0,520,3.15,3 +0,620,3.09,4 +0,540,3.2,1 +1,660,3.47,3 +0,500,3.23,4 +1,560,2.65,3 +0,500,3.95,4 +0,580,3.06,2 +0,520,3.35,3 +0,500,3.03,3 +0,600,3.35,2 +0,580,3.8,2 +0,400,3.36,2 +0,620,2.85,2 +1,780,4,2 +0,620,3.43,3 +1,580,3.12,3 +0,700,3.52,2 +1,540,3.78,2 +1,760,2.81,1 +0,700,3.27,2 +0,720,3.31,1 +1,560,3.69,3 +0,720,3.94,3 +1,520,4,1 +1,540,3.49,1 +0,680,3.14,2 +0,460,3.44,2 +1,560,3.36,1 +0,480,2.78,3 +0,460,2.93,3 +0,620,3.63,3 +0,580,4,1 +0,800,3.89,2 +1,540,3.77,2 +1,680,3.76,3 +1,680,2.42,1 +1,620,3.37,1 +0,560,3.78,2 +0,560,3.49,4 +0,620,3.63,2 +1,800,4,2 +0,640,3.12,3 +0,540,2.7,2 +0,700,3.65,2 +1,540,3.49,2 +0,540,3.51,2 +0,660,4,1 +1,480,2.62,2 +0,420,3.02,1 +1,740,3.86,2 +0,580,3.36,2 +0,640,3.17,2 +0,640,3.51,2 +1,800,3.05,2 +1,660,3.88,2 +1,600,3.38,3 +1,620,3.75,2 +1,460,3.99,3 +0,620,4,2 +0,560,3.04,3 +0,460,2.63,2 +0,700,3.65,2 +0,600,3.89,3 diff --git a/index.Rmd b/index.Rmd index e2ee379..a418a24 100644 --- a/index.Rmd +++ b/index.Rmd @@ -98,7 +98,8 @@ knitr::include_graphics("images/rbook1.png") - 第 \@ref(tidystats) 章介绍方差分析 - 第 \@ref(tests-as-linear) 章介绍统计检验与线性模型的等价性 - 第 \@ref(lmm) 章介绍多层线性模型 - - 第 \@ref(glm) 章介绍广义线性模型 + - 第 \@ref(poisson-regression) 章介绍广义线性模型中的泊松回归 + - 第 \@ref(logistic-regression) 章介绍logistic回归模型 - 第 \@ref(ordinal) 章介绍有序logistic回归模型 - 第 \@ref(tidymodels) 章介绍机器学习 - 应用篇 diff --git a/logistic_regression.Rmd b/logistic_regression.Rmd new file mode 100644 index 0000000..fa8de0d --- /dev/null +++ b/logistic_regression.Rmd @@ -0,0 +1,356 @@ +# logistic回归模型 {#logistic-regression} + +本章讲广义线性模型中的logistic回归模型。 + +## 问题 + +假定这里有一组数据，包含学生GRE成绩和被录取的状态(admit = 1，就是被录取；admit = 0，没有被录取) + +```{r, message=FALSE, warning=FALSE} +library(tidyverse) +gredata <- read_csv("demo_data/gredata.csv") +gredata +``` + +我们能用学生GRE的成绩**预测**录取状态吗？回答这个问题需要用到logistic回归模型， + + +$$ +\begin{align*} +\text{admit}_{i} &\sim \mathrm{Binomial}(1, p_{i}) \\ +logit(p_{i}) &= log\Big(\frac{p_{i}}{1 - p_{i}}\Big) = \alpha + \beta \cdot \text{gre}_{i} +\end{align*} +$$ + + + +这里 $p_i$ 就是被录取的概率。预测因子 $gre$ 的线性组合**模拟**的是 $log\Big(\frac{p_{i}}{1 - p_{i}}\Big)$ ，即对数比率(log-odds). + + +按照上面表达式，用**glm**函数写代码， + +```{r} +model_logit <- glm(admit ~ gre, + data = gredata, + family = binomial(link = "logit") +) + +summary(model_logit) +``` + +得到gre的系数是 + +```{r} +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） + +```{r eval=FALSE, include=FALSE} +library(effects) + +effect_link <- Effect("gre", mod = model_logit) + +plot(effect_link, + type = "link", + main = "gre effect plot\n(log odds scale)" +) +``` + +为了更清楚的理解，我们把log-odds的结果与gre成绩画出来看看 + +```{r} +logit_log_odds <- broom::augment_columns( + model_logit, + data = gredata, + type.predict = c("link") +) %>% + rename(log_odds = .fitted) +``` + + + +```{r} +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` + +```{r} +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成绩画出来看看 + +```{r} +logit_odds <- broom::augment_columns( + model_logit, + data = gredata, + type.predict = c("link") +) %>% + rename(log_odds = .fitted) %>% + mutate(odds_ratio = exp(log_odds)) +``` + + +```{r} +logit_odds %>% + ggplot(aes(x = gre, y = odds_ratio)) + + geom_line(color = "#FB9E07", size = 2) + + labs(title = "Odds (exp(β))", + 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的分数与`录取概率`就不可能是线性关系了，而是呈非线性关系。我们先看下非线性关系长什么样。 + +```{r eval=FALSE, include=FALSE} +library(effects) + +effect_response <- Effect("gre", mod = model) + +plot(effect_response, + type = "response", + main = "gre effect plot\n(probability scale)" +) +``` + + + +```{r} +logit_probs <- broom::augment_columns( + model_logit, + data = gredata, + type.predict = c("response") +) %>% + rename(pred_prob = .fitted) +``` + + + +```{r} +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； + + +```{block, type="danger"} +提请注意的是，以上三种度量的方式中： + +- `log_odds scale`，预测因子与`log_odds()`的关系是一个固定值，具有可加性 +- `odds scale`， 预测因子与`odds()`的关系是一个固定值，具有乘法性 +- `probability`，预测因子与`probability`的关系不再是一个固定值了 +``` + + + +用哪种度量方式来理解模型的输出，取决不同的场景。第一种方式容易计算但理解较为困难，第三种方式最容易理解，但不再是线性关系了。 + + +## 预测 + +### 预测与拟合 + +先认识下两个常用的函数`predict()`和`fitted()`. + + +```{r} +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) )$$ + + + +具体如下 + + +```{r} +gredata %>% + mutate( + pred = predict(model_logit), + fitted = fitted(model_logit), + pred2 = exp(pred) / (1 + exp(pred) ) + ) +``` + +我这样折腾无非是想让大家知道，在glm中`predict()`和`fit()`是不同的。 + +如果想让`predict()`也返回`录取概率`，也可以不用那么麻烦，事实上`predict`的`type = "response")` 选项已经为我们准备好了。 + +```{r} +gredata %>% + mutate( + pred = predict(model_logit, type = "response"), + fitted = fitted(model_logit) + ) +``` + +### 预测 + +有时候，我们需要对给定的GRE分数，用建立的模型**预测**被录取的概率 + +```{r} +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]` + + +```{r} +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) )` + + +```{r} +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]。 + + +```{r} +predict(model_logit, gredata, type = 'terms') %>% + as_tibble() %>% + head() +``` + + +我们复盘一次 + +```{r} +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() +``` + + + + + + +```{r, echo = F} +# remove the objects +# rm(list=ls()) +rm(gredata, logit_log_odds, logit_odds, logit_probs, model_logit, newdata, X) +``` + + +```{r, echo = F, message = F, warning = F, results = "hide"} +pacman::p_unload(pacman::p_loaded(), character.only = TRUE) +``` diff --git a/glm.Rmd b/poisson_regression.Rmd similarity index 99% rename from glm.Rmd rename to poisson_regression.Rmd index 037a156..1b58c02 100644 --- a/glm.Rmd +++ b/poisson_regression.Rmd @@ -1,4 +1,4 @@ -# 广义线性模型 {#glm} +# 广义线性模型 {#poisson-regression} ## 线性回归回顾 @@ -649,7 +649,7 @@ glm(number_of_fish ~ 1 + (1 | pollution_level), knitr::include_graphics(path = "images/One_Picture.png") ``` - +第 \@ref(logistic-regression) 章接着讲广义线性模型中的logistic回归模型。 diff --git a/rsconnect/documents/index.Rmd/bookdown.org/wangminjie/R4DS.dcf b/rsconnect/documents/index.Rmd/bookdown.org/wangminjie/R4DS.dcf index e41fa68..a6bf654 100644 --- a/rsconnect/documents/index.Rmd/bookdown.org/wangminjie/R4DS.dcf +++ b/rsconnect/documents/index.Rmd/bookdown.org/wangminjie/R4DS.dcf @@ -5,6 +5,6 @@ account: wangminjie server: bookdown.org hostUrl: https://bookdown.org/__api__ appId: 3039 -bundleId: 29178 +bundleId: 29384 url: https://bookdown.org/wangminjie/R4DS/ -when: 1602040192.27007 +when: 1602488585.35307