Logistic回归模型、应用建模案例

一、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曲线作为指标

3)相关R应用包

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

因变量多分类 logistic 回归

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

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

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

二、相关应用例子: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=data[order(data$prob),]

n=nrow(data)

tpr=fpr=rep(0,n)

根据不同的临界值threshold来计算TPR和FPR,之后绘制成图

for (i in 1:n){

threshold=data$prob[i]

tp=sum(data$prob>threshold&data$obs==1)

fp=sum(data$prob>threshold&data$obs==0)

tn=sum(data$prob

fn=sum(data$prob

tpr[i]=tp/(tp+fn) #真正率

fpr[i]=fp/(tn+fp) #假正率

}

plot(fpr,tpr,type='l')

abline(a=0,b=1)

R中也有专门绘制ROC曲线的包,如常见的ROCR包,它不仅可以用来画图,还能计算ROC曲线下面面积AUC,以评价分类器的综合性能,该数值取0-1之间,越大越好。

library(ROCR)

pred=prediction(pre,anesthetic$nomove)

performance(pred,'auc')@y.values

perf=performance(pred,'tpr','fpr')

plot(perf)

还可以使用更加强大的pROC包,它可以方便的比较两个分类器,并且能自动标出最优临界点,图形看起来比较漂亮:

install.packages("pROC")

library(pROC)

modelroc=roc(anesthetic$nomove,pre)

plot(modelroc,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=TRUE,auc.polygon.col="blue",print.thres=TRUE)

上面的方法是使用原始的0-1数据进行建模,即每一行数据均表示一个个体,另一种是使用汇总数据进行建模,先将原始数据按下面步骤进行汇总

anestot=aggregate(anesthetic[,c('move','nomove')],by=list(conc=anesthetic$conc),FUN=sum)

结果如下:

conc move nomove

1 0.8 6 1

2 1.0 4 1

3 1.2 2 4

4 1.4 2 4

5 1.6 0 4

6 2.5 0 2

anestot$conc=as.numeric(as.character(anestot$conc))

anestot$total=apply(anestot[,c('move','nomove')],1,sum)

anestot$total

[1] 7 5 6 6 4 2

anestot$prop=anestot$nomove/anestot$total

anestot$prop

[1] 0.1428571 0.2000000 0.6666667 0.6666667 1.0000000 1.0000000

对于汇总数据,有两种方法可以得到同样的结果,一种是将两种结果的向量合并做为因变量,如anes2模型。另一种是将比率做为因变量,总量做为权重进行建模,如anes3模型。这两种建模结果是一样的。

anes2=glm(cbind(nomove,move)~conc,family=binomial(link='logit'),data=anestot)

summary(anes2)

结果显示如下:

Call:

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

data = anestot)

Deviance Residuals:

1 2 3 4 5 6

0.20147 -0.45367 0.56890 -0.70000 0.81838 0.04826

Coefficients:

Estimate Std. Error z value Pr(>z)

(Intercept) -6.469 2.419 -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: 15.4334 on 5 degrees of freedom

Residual deviance: 1.7321 on 4 degrees of freedom

AIC: 13.811

Number of Fisher Scoring iterations: 5

anes3=glm(prop~conc,family=binomial(link='logit'),weights=total,data=anestot)

结果和上面的一样。

根据logistic模型,我们可以使用predict函数来预测结果,下面根据上述模型来绘图

x=seq(from=0,to=3,length.out=30)

y=predict(anes1,data.frame(conc=x),type='response')

plot(prop~conc,pch=16,col='red',data=anestot,xlim=c(0.5,3),main='Logistic回归曲线图',ylab='病人静止概率',xlab='麻醉剂量')

lines(y~x,lty=2,col='blue')

案例二:利用iris数据集,进行逻辑回归二分类测试,该数据集是R语言自带得数据集,包括四个属性,和三个分类。逻辑回归我们用glm函数实现,该函数提供了各种类型的回归,如:提供正态、指数、gamma、逆高斯、Poisson、二项。我们用的logistic回归使用的是二项分布族binomial。

index <- which(iris$Species == 'setosa')

将种类为setosa的数据排除出我们需要的数据集

ir <- iris[- index,]

levels(ir$Species)[1] <- ''

生成训练集

split <- sample(100,100*(2/3))

ir_train <- ir[split,]

生成测试集 ir_test <- ir[-split,]

通过训练集建立模型 model <- glm(Species ~.,family=binomial(link='logit'),data=ir_train)

summary(model)

模型运行结果:

Call:

glm(formula = Species ~ ., family = binomial(link = "logit"),data = ir_train)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.339e-04 -2.100e-08 2.100e-08 2.100e-08 1.059e-04

Coefficients:

Estimate Std. Error z value Pr(>z)

(Intercept) -1502.72 363247.01 -0.004 0.997

Sepal.Length 12.45 66482.13 0.000 1.000

Sepal.Width -285.61 95437.92 -0.003 0.998

Petal.Length 154.76 115968.97 0.001 0.999

Petal.Width 869.60 204513.80 0.004 0.997

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 9.0949e+01 on 65 degrees of freedom

Residual deviance: 4.0575e-08 on 61 degrees of freedom

AIC: 10

Number of Fisher Scoring iterations: 25

通过anova()函数 对模型进行方差分析

anova(model, test="Chisq")

方差分析如下:

Analysis of Deviance Table

Model: binomial, link: logit

Response: Species

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev Pr(>Chi)

NULL 65 90.949

Sepal.Length 1 18.934 64 72.015 1.353e-05 ***

Sepal.Width 1 0.131 63 71.884 0.7176

Petal.Length 1 51.960 62 19.924 5.665e-13 ***

Petal.Width 1 19.924 61 0.000 8.058e-06 ***

---

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

下面通过McFadden R2指标进一步对模型进行分析

install.packages("pscl")

library(pscl)

pR2(model)

llh llhNull G2 McFadden r2ML r2CU

-2.028752e-08 -4.547461e+01 9.094922e+01 1.000000e+00 7.479224e-01 1.000000e+00

为了得到分类结果,需要设定一个阈值p0——当p大于p0时,认为该样本的响应变量为1,否则为0。阈值大小对模型的预测效果有较大影响,需要进一步考虑。首先必须明确模型预测效果的评价指标。

求解训练模型的最佳阀值

对模型做出预测结果

model <- glm(Species ~.,family=binomial(link='logit'),data=ir_train)

pre1=predict(model,type='response')

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

data1=data.frame(prob=pre1,obs=ifelse(ir_train$Species=="virginica",1,0))

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

data1=data1[order(data1$prob),]

n=nrow(data1)

tpr=fpr=rep(0,n)

根据不同的临界值threshold来计算TPR和FPR,之后绘制成图

for (i in 1:n){

threshold=data1$prob[i]

tp=sum(data1$prob>threshold&data1$obs==1)

fp=sum(data1$prob>threshold&data1$obs==0)

tn=sum(data$prob

fn=sum(data$prob

tpr[i]=tp/(tp+fn) #真正率

fpr[i]=fp/(tn+fp) #假正率

}

plot(fpr,tpr,type='l')

abline(a=0,b=1)

下面通过pROC包自动标出最优临界点(0.506)

install.packages("pROC")

library(pROC)

modelroc1=roc(ifelse(ir_train$Species=="virginica",1,0),pre1)

plot(modelroc1,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=TRUE,auc.polygon.col="skyblue",print.thres=TRUE)

评估模型的预测效果

predict <- predict(model,type='response',newdata=ir_test)

predict.results <- ifelse( predict> 0.506,"virginica","versicolor")

misClasificError <- mean(predict.results != ir_test$Species)

print(paste('Accuracy',1-misClasificError))

[1] "Accuracy 1"

最后一步,我们将通过画ROC曲线,并计算其AUC面积,作为评估二类分类效果的一个典型测量

install.packages("ROCR")

library(gplots)

library(ROCR)

p <- predict(model,type='response',newdata=ir_test)

p.results <- ifelse( p> 0.5,1,0)

pr <- prediction(p.results, ifelse(ir_test$Species=="virginica",1,0))

prf <- performance(pr, measure = "tpr", x.measure = "fpr")

plot(prf)

auc <- performance(pr, measure = "auc")

auc <- auc@y.values[[1]]

0.9285714

auc

real <- ir_test$Species

data.frame(real,predict)

res <- data.frame(real,predict =ifelse(predict>0.5,'virginca','versicorlor'))

查看模型效果

plot(res)


原文发布于微信公众号 - 大数据挖掘DT数据分析(datadw)

原文发表时间:2016-10-04

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏AI2ML人工智能to机器学习

一个奇异值的江湖 -- 经典统计观

有江湖的地方就有奇异值, 有时候奇值(outlier)和异值(anomaly)会咬文嚼字的区分一下: outlier是合理的(explainable)小概率事件...

371
来自专栏编程

Python环境下的8种简单线性回归算法

选自Medium 作者:Tirthajyoti Sarkar 机器之心编译 参与:晏奇、刘晓坤 本文中,作者讨论了 8 种在 Python 环境下进行简单线性回...

2009
来自专栏机器之心

入门 | 从结构到性能,一文概述XGBoost、Light GBM和CatBoost的同与不同

选自Medium 机器之心编译 参与:刘天赐、黄小天 尽管近年来神经网络复兴并大为流行,但是 boosting 算法在训练样本量有限、所需训练时间较短、缺乏调参...

3575
来自专栏人工智能头条

LSTM实现详解

1753
来自专栏大数据风控

R中的线性回归分析

回归分析(regression analysis) 回归分析是研究自变量与因变量之间关系形式的分析方法,它主要是通过建立因变量Y与影响它的自变量Xi(i=1,2...

18710
来自专栏智能算法

各大公司广泛使用的在线学习算法FTRL详解

现在做在线学习和CTR常常会用到逻辑回归( Logistic Regression),而传统的批量(batch)算法无法有效地处理超大规模的数据集和在线数据流,...

4606
来自专栏悦思悦读

三个数值告诉你分类器的优劣

预测和测试 通过训练得到机器学习模型后,我们需要用模型来对用户不断输入的语句进行预测(也就是把用户语句输入到模型中让模型吐出一个结果)。 预测肯定能出结果,至于...

3266
来自专栏小樱的经验随笔

MATLAB学习笔记

魔方矩阵(magic(阶数)) 魔方矩阵又称幻方,是有相同的行数和列数,并在每行每列、对角线上的和都相等的矩阵。魔方矩阵中的每个元素不能相同。你能构造任何大小(...

3414
来自专栏我的python

char-rnn 语言模型

项目GitHub地址:https://github.com/ClownW/Char-RNN-Pytorch

2136
来自专栏机器之心

教程 | 如何使用TensorFlow构建、训练和改进循环神经网络

选自SVDS 作者:Matthew Rubashkin、Matt Mollison 机器之心编译 参与:李泽南、吴攀 来自 Silicon Valley Dat...

2849

扫码关注云+社区