如何在R语言中使用Logistic回归模型

在日常学习或工作中经常会使用线性回归模型对某一事物进行预测，例如预测房价、身高、GDP、学生成绩等，发现这些被预测的变量都属于连续型变量。然而有些情况下，被预测变量可能是二元变量，即成功或失败、流失或不流失、涨或跌等，对于这类问题，线性回归将束手无策。这个时候就需要另一种回归方法进行预测，即Logistic回归。

在实际应用中，Logistic模型主要有三大用途：

1）寻找危险因素，找到某些影响因变量的"坏因素"，一般可以通过优势比发现危险因素；

2）用于预测，可以预测某种情况发生的概率或可能性大小；

3）用于判别，判断某个新样本所属的类别。

Logistic模型实际上是一种回归模型，但这种模型又与普通的线性回归模型又有一定的区别：

1）Logistic回归模型的因变量为二分类变量；

2）该模型的因变量和自变量之间不存在线性关系；

3）一般线性回归模型中需要假设独立同分布、方差齐性等，而Logistic回归模型不需要；

4）Logistic回归没有关于自变量分布的假设条件，可以是连续变量、离散变量和虚拟变量；

5）由于因变量和自变量之间不存在线性关系，所以参数(偏回归系数)使用最大似然估计法计算。

logistic回归模型概述

        广义线性回归是探索“响应变量的期望”与“自变量”的关系，以实现对非线性关系的某种拟合。这里面涉及到一个“连接函数”和一个“误差函数”，“响应变量的期望”经过连接函数作用后，与“自变量”存在线性关系。选取不同的“连接函数”与“误差函数”可以构造不同的广义回归模型。当误差函数取“二项分布”而连接函数取“logit函数”时，就是常见的“logistic回归模型”，在0-1响应的问题中得到了大量的应用。​

         Logistic回归主要通过构造一个重要的指标：发生比来判定因变量的类别。在这里我们引入概率的概念，把事件发生定义为Y=1，事件未发生定义为Y=0，那么事件发生的概率为p,事件未发生的概率为1-p，把p看成x的线性函数；

         回归中，最常用的估计是最小二乘估计，因为使得p在[0,1]之间变换，最小二乘估计不太合适，有木有一种估计法能让p在趋近与0和1的时候变换缓慢一些（不敏感），这种变换是我们想要的，于是引入Logit变换,对p/(1-p)也就是发生与不发生的比值取对数，也称对数差异比。经过变换后，p对x就不是线性关系了。

logistic回归的公式可以表示为：

 

​其中P是响应变量取1的概率，在0-1变量的情形中，这个概率就等于响应变量的期望。

​这个公式也可以写成：

 

​可以看出，logistic回归是对0-1响应变量的期望做logit变换，然后与自变量做线性回归。参数估计采用极大似然估计，显著性检验采用似然比检验。

​建立模型并根据AIC准则选择模型后，可以对未知数据集进行预测，从而实现分类。模型预测的结果是得到每一个样本的响应变量取1的概率，为了得到分类结果，需要设定一个阈值p0——当p大于p0时，认为该样本的响应变量为1，否则为0。阈值大小对模型的预测效果有较大影响，需要进一步考虑。首先必须明确模型预测效果的评价指标。

​对于0-1变量的二分类问题，分类的最终结果可以用表格表示为：

​其中，d是“实际为1而预测为1”的样本个数，c是“实际为1而预测为0”的样本个数，其余依此类推。

显然地，主对角线所占的比重越大，则预测效果越佳，这也是一个基本的评价指标——总体准确率(a+d)/(a+b+c+d)。​

准确（分类）率=正确预测的正反例数/总数      Accuracy=(a+d)/(a+b+c+d)​​

误分类率=错误预测的正反例数/总数        Error rate=(b+c)/(a+b+c+d)=1-Accuracy

​正例的覆盖率=正确预测到的正例数/实际正例总数​

Recall(True Positive Rate，or Sensitivity)=d/(c+d)

正例的命中率=正确预测到的正例数/预测正例总数​

Precision(Positive Predicted Value,PV+)=d/(b+d)

负例的命中率=正确预测到的负例个数/预测负例总数

Negative predicted value(PV-)=a/(a+c)

通常将上述矩阵称为“分类矩阵”。一般情况下，我们比较关注响应变量取1的情形，将其称为Positive（正例），而将响应变量取0的情形称为Negative（负例）。常见的例子包括生物实验的响应、营销推广的响应以及信用评分中的违约等等。针对不同的问题与目的，我们通常采用ROC曲线与lift曲线作为评价logistic回归模型的指标。

