【R模型】R语言二元logistic回归 (保姆级教程)

• 💂 个人信息：酷在前行
• 👍 版权: 博文由【酷在前行】原创、需要转载请联系博主
• 👀 如果博文对您有帮助，欢迎点赞、关注、收藏 + 订阅专栏
• 🔖 本文收录于【R模型】，该专栏主要介绍R语言各类型机器学习，如线性回归模型、广义线性模型、混合线性模型、随机森林模型、支持向量机模型、决策树模型等。请大家多多关注点赞和支持，共同进步~ 欢迎大家订阅！

   在任何情况下，如果我们的结果是二进制的，我们实际上是在处理可能性。这些可能性通常不是线性的，因此我们不再有输入与输出结果直接呈线性相关。直接的线性回归方法，如普通最小二乘法回归，并不适合于这种类型数据的分析。相反，线性关系可以通过结果变量的转换来推断，这给我们提供了一条建立可解释模型的途径。因此，二项式逻辑回归被说成是属于广义线性模型或GLM的一个类别。理解逻辑回归并在实践中可靠地使用它并不简单，但在人员分析领域，它是一项宝贵的技能。

🐣 一、二项逻辑回归的介绍

逻辑函数最初是由比利时数学家Pierre François Verhulst在19世纪中期引入的，作为人类、动物和某些植物和水果物种的人口增长建模工具。那时，人们普遍认为，人口增长不可能永远呈指数级增长，环境和资源的限制对人口的规模施加了最大的限制。Verhulst函数的公式为:

其中e是指数常数，x0是x在中点的值，L是y的最大值(称为“承载力”)，k是曲线的最大斜率。

   如图1所示，logistic函数可以准确地捕捉到人口增长的理论阶段，在初始阶段增长缓慢，在中间阶段转向指数增长，然后随着人口接近其承载能力而转向缓慢增长。

 图1 Verhulst的logistic函数同时模拟了人口增长的指数特性和自然极限

20世纪初，从经济学和化学的应用开始，logistic函数作为一种有用的现象建模工具被广泛应用于许多领域。在统计学中，我们可以观察到，logistic函数与概率的累积正态分布具有类似的s型(或sigmoid)，如图2 所示，其中x刻度表示均值附近的标准差。logistic函数给出了一个数学模型，其中的系数很容易根据结果的可能性进行解释。因此，不出所料，logistic模型很快成为建模概率现象的常用方法。

 图2 logistic函数(蓝色虚线)与累积正态分布(红色实线)非常相似

🐤 二、二项逻辑回归的R语言实例

当感兴趣的结果本质上是二值或二值的时候，可以使用二项逻辑回归。也就是说，它需要两个值中的一个。例如，1或0，真或假，是或否。这些类通常被描述为“正类”和“负类”。这里有一个潜在的假设，即结果的累积概率形状类似于累积正态分布。下面是一些可以使用二项式逻辑回归解决的问题。

• 给定一组关于组织中销售经理的数据，包括与目标相比的绩效、团队规模、在组织中的任期和其他因素，这些因素对个人获得高绩效评级的可能性有什么影响?
• 给定一组人口统计、收入和位置数据，它们各自对个人在选举中投票的可能性有什么影响?
• 给定一组关于足球运动员在游戏中的活动的统计数据，每个统计数据与球员进球的可能性有什么关系?

案例 : 你是一家大公司的分析师，由全国各地的区域销售团队组成。这家公司每年两次提拔一些销售人员。晋升由每个区域销售团队的负责人自行决定，考虑财务业绩、客户满意度评级、近期业绩评级和个人判断。公司管理层要求您进行分析，以确定财务业绩、客户评级和业绩评级等因素如何影响给定销售人员获得晋升的可能性。你得到了一个数据集，其中包含最近三年考虑晋升的销售人员的数据。销售员数据集包含以下字段。

• promoted : 二进制值，表示晋升为1，未晋升为0
• sales : 在促销期间属于个人的销售额(千元)
• customer_rate : 在促销期间对个人客户的平均满意度调查
• performance : 晋升前最近的业绩评级，从1(最低)到4(最高)

  🍎 1. 数据读取及探索

