本文主要内容：
- 生成Logistic 回归模型结果
- 绘制Logistic回归曲线
- 绘制带有数据分布的Logistic回归曲线
当你拟合逻辑回归模型时，有很多方法可以显示结果。最为传统的方法是在表格中展示系数的摘要。但是由于逻辑回归曲线的弯曲性质，通过绘图显示拟合线，通常是一种更好的展示方式。典型的逻辑回归线图的主要缺点是它们通常不显示数据，因为会出现在 y 轴上过度绘图的现象。但是ggdist 包支持在绘制逻辑回归曲线时可以轻松地显示数据。
阅读这篇文章，我假设你了解一些背景知识：
- 熟悉Logistic回归模型的基本原理和结果分析。
- 熟悉R语言的基本操作。数据整理和绘图是在
tidyverse包 和broom包的帮助下完成的。数据来自包fivethirtyeight。
1. 导入所需包和数据
- 加载包
if(!require(pacman))install.packages("pacman")
pacman::p_load(
rio, # to import data
here, # to locate files
tidyverse, # to clean, handle, and plot the data (includes ggplot2 package)
fivethirtyeight,
broom,
ggdist,
showtext # 支持ggplot显示中文
)
#Set the default theme for ggplot objects to theme_bw()
theme_set(theme_bw())
theme_update(panel.grid = element_blank())- 数据集
在这篇文章中，我们将使用 bechdel 数据集。从文档中，该数据是‘反对好莱坞排斥女性的美元和美分案’故事背后的原始数据。
data("bechdel")
skimr::skim(bechdel)Name | bechdel |
Number of rows | 1794 |
Number of columns | 15 |
_______________________ | |
Column type frequency: | |
character | 5 |
factor | 1 |
numeric | 9 |
________________________ | |
Group variables | None |
Variable type: character
skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
imdb | 0 | 1 | 8 | 10 | 0 | 1794 | 0 |
title | 0 | 1 | 1 | 83 | 0 | 1768 | 0 |
test | 0 | 1 | 2 | 16 | 0 | 10 | 0 |
binary | 0 | 1 | 4 | 4 | 0 | 2 | 0 |
code | 0 | 1 | 8 | 8 | 0 | 85 | 0 |
Variable type: factor
skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
clean_test | 0 | 1 | TRUE | 5 | ok: 803, not: 514, men: 194, dub: 142 |
Variable type: numeric
skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
year | 0 | 1.00 | 2002.55 | 8.98 | 1970 | 1998 | 2005 | 2009 | 2013 | ▁▁▂▅▇ |
budget | 0 | 1.00 | 44826462.61 | 48186026.12 | 7000 | 12000000 | 28000000 | 60000000 | 425000000 | ▇▁▁▁▁ |
domgross | 17 | 0.99 | 69132048.28 | 80367309.51 | 0 | 16311571 | 42194060 | 93354918 | 760507625 | ▇▁▁▁▁ |
intgross | 11 | 0.99 | 150385700.05 | 210335267.83 | 828 | 26129470 | 76482461 | 189850852 | 2783918982 | ▇▁▁▁▁ |
budget_2013 | 0 | 1.00 | 55464608.45 | 54918635.60 | 8632 | 16068918 | 36995786 | 78337905 | 461435929 | ▇▂▁▁▁ |
domgross_2013 | 18 | 0.99 | 95174783.58 | 125965348.89 | 899 | 20546594 | 55993640 | 121678352 | 1771682790 | ▇▁▁▁▁ |
intgross_2013 | 11 | 0.99 | 197837984.97 | 283507948.20 | 899 | 33232604 | 96239640 | 241478970 | 3171930973 | ▇▁▁▁▁ |
period_code | 179 | 0.90 | 2.42 | 1.19 | 1 | 1 | 2 | 3 | 5 | ▇▇▆▅▂ |
decade_code | 179 | 0.90 | 1.94 | 0.69 | 1 | 1 | 2 | 2 | 3 | ▅▁▇▁▃ |
我们的因变量是二进制的，它指示给定电影是否通过了 Bechdel 测试。在数据集中的 1,794 部电影中，只有不到一半通过。
bechdel %>%
count(binary) %>%
mutate(percent = 100 * n/ sum(n))# A tibble: 2 × 3
binary n percent
<chr> <int> <dbl>
1 FAIL 991 55.2
2 PASS 803 44.8我们唯一的预测变量是budget_2013，每部电影2013年的的预算（美元）。
bechdel %>%
ggplot(aes(x = budget_2013)) +
geom_histogram() +
facet_wrap(~ binary, ncol = 1)
为了方便，我们将字符变量 binary 转换为一个传统的 0/1 数值变量，称为 pass。
bechdel = bechdel %>%
mutate(pass = ifelse(binary == "FAIL", 0, 1))
bechdel %>%
select(binary, pass) %>%
head()# A tibble: 6 × 2
binary pass
<chr> <dbl>
1 FAIL 0
2 PASS 1
3 FAIL 0
4 FAIL 0
5 FAIL 0
6 FAIL 02. 建模
我们的Logistic统计模型公式如下：
我们使用传统的 logit 连接函数来确保二项式概率被限制在0和1的范围内。我们可以像这样使用基本的 R glm() 函数来拟合这样的模型。
fit = glm(
data = bechdel,
family = binomial,
pass ~ 1 + budget_2013
)传统方式是在系数表中呈现结果，见下表。
tidy(fit) %>%
knitr::kable()term | estimate | std.error | statistic | p.value |
(Intercept) | 0.1113148 | 0.0689661 | 1.614051 | 0.1065163 |
budget_2013 | 0.0000000 | 0.0000000 | -6.249724 | 0.0000000 |
由于budget_2013 变量的规模，它的点估计和标准误差都非常小。
为了更清晰的展示自变量的影响效应，下面是预算增加 100,000,000 美元时预期的对数几率减少。
c(coef(fit)[2], confint(fit)[2,]) * 1e8budget_2013 2.5 % 97.5 %
-0.5972374 -0.7875178 -0.4126709请注意如何添加 95% 置信区间以获得良好的衡量标准。
3. Logistic回归曲线
现在已经用表格的方式解释了模型。为了更好的理解和解释模型结果，我们将使用广泛使用的仅显示拟合线的方法。
# define the new data
nd <- tibble(budget_2013 = seq(from = 0, to = 500000000, length.out = 100))
p <-
# 计算拟合线和 标准差SE 的
predict(fit,
newdata = nd,
type = "link",
se.fit = TRUE) %>%
data.frame() %>%
mutate(ll = fit - 1.96 * se.fit,
ul = fit + 1.96 * se.fit) %>%
select(-residual.scale, -se.fit) %>%
mutate_all(plogis) %>%
bind_cols(nd)
glimpse(p)Rows: 100
Columns: 4
$ fit <dbl> 0.5278000, 0.5202767, 0.5127442, 0.5052059, 0.4976652, 0.4…
$ ll <dbl> 0.4940356, 0.4881515, 0.4821772, 0.4760998, 0.4699043, 0.4…
$ ul <dbl> 0.5613120, 0.5522351, 0.5432161, 0.5342767, 0.5254405, 0.5…
$ budget_2013 <dbl> 0, 5050505, 10101010, 15151515, 20202020, 25252525, 303030…- 逻辑回归模型的常规线图。
p %>%
ggplot(aes(x = budget_2013, y = fit)) +
geom_ribbon(aes(ymin = ll, ymax = ul),
alpha = 1/2) +
geom_line() +
scale_y_continuous("probability of passing",
expand = c(0, 0), limits = 0:1)
拟合线为黑色和 95% 置信区间的半透明灰色丝带标记。该曲线很好地展示了预算较大的电影在通过贝克德尔测试时往往做得更差。
3. 通过添加数据改进可视化
如果想将数据添加到我们的绘图中，最简单的方法是使用 geom_point()函数。
p %>%
ggplot(aes(x = budget_2013, y = fit)) +
geom_ribbon(aes(ymin = ll, ymax = ul),
alpha = 1/2) +
geom_line() +
geom_point(data = bechdel,
aes(y = pass),
alpha = 1/2) +
scale_y_continuous("probability of passing",
expand = expansion(mult = 0.01))
即使通过使用 alpha 参数使点半透明，过度绘图问题也使得很难理解数据。
- 另一种方法是添加垂直抖动。
我们可以用 geom_jitter() 来做到这一点。
p %>%
ggplot(aes(x = budget_2013, y = fit)) +
geom_ribbon(aes(ymin = ll, ymax = ul),
alpha = 1/2) +
geom_line() +
geom_jitter(data = bechdel,
aes(y = pass),
size = 1/4, alpha = 1/2, height = 0.05) +
scale_y_continuous("probability of passing",
expand = c(0, 0))
尽管有了很大的改进，但这种方法仍然不能最好地描述budget_2013 值的分布。
- 添加自变量的分布
如果可能的话，最好用更像直方图的方式明确描述 budget_2013 分布。
p %>%
ggplot(aes(x = budget_2013)) +
geom_ribbon(aes(ymin = ll, ymax = ul),
alpha = 1/2) +
geom_line(aes(y = fit)) +
stat_dots(data = bechdel,
aes(y = pass, side = ifelse(pass == 0, "top", "bottom")),
scale = 1/3) +
scale_y_continuous("probability of passing",
expand = c(0, 0))
使用 stat_dots() 函数，我们添加了点图，这是直方图的绝妙替代品，直方图将每个数据值显示为一个单独的点。对于 side 参数，我们使用条件语句告诉 stat_dots() 我们希望将budget_2013 的一些值显示在底部，而将其他值显示在顶部。使用 scale 参数，我们指出了我们希望点图分布占据 y 轴范围内的总空间的多少。
- 更精致的版本
p %>%
ggplot(aes(x = budget_2013)) +
geom_ribbon(aes(ymin = ll, ymax = ul),
alpha = 1/2) +
geom_line(aes(y = fit)) +
stat_dots(data = bechdel %>%
mutate(binary = factor(binary, levels = c("PASS", "FAIL"))),
aes(y = pass,
side = ifelse(pass == 0, "top", "bottom"),
color = binary),
scale = 0.4, shape = 19) +
scale_color_manual("Bechdel test", values = c("#009E73", "#D55E00")) +
scale_x_continuous("budget (in 2013 dollars)",
breaks = c(0, 1e8, 2e8, 3e8, 4e8),
labels = c(0, str_c(1:4 * 100, " mill")),
expand = c(0, 0), limits = c(0, 48e7)) +
scale_y_continuous("probability of passing",
expand = c(0, 0)) +
theme(panel.grid = element_blank())
- 添加自变量的直方图分布
其他分布形式也是可能的。例如，这里我们在 stat_slab() 函数中设置slab_type = "histogram" 以将点图换成直方图。
p %>%
ggplot(aes(x = budget_2013)) +
geom_ribbon(aes(ymin = ll, ymax = ul),
alpha = 1/2) +
geom_line(aes(y = fit)) +
# the magic lives here
stat_slab(data = bechdel %>%
mutate(binary = factor(binary, levels = c("PASS", "FAIL"))),
aes(y = pass,
side = ifelse(pass == 0, "top", "bottom"),
fill = binary, color = binary),
slab_type = "histogram",
scale = 0.4, breaks = 40, size = 1/2) +
scale_fill_manual("Bechdel test", values = c(alpha("#009E73", .7), alpha("#D55E00", .7))) +
scale_color_manual("Bechdel test", values = c("#009E73", "#D55E00")) +
scale_x_continuous("budget (in 2013 dollars)",
breaks = c(0, 1e8, 2e8, 3e8, 4e8),
labels = c(0, str_c(1:4 * 100, " mill")),
expand = c(0, 0), limits = c(0, 48e7)) +
scale_y_continuous("probability of passing",
expand = c(0, 0)) +
theme(panel.grid = element_blank())