​1）ROC曲线

设置了两个相应的指标：TPR与FPR。

TPR：True Positive Rate（正例覆盖率），将实际的1正确地预测为1的概率，d/(c+d)。

FPR：False Positive Rate，将实际的0错误地预测为1的概率，b/(a+b)。

TPR也称为Sensitivity（即生物统计学中的敏感度），也可以称为“正例的覆盖率”——将实际为1的样本数找出来的概率。覆盖率是重要的指标，例如若分类的目标是找出潜在的劣质客户（响应变量取值为1），则覆盖率越大表示越多的劣质客户被找出。

类似地，1-FPR其实就是“负例的覆盖率”，也就是把负例正确地识别为负例的概率。

​TPR与FPR相互影响，而我们希望能够使TPR尽量地大，而FPR尽量地小。影响TPR与FPR的重要因素就是上文提到的“阈值”。当阈值为0时，所有的样本都被预测为正例，因此TPR=1，而FPR=1。此时的FPR过大，无法实现分类的效果。随着阈值逐渐增大，被预测为正例的样本数逐渐减少，TPR和FPR各自减小，当阈值增大至1时，没有样本被预测为正例，此时TPR=0，FPR=0。

由上述变化过程可以看出，TPR与FPR存在同方向变化的关系（这种关系一般是非线性的），即，为了提升TPR（通过降低阈值），意味着FPR也将得到提升，两者之间存在类似相互制约的关系。我们希望能够在牺牲较少FPR的基础上尽可能地提高TPR，由此画出了ROC曲线。

ROC曲线的全称为“接受者操作特性曲线”（receiver operating characteristic），其基本形式为：

ROC曲线  

​当预测效果较好时，ROC曲线凸向左上角的顶点。平移图中对角线，与ROC曲线相切，可以得到TPR较大而FPR较小的点。模型效果越好，则ROC曲线越远离对角线，极端的情形是ROC曲线经过（0，1）点，即将正例全部预测为正例而将负例全部预测为负例。ROC曲线下的面积可以定量地评价模型的效果，记作AUC，AUC越大则模型效果越好。

当我们分类的目标是将正例识别出来时（例如识别有违约倾向的信用卡客户），我们关注TPR，此时ROC曲线是评价模型效果的准绳。

​2）lift曲线

在营销推广活动中，我们的首要目标并不是尽可能多地找出那些潜在客户，而是提高客户的响应率。客户响应率是影响投入产出比的重要因素。此时，我们关注的不再是TPR（覆盖率），而是另一个指标：命中率。

回顾前面介绍的分类矩阵，正例的命中率是指预测为正例的样本中的真实正例的比例，即d/(b+d)，一般记作PV。

在不使用模型的情况下，我们用先验概率估计正例的比例，即(c+d)/(a+b+c+d)，可以记为k。

定义提升值lift=PV/k。

lift揭示了logistic模型的效果。例如，若经验告诉我们10000个消费者中有1000个是我们的潜在客户，则我们向这10000个消费者发放传单的效率是10%（即客户的响应率是10%），k=(c+d)/(a+b+c+d)=10%。通过对这10000个消费者进行研究，建立logistic回归模型进行分类，我们得到有可能比较积极的1000个消费者，b+d=1000。如果此时这1000个消费者中有300个是我们的潜在客户，d=300，则命中率PV为30%。此时，我们的提升值lift=30%/10%=3，客户的响应率提升至原先的三倍，提高了投入产出比。

为了画lift图，需要定义一个新的概念depth深度，这是预测为正例的比例，(b+d)/(a+b+c+d)。

与ROC曲线中的TPR和FPR相同，lift和depth也都受到阈值的影响。

当阈值为0时，所有的样本都被预测为正例，因此depth=1，而PV=d/(b+d)=(0+d)/(0+b+0+d)=k，于是lift=1，模型未起提升作用。随着阈值逐渐增大，被预测为正例的样本数逐渐减少，depth减小，而较少的预测正例样本中的真实正例比例逐渐增大。当阈值增大至1时，没有样本被预测为正例，此时depth=0，而lift=0/0。

由此可见，lift与depth存在相反方向变化的关系。在此基础上作出lift图：

lift 曲线

​ 与ROC曲线不同，lift曲线凸向（0，1）点。我们希望在尽量大的depth下得到尽量大的lift（当然要大于1），也就是说这条曲线的右半部分应该尽量陡峭。

至此，我们对ROC曲线和lift曲线进行了描述。这两个指标都能够评价logistic回归模型的效果，只是分别适用于不同的问题：