url <- "http://peopleanalytics-regression-book.org/data/salespeople.csv"
salespeople <- read.csv(url)
head(salespeople)
##   promoted sales customer_rate performance
## 1        0   594          3.94           2
## 2        0   446          4.06           3
## 3        1   674          3.83           4
## 4        0   525          3.62           2
## 5        1   657          4.40           3
## 6        1   918          4.54           2

从前六行数据来看与预期符合，查看下数据摘要

summary(salespeople)
##     promoted          sales       customer_rate    performance 
##  Min.   :0.0000   Min.   :151.0   Min.   :1.000   Min.   :1.0  
##  1st Qu.:0.0000   1st Qu.:389.2   1st Qu.:3.000   1st Qu.:2.0  
##  Median :0.0000   Median :475.0   Median :3.620   Median :3.0  
##  Mean   :0.3219   Mean   :527.0   Mean   :3.608   Mean   :2.5  
##  3rd Qu.:1.0000   3rd Qu.:667.2   3rd Qu.:4.290   3rd Qu.:3.0  
##  Max.   :1.0000   Max.   :945.0   Max.   :5.000   Max.   :4.0  
##                   NA's   :1       NA's   :1       NA's   :1

首先，我们看到少量缺失值，我们应该删除这些观测值。我们看到大约三分之一的人得到了提升，销售额从151k美元到945k美元，正如预期的那样，平均满意度评级从1到5之间，最后我们看到了四个绩效评级，尽管绩效类别应该是有序因子时是数字的，而升职是数字的，而应该是类别的。我们先对这些数据进行转换，然后绘制一个配对图来快速了解一些潜在的关系，如图3 所示。

library(GGally)
# 剔除NA值
salespeople <- salespeople[complete.cases(salespeople), ]
salespeople$performance <- ordered(salespeople$performance, 
                                   levels = 1:4)
salespeople$promoted <- as.factor(salespeople$promoted)
ggpairs(salespeople)

 图3 salespeople数据集的配对图

从这幅图中我们可以看到，升职的人的销售额明显高于没有升职的人。我们还看到了客户评分和销售额之间的中度关系，这很直观(如果客户不怎么看好你，销售额可能不会很高)。

因此，我们可以看到这里可能存在一些与我们的结果的关系，但不清楚如何梳理它们并量化它们彼此之间的关系。让我们来看看二项逻辑回归是如何帮助我们实现这一目标的。

  🍎 2. 构建逻辑回归模型

假设你有一个结果事件y，它要么发生，要么不发生。显然，y发生的概率，即P (y = 1)的值在0到1之间。现在假设某个输入变量x对事件发生的概率有正向影响。那么你自然会期望P (y = 1)随着x的增加而增加。在我们的销售人员数据集中，让我们根据销售输入绘制促销结果。这可以在图4中看到。

 图4 salespeople数据集中促销与销售额的关系图

很明显，销售水平越高，促销的可能性越大。当我们沿着x轴从左到右移动，并逐渐包含越来越多的销售额较高的个人时，我们知道整体上推广的概率是逐渐增加的。

在此数据上，我们将x0设置为销售的均值，k设置为某个最大梯度值。在图5中，我们可以看到不同k值下的这些logistic函数。所有这些似乎都在某种程度上反映了我们正在观察的模式，但如何确定最佳拟合的逻辑函数呢?

 图5 在前一幅图上叠加了不同梯度的逻辑函数

  🍎 3. 对数优势比（odds）的概念

   我们来仔细看看指数常数e在logistic函数分母上的下标。注意，因为x0是一个常数，我们有

其中β0 = - kx0和β1 = k。因此,

这个方程很直观。随着x的增大，e−(β0 + β1x)越来越小，趋近于0，因此P (y = 1)接近其理论最大值1。当x的值趋于0时，我们看到P(y = 1)的值趋于最小值 $\frac{1}{1+e^{-\beta 0}}$ 。回到销售人员的例子，我们可以看到β0有助于确定在根本没有销售的情况下促销的基线概率。如果β0具有一个非常负的值，则基线概率将接近其理论最小值0。

