2020-03-06 如何在R中进行Logistic回归

今天被问及到相关的问题,回答不上来,所以还是找来学习资料学习一下,数学模型什么的太抽象了。

回归正题:
Logistic回归是一种拟合回归曲线的方法(y = f(x),),当y是分类变量时。该模型的典型应用是在给定一组预测值x的情况下预测y。预测值可以是连续的变量,也可以是离散的分类变量,或者两者兼而有之。
分类变量y通常可以采用不同的值。在最简单的情况下,y是个逻辑值,这意味着它可以取值1或0。机器学习中使用的一个经典例子是电子邮件分类:给定每个电子邮件的一组属性,如单词、链接和图片的数量,算法应该确定该电子邮件是否为垃圾邮件(1)或不是(0)。在这篇文章中,我们称该模型为“二项Logistic回归”,因为要预测的变量是二元的,但是,Logistic回归也可以用来预测一个可以假设两个以上值的因变量。在第二种情况下,我们称该模型为“多项Logistic回归”。例如,一个典型的例子是将电影分类为“娱乐性的”、“中性的”或“无聊的”( “Entertaining”, “borderline” or “boring”)。

使用R进行拟合Logistic回归模型非常容易。要调用的函数是glm(),拟合过程与线性回归中使用的过程没有太大不同。接下来,我将建立一个二元Logistic回归模型,并解释每个步骤。

获取示例数据集

我们将研究泰坦尼克号的数据集。这个数据集在网上有不同的免费版本,但是我建议使用Kaggle上提供的版本,这个比较便利(要下载它,您需要注册Kaggle)。Kaggle网站上有很多数据集,包括现在特别被关注的2019新冠肺炎相关的数据,都可以注册后获取。


image.png

image.png

首先进行数据集(training),这个数据集是关于一些乘客(准确地说是889人)的数据集合,training目标是基于一些特征(如服务等级、性别、年龄等)预测其生存状况(如果乘客幸存,则为1;如果死亡,则为0)。正如您所看到的,我们将使用分类变量和连续变量。

数据清洗

当使用真实的数据集时,我们需要考虑到某些数据可能丢失或损坏的事实,因此我们需要为我们要分析的数据进行清洗。第一步,使用read.csv()函数读取数据。
确保参数na.string等于c(""),以便将每个缺少的值编码为NA。这将有助于下一步的工作。

> rm(list = ls())
> options(stringsAsFactors = F)
> training.data.raw <- read.csv('train.csv',header=T,na.strings=c(""))
> dim(training.data.raw)
[1] 891  12

现在,我们需要检查缺少的值,并使用sapply()函数查看每个变量有多少个唯一值,该函数将作为参数传递的函数应用到数据集的每一列。

> sapply(training.data.raw,function(x) sum(is.na(x)))
PassengerId    Survived      Pclass        Name         Sex         Age       SibSp       Parch 
          0           0           0           0           0         177           0           0 
     Ticket        Fare       Cabin    Embarked 
          0           0         687           2 
> sapply(training.data.raw, function(x) length(unique(x)))
PassengerId    Survived      Pclass        Name         Sex         Age       SibSp       Parch 
        891           2           3         891           2          89           7           7 
     Ticket        Fare       Cabin    Embarked 
        681         248         148           4 

直观地了解缺少的值可能会有所帮助:Amelia包有一个特殊的打印函数missmap(),该函数将绘制您的数据集并突出显示缺少的值:

> library(Amelia)
Loading required package: Rcpp
## 
## Amelia II: Multiple Imputation
## (Version 1.7.6, built: 2019-11-24)
## Copyright (C) 2005-2020 James Honaker, Gary King and Matthew Blackwell
## Refer to http://gking.harvard.edu/amelia/ for more information
## 
> missmap(training.data.raw, main = "Missing values vs observed")

出图如下:


空白处则为缺失值

变量cabin缺失值太多,我们将不使用它。我们还将删除PassengerId,因为它只是一个索引和票证。
使用subset()函数,我们将原始数据集通过仅选择相关列取子集。