如果是类似信用评分的问题，希望能够尽可能完全地识别出那些有违约风险的客户（不使一人漏网），我们需要考虑尽量增大TPR（覆盖率），同时减小FPR（减少误杀），因此选择ROC曲线及相应的AUC作为指标；

如果是做类似数据库精确营销的项目，希望能够通过对全体消费者的分类而得到具有较高响应率的客户群，从而提高投入产出比，我们需要考虑尽量提高lift（提升度），同时depth不能太小（如果只给一个消费者发放传单，虽然响应率较大，却无法得到足够多的响应），因此选择lift曲线作为指标。

相关R应用包

普通二分类 logistic 回归 用系统的 glm

因变量多分类 logistic 回归

有序分类因变量：用 MASS 包里的 polrb

无序分类因变量：用 nnet 包里的 multinom

条件logistic回归，用 survival 包里的 clogit

R运行逻辑回归

R可以让逻辑回归建模过程变得很简单。我们可以使用glm()函数进行逻辑回归建模，而且，它的拟合过程和我们之前所做的线性回归没有很大的区别。在这篇博客中，我们将介绍如何一步一步的进行逻辑回归建模。

我们用r来看一下曲线的形状：

x <- seq(from = -10, to = 10, by = 0.01)
y = exp(x)/(1+exp(x))
library(ggplot2)
p <- ggplot(data = NULL, mapping = aes(x = x,y = y))
p + geom_line(colour = 'blue')
+ annotate('text', x = 1, y = 0.3, label ='y==e^x / 1-e^x', parse = TRUE)
+ ggtitle('Logistic曲线')

p在0和1附近变化不敏感，ps:这条曲线帮了我一个大忙，最近做RTB竞价算法时，需要按照CTR设定一个CTR指导价，然后以CPC做曲线的修正，出于成本考虑，我们的出价不会高于某个固定值K，所以把这条曲线变换一下就派上大用场了。（扯远了）

这条曲线的主要特性是θ函数取值可以在-∞变到+∞，p的取值在[0,1],且极值端变化不敏感。这就是我们想要的。值得注意的是，经过logit变换，已经不是线性模型。

二、相关应用例子​：Binary Logistic（因变量只能取两个值1和0虚拟因变量）​

案例一：本文用例来自于John Maindonald所著的《Data Analysis and Graphics Using R》一书，其中所用的数据集是anesthetic，数据集来自于一组医学数据，其中变量conc表示麻醉剂的用量，move则表示手术病人是否有所移动，而我们用nomove做为因变量，因为研究的重点在于conc的增加是否会使nomove的概率增加。

​首先载入数据集并读取部分文件，为了观察两个变量之间关系，我们可以利cdplot函数来绘制条件密度图

install.packages("DAAG")​

library(lattice)​

library(DAAG)​

head(anesthetic)​

    move conc   logconc nomove

1    0  1.0 0.0000000      1

2    1  1.2 0.1823216      0

3    0  1.4 0.3364722      1

4    1  1.4 0.3364722      0

5    1  1.2 0.1823216      0

6    0  2.5 0.9162907      1

cdplot(factor(nomove)~conc,data=anesthetic,main='条件密度图',ylab='病人移动',xlab='麻醉剂量')​

从图中可见，随着麻醉剂量加大，手术病人倾向于静止。下面利用logistic回归进行建模，得到intercept和conc的系数为-6.47和5.57，由此可见麻醉剂量超过1.16(6.47/5.57)时，病人静止概率超过50%。​

anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic)​

summary(anes1)

结果显示：

Call:

glm(formula = nomove ~ conc, family = binomial(link = "logit"), 

    data = anesthetic)

Deviance Residuals: 

     Min        1Q    Median        3Q       Max  

-1.76666  -0.74407   0.03413   0.68666   2.06900  

Coefficients:

                   Estimate Std. Error    z value        Pr(>|z|)   

(Intercept)   -6.469      2.418       -2.675         0.00748 **

conc             5.567      2.044        2.724          0.00645 **

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 41.455  on 29  degrees of freedom

Residual deviance: 27.754  on 28  degrees of freedom

AIC: 31.754

Number of Fisher Scoring iterations: 5​

 下面做出模型的ROC曲线

anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic)​

对模型做出预测结果​

pre=predict(anes1,type='response')​

将预测概率pre和实际结果放在一个数据框中​

data=data.frame(prob=pre,obs=anesthetic$nomove)​

将预测概率按照从低到高排序​​

data=d