让我们将β0和β1在积极结果可能性中的作用正式化。我们知道，对于任何二元事件y, P (y = 0)等于1−P (y = 1)，因此

把这些放在一起，我们发现

或者，如果两边同时取自然对数

P (y = 1)是事件发生的概率，P (y = 0)是事件不发生的概率。你可能在赛马或其他赌博运动中很熟悉，这两者的比率代表一个事件的几率。例如，如果某匹马的赔率为1:4，这意味着他们获胜的概率为20%，不获胜的概率为80%。

因此，我们可以得出结论，y的优势比的自然对数(通常称为y的对数优势比)在x中是线性的。

  🍎 4. 解释优势比及模型结果

sales_model <- glm(formula = promoted ~ sales, 
                   data = salespeople, family = "binomial")
sales_model$coefficients
##  (Intercept)        sales 
## -21.77642020   0.03675848

我们可以将系数解释为:

1. Intercept系数是x输入值为零时的对数优势比值——它是没有销售时促销的对数优势比值。
2. sales系数表示每增加一单位销售额，促销的对数几率的增加。

我们可以通过指数函数将这些系数从对数优势系数转换为优势系数，从而恢复到我们之前得到的恒等式

由此，我们可以理解为$e^{\beta 1}$代表在没有销售的情况下促销的基本几率，每增加一个单位的销售，这些基本几率乘以$e^{\beta 1}$。鉴于$e^{\beta 1}$对概率的这种乘法效应，它被称为 优势比(odds ratio) 。

# 将对数比值转换为基数比值和比值比
exp(sales_model$coefficients)
##  (Intercept)        sales 
## 3.488357e-10 1.037442e+00

因此，我们可以看到，销售为零的促销的基本几率非常接近于零，这是有道理的。请注意，只有在不可能属于正类(即没有人得到提升)的情况下，概率才能精确为零。我们还可以看到，每卖出一个单位(即每1000美元)，基本几率就会增加大约1.04倍。换句话说，促销几率增加了4%。

  🍎 5. 优势比与概率

有必要花点时间来理解优势比的概念以及它与概率的关系。这两个术语经常被当作同义词使用，这在解释逻辑回归模型时可能会导致严重的误解。

如果某个事件的概率是0.1，那么这意味着它的概率是1:9，即0.111。如果概率是0.5，那么优势比就是1; 如果概率是0.9，那么优势比就是9; 如果概率是0.99，那么优势比就是99。当概率接近1时，优势比会呈指数级增长，如图6所示:

 图6 优势比与概率的关系图

这样做的结果是，优势比的给定增加可能会对概率产生不同的影响，这取决于最初的概率是什么。如果概率已经很低了，例如0.1，那么4%的优势比增加将转化为0.116的优势比，这将转化为0.103586的新概率，代表3.59%的优势比增加，这非常接近于优势比的增加。如果该概率已经很高，例如0.9，那么4%的优势比增加将转化为9.36的优势比，这将转化为0.903475的新概率，代表0.39%的优势比增加，这与优势比的增加非常不同。图7 显示了根据事件的原始概率增加4%的概率所产生的影响。

 图7 优势比增加4%对原始概率的影响

我们可以看到，基本概率越接近零，增加对优势比和概率的影响就越相似。然而，事件发生的概率越高，优势比增加的影响就越小。无论如何，记住将优势比转换为概率的公式和将优势比转换为概率的公式是有用的。如果O表示概率，P表示概率，则有:

  🍎 6. 多元二项逻辑回归模型的介绍

设y为二分类结果，设x1, x1，…，xp为输入变量。然后

对于系数β0， β1，…，βp。像以前一样:

• β0表示所有输入为0时结果的对数概率
• 假设其他输入不变，每个βi代表xi单位变化时，结果的对数概率的增加。

我们得到，

因此我们可以得出结论:

• eβ0表示所有输入均为0时的结果概率。
• 每个eβi代表在其他输入不变的情况下，xi增加一单位所对应的优势比(即x i增加一单位，我们的结果的优势比乘以eβi)。

  🍎 7. 多元二项逻辑回归模型解释