> data <- subset(training.data.raw,select=c(2,3,5,6,7,8,10,12))

处理缺失值

现在我们需要说明其他缺失值。在拟合广义线性模型时,R通过在拟合函数中设置一个参数,可以很容易地处理这些问题。不过,就我个人而言,在可能的情况下,我更喜欢“手动”更换NAs。有不同的方法可以做到这一点,典型的方法是用现有的平均值、中位数或模式替换缺失的值。我会用平均数来替换。

> data$Age[is.na(data$Age)] <- mean(data$Age,na.rm=T)

就分类变量而言,默认情况下使read.table()read.csv()会将分类变量编码为因子。因子是R如何处理分类变量的依据。
我们可以使用以下代码行检查编码

> is.factor(data$Sex)
[1] FALSE
> is.factor(data$Embarked)
[1] FALSE

这个结果是由于一开始我启用了options(stringsAsFactors = F)。回去重新修改一下参数为T。
为了更好地理解R将如何处理类别变量,我们可以使用conflsts()函数。这个函数将向我们展示变量是如何被R虚化的,以及如何在模型中解释它们。

> contrasts(as.factor(data$Sex))
       male
female    0
male      1
> contrasts(as.factor(data$Embarked))
  Q S
C 0 0
Q 1 0
S 0 1

比如你可以看到,在这个分类变量性别中,会用女性作为参照物。至于Emploked中缺失的值,因为只有两个,我们将丢弃这两行(我们也可以用模式替换缺失的值,并保留数据点)。

> data <- data[!is.na(data$Embarked),]
> rownames(data) <- NULL

在进行拟合过程之前,数据的清理和格式化非常重要。这一预处理步骤对于获得良好的模型拟合和更好的预测能力通常是至关重要的。

建模

我们将数据分成两个块:训练集和测试集。训练集将用于拟合我们的模型,我们将在测试集上测试该模型。

> train <- data[1:800,]
> test <- data[801:889,]

现在,让我们来试一下构建模型。确保在函数glm()中指定参数family=binomial。用函数summary()获取模型具体信息

> model <- glm(Survived ~.,family=binomial(link='logit'),data=train)
> summary(model)

Call:
glm(formula = Survived ~ ., family = binomial(link = "logit"), 
    data = train)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6064  -0.5954  -0.4254   0.6220   2.4165  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.137627   0.594998   8.635  < 2e-16 ***
Pclass      -1.087156   0.151168  -7.192 6.40e-13 ***
Sexmale     -2.756819   0.212026 -13.002  < 2e-16 ***
Age         -0.037267   0.008195  -4.547 5.43e-06 ***
SibSp       -0.292920   0.114642  -2.555   0.0106 *  
Parch       -0.116576   0.128127  -0.910   0.3629    
Fare         0.001528   0.002353   0.649   0.5160    
EmbarkedQ   -0.002656   0.400882  -0.007   0.9947    
EmbarkedS   -0.318786   0.252960  -1.260   0.2076    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1065.39  on 799  degrees of freedom
Residual deviance:  709.39  on 791  degrees of freedom
AIC: 727.39

Number of Fisher Scoring iterations: 5

解释我们Logistic回归模型的结果

现在我们可以分析拟合并解释模型告诉我们的内容。
首先,我们可以看到,SibSp, Fare 和 Embarked 没有统计学意义。至于有统计学意义的变量,性别的p值最低,表明乘客的性别与生存概率有很强的关联。这个预测值的负系数表明,在所有其他变量相同的情况下,男性乘客生还的可能性较小。请记住,在LOGIT模型中,响应变量是log oddsln(odds) = ln(p/(1-p)) = a*x1 + b*x2 + … + z*xn。由于男性是一个虚拟变量,作为男性会使log odds降低2.75%,而单位年龄增加则会使log odds降低0.037。
现在,我们可以在模型上运行函数anova()来分析偏差:

> anova(model, test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: Survived

Terms added sequentially (first to last)


         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                       799    1065.39              
Pclass    1   83.607       798     981.79 < 2.2e-16 ***
Sex       1  240.014       797     741.77 < 2.2e-16 ***
Age       1   17.495       796     724.28 2.881e-05 ***
SibSp     1   10.842       795     713.43  0.000992 ***
Parch     1    0.863       794     712.57  0.352873    
Fare      1    0.994       793     711.58  0.318717    
Embarked  2    2.187       791     709.39  0.334990    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

null偏差和residual偏差之间的差异显示了我们的模型相对于null模型(只有截距的模型)的表现。这个差距越大越好。分析表格,我们可以看到每次添加一个变量时偏差的下降。再次,增加PclassSexAge显著降低残差。其他变量似乎对模型的改善较小,即使SibSpp值较低。这里的p值很大,表示没有变量的模型解释的变异量大致相同。归根结底,你希望看到的是偏差和AIC的显着下降。
在不存在与线性回归的R^2完全等价的情况下,可以使用McFadden R^2指数来评估模型的拟合程度。

> library(pscl)
Classes and Methods for R developed in the
Political Science Computational Laboratory
Department of Political Science
Stanford University
Simon Jackman
hurdle and zeroinfl functions by Achim Zeileis
> pR2(model)
         llh      llhNull           G2     McFadden         r2ML         r2CU 
-354.6950111 -532.6961008  356.0021794    0.3341513    0.3591775    0.4880244 

评估模型的预测能力

在上面的步骤中,我们简要地评估了模型的拟合情况,现在我们想看看模型在预测一组新数据时的表现。通过设置参数type='response',R将以P(y=1|X)的形式输出概率。我们的决定边界是0.5。如果P(y=1|X) > 0.55,则y=1,否则y=0。请注意,对于某些应用程序,不同的阈值可能是更好的选择。

> fitted.results <- predict(model,newdata=subset(test,select=c(2,3,4,5,6,7,8)),type='response')
> fitted.results <- ifelse(fitted.results > 0.5,1,0)
> misClasificError <- mean(fitted.results != test$Survived)
> print(paste('Accuracy',1-misClasificError))
[1] "Accuracy 0.842696629213483"

在测试集上达到0.84的准确率是一个相当好的结果。但是,请记住,此结果在某种程度上依赖于我前面所做的手动数据拆分,因此,如果您希望获得更精确的分数,最好运行某种交叉验证,如k-折交叉验证。
作为最后一步,我们将绘制ROC曲线并计算AUC(曲线下面积),这是二元分类器的典型性能测量。
ROC是在各种阈值设置下绘制真阳性率(TPR)与假阳性率(FPR)的曲线,而AUC是ROC曲线下的面积。根据经验,具有良好预测能力的模型的AUC应该更接近1(1是理想的),而不是0.5。

> library(ROCR)
Loading required package: gplots

Attaching package: 'gplots'

The following object is masked from 'package:stats':

    lowess

> p <- predict(model, newdata=subset(test,select=c(2,3,4,5,6,7,8)), type="response")
> pr <- prediction(p, test$Survived)
> prf <- performance(pr, measure = "tpr", x.measure = "fpr")
> plot(prf)
> auc <- performance(pr, measure = "auc")
> auc <- auc@y.values[[1]]
> auc
[1] 0.8647186
ROC plot

这个教程还是非常的有用的,终于有点理解了什么是逻辑斯蒂回归和AUC曲线。

参考学习资料:https://www.r-bloggers.com/how-to-perform-a-logistic-regression-in-r/

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 195,898评论 5 462
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 82,401评论 2 373
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 143,058评论 0 325
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,539评论 1 267
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,382评论 5 358
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,319评论 1 273
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,706评论 3 386
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,370评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,664评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,715评论 2 312
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,476评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,326评论 3 313
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,730评论 3 299
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,003评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,275评论 1 251
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,683评论 2 342
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,877评论 2 335

推荐阅读更多精彩内容