Logistic回归详解(一)

发表: 2017-03-29 浏览: 3350

前言

   Logistic回归应该是大家再熟悉不过的针对二分类或多分类因变量的回归模型之一。由于其概率表达式的显性特点，模型的可解释性强，是社会学、生物统计学、临床、数量心理学、市场营销等统计实证分析的常用方法。但是在模型的运用方面，由于不理解模型的应用特点，出现非常多的“生搬硬套”嫌疑。

    为了更深入学习Logistic回归，我将从常规Logistic回归、确切Logistic回归（Exact Logistic Regression）、多分类Logistic回归、有序Logistic回归及Probit 回归这些回归模型进行简要的介绍及R实现。由于涉及内容较多，自身也在学习当中，就学到哪写到哪，分几期来完成它。这一期让我们来看看常规的Logistic回归及R实现。

Logistic回归概述

    Logistic回归最常用的是解决二分类问题，即预测或判别y=0或1。这就期望我们能构建一个分类器h(x)在0-1之间，于是就有了Logistic函数或Logistic曲线，又叫 sigmoid曲线（S型曲线）。我们尽可能的少用公式展示，一方面是不好编辑，另一方面不希望大家有望而生畏的错觉。

  由上图可以看到，当z<0时，h(x)<0.5;z>= 0时，h(x)>=0.5,且h(x)在0-1 之间 ，很好的构建了一个符合我们要求的二分类器。为更好的理解Logistic回归，我们再引入一个例子：

    研究者希望通过查看学生的GRE成绩，平均绩点GPA及毕业学校的威望排名来评估是否录取该学生，即因变量y是录取（y=1）或不录取（y=0），二分类问题。

#读入数据mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")head(mydata)

##   admit gre  gpa rank

## 1     0 380 3.61    3

## 2     1 660 3.67    3

## 3     1 800 4.00    1

## 4     1 640 3.19    4

## 5     0 520 2.93    4

## 6     1 760 3.00    2

    我们的目的是通过已知某个学生的GRE成绩、GPA及毕业学校来预测该生被录取的概率，也就是h(x)为纳入X自变量后，y=1时的概率值。当这个概率值h(x)<0.5时，认为不会被录取（y=0），反之则被录取。总之一点，我们求得的h(x)是y=1的概率值。

Logistic回归模型的适用条件

这部分是百度上的，适用条件其实和线性回归模型差不多，应该不难理解。

因变量y为二分类的分类变量或某事件的发生率
残差和因变量都要服从二项分布,运用最大似然法来解决方程估计和检验问题。（至于什么是二项分布，什么是最大似然法，看不懂的可以略过）
自变量和Logistic概率是线性关系
各观测对象间相互独立

模型实现

运用广义线性模型glm实现Logistic回归.

library(aod)

library(ggplot2)
library(Rcpp)
mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")summary(mylogit)

## 

## Call:

## glm(formula = admit ~ gre + gpa + rank, family = "binomial", 

##     data = mydata)

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.6268  -0.8662  -0.6388   1.1490   2.0790  

## 

## Coefficients:

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

## (Intercept) -3.989979   1.139951  -3.500 0.000465 ***

## gre          0.002264   0.001094   2.070 0.038465 *  

## gpa          0.804038   0.331819   2.423 0.015388 *  

## rank2       -0.675443   0.316490  -2.134 0.032829 *  

## rank3       -1.340204   0.345306  -3.881 0.000104 ***

## rank4       -1.551464   0.417832  -3.713 0.000205 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for binomial family taken to be 1)

## 

##     Null deviance: 499.98  on 399  degrees of freedom

## Residual deviance: 458.52  on 394  degrees of freedom

## AIC: 470.52

## 

## Number of Fisher Scoring iterations: 4

#系数95%CI
confint(mylogit)

##                     2.5 %       97.5 %

## (Intercept) -6.2716202334 -1.792547080

## gre          0.0001375921  0.004435874

## gpa          0.1602959439  1.464142727

## rank2       -1.3008888002 -0.056745722

## rank3       -2.0276713127 -0.670372346

## rank4       -2.4000265384 -0.753542605

## odds ratios（OR值）及95%CI
exp(cbind(OR = coef(mylogit), confint(mylogit)))

##                    OR       2.5 %    97.5 %

## (Intercept) 0.0185001 0.001889165 0.1665354

## gre         1.0022670 1.000137602 1.0044457

## gpa         2.2345448 1.173858216 4.3238349

## rank2       0.5089310 0.272289674 0.9448343

## rank3       0.2617923 0.131641717 0.5115181

## rank4       0.2119375 0.090715546 0.4706961

    模型的结果给出了模型的结构，Deviance Residuals、beta系数 及P值和模型 评价指标：Null deviance、Residual deviance、AIC。

    模型解释：对于GRE每多考一分，被录取的概率的log值（相对于未被录取）增加0.002；GPA每多考一分，被录取的概率的log值增加0.804；毕业学威望rank是2的（与威望rank是1比较）被录取的概率的log值降低0.675。由于有log的变换，解释会比较奇怪，因此取e变换后即为我们常见的OR值。