library(makedummies)
# 转换虚拟变量
salespeople_dummies <- makedummies::makedummies(salespeople)
# 查看前六行
head(salespeople_dummies)
##   promoted sales customer_rate performance_2 performance_3 performance_4
## 1        0   594          3.94             1             0             0
## 2        0   446          4.06             0             1             0
## 3        1   674          3.83             0             0             1
## 4        0   525          3.62             1             0             0
## 5        1   657          4.40             0             1             0
## 6        1   918          4.54             1             0             0

现在我们可以运行我们的模型( 公式promoted ~ .表示逆序回归)并查看系数。

# 运行广义线性模型 选择二分类（逻辑）
full_model <- glm(formula = "promoted ~ .",
                  family = "binomial",
                  data = salespeople_dummies)
# 获得系数摘要
(coefs <- summary(full_model)$coefficients)
##                   Estimate  Std. Error    z value     Pr(>|z|)
## (Intercept)   -19.85893195 3.444078811 -5.7661085 8.112287e-09
## sales           0.04012425 0.006576429  6.1012212 1.052611e-09
## customer_rate  -1.11213130 0.482681585 -2.3040682 2.121881e-02
## performance_2   0.26299953 1.021980179  0.2573431 7.969139e-01
## performance_3   0.68495453 0.982166998  0.6973911 4.855581e-01
## performance_4   0.73449340 1.071963758  0.6851849 4.932272e-01

请注意，只有三个性能显示。这是因为每个人都属于四个性能类别之一，因此模型使用performance_1作为参考案例。我们可以将每个性能系数解释为从performance_1移动到该性能类别的影响。

从系数汇总的最后一列(系数p值)可以看出，只有sales和customer_rate满足小于0.05的显著性阈值。有趣的是，从估计列中可以看出customer_rate对促销的对数概率有负面影响。为方便起见，我们可以在系数汇总中增加一列来创建估计系数的指数，以便查看优势比。如果有需要，也可以删除不太有用的列。

# 创建包含估计值、p值和优势比的系数表
(full_coefs <- cbind(coefs[ ,c("Estimate", "Pr(>|z|)")], 
                     odds_ratio = exp(full_model$coefficients))) 
##                   Estimate     Pr(>|z|)   odds_ratio
## (Intercept)   -19.85893195 8.112287e-09 2.373425e-09
## sales           0.04012425 1.052611e-09 1.040940e+00
## customer_rate  -1.11213130 2.121881e-02 3.288573e-01
## performance_2   0.26299953 7.969139e-01 1.300826e+00
## performance_3   0.68495453 4.855581e-01 1.983682e+00
## performance_4   0.73449340 4.932272e-01 2.084426e+00

现在我们可以这样解释我们的模型:

• 在其他条件相同的情况下，销售额对促销的可能性有显著的积极影响，每增加1000美元的销售额，促销的可能性就增加4%
• 在其他条件相同的情况下，客户评分对促销的可能性有显著的负面影响，一个完整的评分越高，促销的可能性就越低67%
• 在其他条件相同的情况下，绩效评分对晋升的可能性没有显著影响

第二个结论可能有悖直觉，销售和客户评级之间已经存在中度相关性，而这个模型将控制这种关系。我们的优势比假设所有其他变量都相同。因此，如果两个人具有相同的销售和业绩评级，那么客户评级较低的那个人更有可能获得晋升。同样，如果两个人有相同的销售水平和相同的客户评级，他们的业绩评级将对晋升的可能性没有显著影响。

此外，我们可以计算指数置信区间来为优势比提供一个范围。

exp(confint(full_model))
##                      2.5 %       97.5 %
## (Intercept)   7.879943e-13 7.385387e-07
## sales         1.029762e+00 1.057214e+00
## customer_rate 1.141645e-01 7.793018e-01
## performance_2 1.800447e-01 1.061602e+01
## performance_3 3.060299e-01 1.547188e+01
## performance_4 2.614852e-01 1.870827e+01

因此，我们可以说，在其他条件相同的情况下，每增加一个销售单位，促销的几率就会增加3.0%到5.7%，客户评分每增加一个点，促销的几率就会减少22%到89%。

与其他回归模型类似，在解释时需要考虑单位尺度。优势比下降89%似乎比几率上升5.7%要重要得多。然而，高达5.7%的增长是针对数千个销售单元中的一个单元(1000美元)，超过10或100个额外单元可以对促销几率产生显著的复合效应。在只有4个满分的范围内，减少了多达89%的全客户评分点。

  🍎 8. 拟合优度（ R 2 R^{2} R2)

simpler_model <- glm(formula = promoted ~ sales + customer_rate,
                     family = "binomial",
                     data = salespeople)

展示了一条“扭曲”的3D sigmoid曲线，反映了销售额和客户率对结果的相对影响。

 图8 simpler_model拟合结果的3D可视化

查看模型摘要:

summary(simpler_model)
## 
## Call:
## glm(formula = promoted ~ sales + customer_rate, family = "binomial", 
##     data = salespeople)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.02984  -0.09256  -0.02070   0.00874   3.06380  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -19.517689   3.346762  -5.832 5.48e-09 ***
## sales           0.040389   0.006525   6.190 6.03e-10 ***
## customer_rate  -1.122064   0.466958  -2.403   0.0163 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 440.303  on 349  degrees of freedom
## Residual deviance:  65.131  on 347  degrees of freedom
## AIC: 71.131
## 
## Number of Fisher Scoring iterations: 8

我们的摘要没有提供整体模型拟合或拟合优度的统计量。主要原因是在统计学界，对于logistic回归的模型拟合没有一个统一的衡量标准。然而，有很多方法可以用来评估这些模型的拟合和拟合优度。

library(DescTools)
DescTools::PseudoR2(
  simpler_model, 
  which = c("McFadden", "CoxSnell", "Nagelkerke", "Tjur")
)
##   McFadden   CoxSnell Nagelkerke       Tjur 
##  0.8520759  0.6576490  0.9187858  0.8784834

我们看到Cox和Snell变体明显低于其他估计。然而，其他估计是合理对齐的，表明拟合性很强。

logistic回归模型的拟合优度检验将预测结果与观测结果进行比较，并检验它们相似的零假设。这意味着，与线性回归不同，p值越低，拟合越差。一种常用的方法是Hosmer-Lemeshow检验(Hosmer-Lemeshow test)，它根据拟合概率将观测值分成若干组(通常是10组)，计算每组中为正的比例，然后使用卡方检验将其与模型预测的期望比例进行比较。然而，这种方法有局限性。在样本量较小的情况下，根据所使用的组的数量，它可能返回高度不同的结果，这是个重要的问题。因此，我们建议使用一系列的拟合优度测试，不要完全依赖于任何一种特定的方法。

在R中，LogisticDx包为逻辑回归模型提供了一系列诊断工具，推荐用于探索。下面是一个使用gof()函数评估拟合优度的例子。

library(LogisticDx)
simpler_model_diagnostics <- LogisticDx::gof(simpler_model, 
                                             plotROC = FALSE)
names(simpler_model_diagnostics)
## [1] "ct"    "chiSq" "ctHL"  "gof"   "R2"    "auc"

这个列表中的gof对象提供了一系列拟合优度统计量。

simpler_model_diagnostics$gof
##          test  stat         val df        pVal
## 1:         HL chiSq  3.44576058  8 0.903358158
## 2:        mHL     F  2.74709957  8 0.005971045
## 3:       OsRo     Z -0.02415249 NA 0.980730971
## 4: SstPgeq0.5     Z  0.88656856 NA 0.375311227
## 5:   SstPl0.5     Z  0.96819352 NA 0.332947728
## 6:    SstBoth chiSq  1.72340251  2 0.422442787
## 7: SllPgeq0.5 chiSq  1.85473814  1 0.173233325
## 8:   SllPl0.5 chiSq  0.68570870  1 0.407627859
## 9:    SllBoth chiSq  1.86640617  2 0.393291943

这证实了几乎所有的测试，包括列表中的第一个Hosmer-Lemeshow测试，都表明我们的模型适合。