内部变量的比较

    运用wald.test（aod包）对rank变量的总效应进行测试。b为模型的系数，Sigma为变量协方差矩阵，Terms为需要比较的内部变量，这里指Rank变量的3个水平，位置4-6。卡方为20.9，p<0.001,说明rank变量有统计学意义。

wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)

## Wald test:

## ----------

## Chi-squared test:

## X2 = 20.9, df = 3, P(> X2) = 0.00011

    假如我们想比较rank=2与rank=3的系数值是否相等，可以用如下操作：

卡方为5.5，p<0.05,有统计学意义，说明rank=2与rank=3的系数不相等，差异有意义。

l <- cbind(0,0,0,1,-1,0)
wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)

## Wald test:

## ----------

## Chi-squared test:

## X2 = 5.5, df = 1, P(> X2) = 0.019

模型预测

我们构建一个新的数据框，来预测下新数据的概率值：

newdata1 <- with(mydata,data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))
## view data frame
newdata1

##     gre    gpa rank

## 1 587.7 3.3899    1

## 2 587.7 3.3899    2

## 3 587.7 3.3899    3

## 4 587.7 3.3899    4

newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response")
newdata1

##     gre    gpa rank     rankP

## 1 587.7 3.3899    1 0.5166016

## 2 587.7 3.3899    2 0.3522846

## 3 587.7 3.3899    3 0.2186120

## 4 587.7 3.3899    4 0.1846684

    在以上的结果中，取平均GRE和GPA,我们看到威望rank=1的毕业学校的学生被录取的概率为0.517，而rank=4的学生被录取的概率只有0.185（说明教育背景真的很重要）。

    我们还可以创建一个图表，查看GRE得分与rank在录取概率上的关系。

newdata2 <- with(mydata,  data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100), 4),  gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))
newdata3 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type="link", se=TRUE))
newdata3 <- within(newdata3, {

  PredictedProb <- plogis(fit)

  LL <- plogis(fit - (1.96 * se.fit))

  UL <- plogis(fit + (1.96 * se.fit))})
head(newdata3)

##        gre    gpa rank        fit    se.fit residual.scale        UL

## 1 200.0000 3.3899    1 -0.8114870 0.5147714              1 0.5492064

## 2 206.0606 3.3899    1 -0.7977632 0.5090986              1 0.5498513

## 3 212.1212 3.3899    1 -0.7840394 0.5034491              1 0.5505074

## 4 218.1818 3.3899    1 -0.7703156 0.4978239              1 0.5511750

## 5 224.2424 3.3899    1 -0.7565919 0.4922237              1 0.5518545

## 6 230.3030 3.3899    1 -0.7428681 0.4866494              1 0.5525464

##          LL PredictedProb

## 1 0.1393812     0.3075737

## 2 0.1423880     0.3105042

## 3 0.1454429     0.3134499

## 4 0.1485460     0.3164108

## 5 0.1516973     0.3193867

## 6 0.1548966     0.3223773

ggplot(newdata3, aes(x = gre, y = PredictedProb)) +

  geom_ribbon(aes(ymin = LL, ymax = UL, fill = rank), alpha = .2) +

  geom_line(aes(colour = rank), size=1)

模型评价

    最后再来看看模型的拟合好坏，通过summary(mylogit)已经看到了模型评价的指标，现在来做下模型的似然估计测试。结果p<0.001，说明模型拟合有意义。

with(mylogit, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))

## [1] 7.578194e-08

####模型ROC曲线
pre <- predict(mylogit,data = mydata[,2:4])
library(pROC)

modelroc <- roc(mydata$admit,pre,plot = TRUE,print.thres=TRUE, print.auc=TRUE)

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="skyblue", print.thres=TRUE)

## 

## Call:

## roc.default(response = mydata$admit, predictor = pre, plot = TRUE,     print.thres = TRUE, print.auc = TRUE)

## 

## Data: pre in 273 controls (mydata$admit 0) < 127 cases (mydata$admit 1).

## Area under the curve: 0.6928

可以看出，模型的AUC为0.693，在截断点为-0.603时，取值模型敏感性0.744，特异性0.575.

#模型评价
pres <- ifelse(pre >= -0.603,1,0)pres <- as.factor(pres)
#模型混淆矩阵
xtabs(~pres + mydata$admit)

##     mydata$admit

## pres   0   1

##    0 203  54

##    1  70  73

参考文献

idre institute for digital research and education R Data Analysis Examples: Logit Regressionhttp://www.ats.ucla.edu/stat/r/dae/logit.htm

coursera 机器学习 Logistic regression

logistic回归http://baike.baidu.com/link?url=VjGTrRmeLGHxVxDSrNQg_ZuzV_tuutEjB_3M2EiULyC22lMDejqhRot2Pj5Ji0IrX0EB8nZJo1m5apqMKGaJea

维基百科 Logistic regression https://en.wikipedia.org/wiki/Logistic_regression

数据科学中的R和Python: 浅谈ROC曲线http://www.360doc.com/content/15/0326/21/17553313_458307034.shtml

0 个评论

要回复文章请先登录或注册