二项逻辑回归模型也可用于预测环境下的准确度评估，如精确率、召回率和roc曲线分析。这些模型特别适用于将逻辑回归模型作为机器学习的预测分类器来实现。

  🍎 9. 模型简约（AIC）

我们看到，在线性回归和逻辑回归方法中，当确定变量对结果没有显著影响时，我们决定将其从模型中删除。在其他条件相同的情况下，最简单的解释是最好的。从这个意义上说，一个包含对其主要推理目标没有贡献的信息的模型比它需要的更复杂。这样的模型在向其他人解释其结果时增加了沟通负担，而没有显著的分析好处。

简约指的是小心对待资源或信息的概念。如果一个模型可以用更少的输入数量达到相同(或非常接近相同)的拟合，那么它可以被描述为更优。赤池信息准则(Akaike Information Criterion，简称AIC)是一种衡量模型简约性的标准，适用于对数似然模型(如logistic回归模型)，AIC越低表明模型越简约。AIC常作为逻辑回归模型汇总报告的标准计算，但也可以独立计算。我们使用AIC来比较本章中模型的不同迭代次数。

AIC(sales_model)
## [1] 76.49508

AIC(simpler_model)
## [1] 71.13145

AIC(full_model)
AIC(full_model)

我们可以看到，根据AIC，受限于我们的两个重要输入(销售和客户评级)的模型被确定为最精简的模型。注意，AIC不应该用来解释模型的质量或可信度——即使AIC最小，也可能仍然非常不合适。

当有大量输入变量时，模型精简成为一个重要问题。一般来说，模型中的输入变量越多，就越有可能难以清楚地解释该模型，也就越有可能出现测量问题，例如多重共线性。急于取悦客户、客户、教授或老板的分析师很容易被诱惑，为他们的模型想出新的潜在输入，这些输入通常是从模型中已经输入的指标中数学推导出来的。过不了多久，模型就会变得过于复杂，在极端情况下，输入的数量比观测的数量还要多。管理模型复杂性的主要方法是谨慎选择模型输入。当无法避免大量输入时，系数正则化方法(如LASSO回归)可以帮助模型简化。

  🍎 10. 注意事项

要对新数据进行预测，使用predict()函数。该函数识别所使用的模型类型(广义线性模型)，并相应地调整其预测方法。特别是，如果你想返回新观测被提升的概率，则需要使用type = "response"作为参数。

(new_data <- data.frame(sales = c(420, 510, 710), 
                        customer_rate = c(3.4, 2.3, 4.2)))
                        ##   sales customer_rate
## 1   420           3.4
## 2   510           2.3
## 3   710           4.2

预测晋升概率

predict(simpler_model, new_data, type = "response")
##          1          2          3 
## 0.00171007 0.18238565 0.98840506

线性回归的许多原则在逻辑回归中也同样重要。例如，应该以类似的方式管理输入变量。应该关注共线性和多重共线性。输入变量的交互可以建模。在大多数情况下，分析人员在考虑这些问题时应该了解逻辑回归模型中发生的基本数学变换。例如，在线性回归中，系数对y有直接的加性影响，而在逻辑回归中，系数对y的对数优势比有直接的加性影响，或者它们的指数对y的优势比有直接的乘法影响。因此，系数过高估计(如未通过共线性时可能发生的情况)可能导致推论严重夸大输入变量的重要性或效果。

由于我们的结果变量是二元的，logistic回归模型的残差限制了其对所研究问题的直接应用。在实际环境中，logistic回归模型的残差很少被检验，但它们在识别离群点或特别有影响力的观测值以及评估拟合优度方面很有用。在对残差进行检测时，需要对残差进行变换，以便进行适当的分析。例如，皮尔逊残差是逻辑回归残差的标准化形式，可以期望在足够大的样本上具有正态分布。在图9中，我们可以看到simpler_model的情况就是如此，但模型中也有少量数据被严重低估。

d <- density(residuals(simpler_model, "pearson"))
plot(d, main= "")

 图9 simpler_model的皮尔逊残差分布文章来源地址https://www.toymoban.com/news/detail-439263.html

到了这里，关于【R模型】R语言二元logistic回归 (保姆级教程)的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！