R语言实现逻辑回归模型

2021-04-09 11:07:04 浏览数 (1)

首先,本章节使用到的数据集是ISLR包中的Default数据集,数据包含客户信息的模拟数据集。这里的目的是预测哪些客户将拖欠他们的信用卡债务,这个数据集有1w条数据,3个特征:

代码语言:javascript复制
library("ISLR")
library("tibble")
as_tibble(Default)
## # A tibble: 10,000 x 4
##    default student balance income
##    <fct>   <fct>     <dbl>  <dbl>
##  1 No      No         730. 44362.
##  2 No      Yes        817. 12106.
##  3 No      No        1074. 31767.

这里需要分析的是学生身份,信用卡余额,收入这三个特征对违约进行适当的一个分类。

数据探索

拿到数据的第一步还是需要对于数据进行了解。了解的方面包括数据的质量,数据的分布,以及数据之间的关系。密度图可用于识别预测变量相对于彼此的分布以及响应变量,使用ggplot2绘制关于balance特征密度直方图,如图1。

代码语言:javascript复制
library(ggplot2)
ggplot(data = Default,aes(x = balance,color = default)) geom_density() labs(title = "default with balance")

图1 balance的分布

图1描述的违约,不违约两种情况下信用卡余额的分布,从图中可以看出,这两种情况下收入的分布是不一样的。说明这两群人的信用卡余额是不一样的,也就是说,收入对于是否是违约者有很好的区分能力。然后绘制income特征的密度直方图,结果如图2。

代码语言:javascript复制
ggplot(data = Default,aes(x = income,color = default)) geom_density() labs(title = "default with income")

图2 income的分布

从图2中,观察到的是否违约两个群体间的收入分布差异不大,可能认为收入对我们的模型不是特别有用。另一方面,对于是否违约,平均收入在1400左右的值上分布似乎存在很大差异。绘制incomestudent特征的密度直方图,结果如图3。

代码语言:javascript复制
library(ggplot2)
ggplot(data = Default,aes(x = income,color = student)) geom_density() labs(title = "student with income")

图3 income与是否是学生之间的关系

从图3中还可以观察收入和平衡与其他预测学生的分布 这些密度图表明学生的收入远低于其他人口。

构建逻辑回归模型

虽然我们可以做更多的数据探索,但是对于这份数据已经足够了,接下来可以创建逻辑回归模型。为了实现良好的建模实践,将创建训练和测试拆分,以避免在执行回归时过度拟合,下面的代码首先划分了数据集合,一半的数据集为训练集合,一般的结合为测试集合,然后构建逻辑回归模型,使用的是glm构建逻辑回归模型,逻辑回归模型中,使用default特征作为因变量,数据集中所有的其他特征作为自变量。

代码语言:javascript复制
# Split into train/test splits first.
set.seed(42)
default_idx <- sample(nrow(Default), ceiling(nrow(Default) / 2))
default_trn <-  Default[default_idx, ]
default_tst <- Default[-default_idx, ]

# Create the model.
model_glm <- glm(default ~ ., data = default_trn, family = "binomial")

创建逻辑回归模型应该与创建线性回归模型非常相似。但是,我们使用glm()代替lm()。另请注意,我们必须为二进制分类上指定family =“binomial”。(实际上,用family =“gaussian”调用glm()将等同于lm()然后使用summary函数用于查看逻辑回归的详细信息。

代码语言:javascript复制
summary(model_glm)
## 
## Call:
## glm(formula = default ~ ., family = "binomial", data = default_trn)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4137  -0.1496  -0.0596  -0.0214   3.7295  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.026e 01  6.852e-01 -14.976  < 2e-16 ***
## studentYes  -1.010e 00  3.248e-01  -3.109  0.00188 ** 
## balance      5.663e-03  3.252e-04  17.412  < 2e-16 ***
## income      -8.386e-06  1.139e-05  -0.736  0.46152    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1470.42  on 4999  degrees of freedom
## Residual deviance:  813.19  on 4996  degrees of freedom
## AIC: 821.19
## 
## Number of Fisher Scoring iterations: 8

与为线性回归模型计算的summary结果一样,将获得有关残差相关的信息,以及预测变量的显着性估计,logisitic回归框架中p值的解释与线性回归模型的p值相同。

但请注意,逻辑回归模型得到的是z值而不是t值。在没有过多地理解这种差异的理论的情况下,应该理解这个值的这个含义类似于t值的含义。另外,从结果中可以看到看到Null偏差(Null deviance),AIC和Fisher Scoring迭代次数,而不是剩余标准误差,Multipe R平方,调整R平方和F统计量。

summary() 调用生成的逻辑回归诊断值通常不直接用于解释模型的“拟合优度”。

在进行任何预测之前,让我们用summary()简要检查模型。除其他外,重要的是要看看我们的模型估计了哪些系数值。

逻辑回归进行预测

但是,在更仔细地研究更适合于逻辑回归的模型诊断之前,首先应该了解如何使用带有glm()的predict()函数。为了返回概率,我们必须指定type =“response”。

代码语言:javascript复制
head(predict(model_glm, type = "response"))
##         9149         9370         2861         8302         6415 
## 9.572703e-04 4.550820e-01 9.532154e-03 3.281078e-05 1.214581e-04 
##         5189 
## 2.968213e-04

如前所述,这些预测值是可能性,而不是分类。我们必须“手动”将概率转换为分类。传统上,诸如0.5的中点值用于“分类”概率。(这实际上等于指定type =“link”并使用阈值0),使用ifesle将所有预测值大于0.5的结果记为Yes,小于等于0.5的结果记为No。

代码语言:javascript复制
trn_pred <- ifelse(predict(model_glm, type = "response") > 0.5, "Yes", "No")
head(trn_pred)
## 9149 9370 2861 8302 6415 5189 
## "No" "No" "No" "No" "No" "No"

逻辑回归模型评估

评估分类模型最常见的事情可能是使用交叉表将实际响应值与预测响应值进行比较,交叉表通常称为混淆矩阵。可以使用basetable()函数生成此矩阵。

代码语言:javascript复制
trn_tab <- table(predicted = trn_pred, actual = default_trn$default)
trn_tab
##          actual
## predicted   No  Yes
##       No  4814  121
##       Yes   18   47
# Making predictions on the test set.
tst_pred <- ifelse(predict(model_glm, newdata = default_tst, type = "response") > 0.5, "Yes", "No")
tst_tab <- table(predicted = tst_pred, actual = default_tst$default)
tst_tab
##          actual
## predicted   No  Yes
##       No  4815  111
##       Yes   20   54

也许不足为奇的是,评估逻辑回归模型的最常见指标是错误率和准确度(这只是错误率的加性倒数),可以直接从confustion矩阵计算这些指标,下面编写了一个函数,用于计算模型的错误率。

代码语言:javascript复制
calc_class_err <- function(actual, predicted) {
  mean(actual != predicted)
}
calc_class_err(actual = default_trn$default, predicted = trn_pred)
## [1] 0.0278
calc_class_err(actual = default_tst$default, predicted = tst_pred)
## [1] 0.0262

表1 混淆矩阵

现在,我们现在更详细地考虑混淆矩阵。名称真阳性(TP),真阴性(TN),假阳性(FP)和假阴性(FN)通常用于参考燃烧矩阵的四个细胞。

text {Sens}=text {TruePositiveRate}=frac{T P}{P}=frac{T P}{T P F N}
text { Spec }=text { TrueNegativeRate }=frac{T N}{N}=frac{T N}{T N F P}
text {Prev}=frac{P}{text {Totalobs}}=frac{T P F N}{text {Totalobs}}

从混淆矩阵导出诸如灵敏度,特异性和普遍性的度量的计算。这些(和其他)度量的重要性取决于数据的性质(例如,如果认为数据难以预测,则较低的值可能是可接受的),以及对错误分类类型的容忍度。例如,我们可能希望偏向我们对默认值进行分类的预测,以便我们更有可能在未发生默认值时预测默认值。我们必须仔细确定我们是否要优先考虑敏感性或特异性。

我们可以使用caret包中的confusionMatrix()函数轻松获得灵敏度,特异性等值。

代码语言:javascript复制
library("caret")
## Loading required package: lattice
confusionMatrix(trn_tab, positive = "Yes")
## Confusion Matrix and Statistics
## 
##          actual
## predicted   No  Yes
##       No  4814  121
##       Yes   18   47
##                                           
##                Accuracy : 0.9722          
##                  95% CI : (0.9673, 0.9766)
##     No Information Rate : 0.9664          
##     P-Value [Acc > NIR] : 0.01101         
##                                           
##                   Kappa : 0.392           
##  Mcnemar's Test P-Value : < 2e-16         
##                                           
##             Sensitivity : 0.2798          
##             Specificity : 0.9963          
##          Pos Pred Value : 0.7231   
##          Neg Pred Value : 0.9755          
##              Prevalence : 0.0336          
##          Detection Rate : 0.0094          
##    Detection Prevalence : 0.0130          
##       Balanced Accuracy : 0.6380          
##                                           
##        'Positive' Class : Yes             
## 

现在,让我们考虑另一个与混淆矩阵无关的度量。还记得我们选择0.5作为分类门槛的地方吗?我们怎么知道0.5值是准确度的“最佳”值。实际上,其他门槛值可能更好(如果所有模型假设都为真并且样本量相当大,则0.5将倾向于最佳值)。

ROC曲线说明了所有可能的门槛值的灵敏度和特异性。我们可以使用pROC包中的roc()函数为的预测生成ROC曲线,roc()函数的第一个参数是数据集的真实标签,第二个参数是模型的预测结果,第三个参数plot需要输入一个逻辑值,用以表明是否需要绘制ROC曲线图。图3是模型的ROC曲线。

代码语言:javascript复制
library("pROC")
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
test_prob <- predict(model_glm, newdata = default_tst, type = "response")
test_roc <- roc(default_tst$default , test_prob, plot = TRUE, print.auc = TRUE)

图4 ROC曲线图

代码语言:javascript复制
as.numeric(test_roc$auc)
## [1] 0.9519489

一般来说,希望曲线靠向左边界和上边界(表明高灵敏度和特异性)。AUC(曲线下面积)用于量化ROC的轮廓,从图4中可以看到,AUC的值为0.952,模型效果很不错。

注: 本文选自于清华大学出版社出版的《深入浅出R语言数据分析》一书的小节,略有改动。经出版社授权刊登于此

本推文仅供学习,不支持转载,转载需找出版社授权。

内容简介:《深入浅出R语言数据分析/新时代·技术新未来》首先介绍数据分析的方法论,然后介绍数据分析的相关模型方法,并进一步通过数据分析案例,讲解数据分析的思维、方法及模型实现过程。该书重点介绍R语言在数据分析方面的应用,让读者能够快速地使用R语言进行数据分析、构建模型。全书分为17章,内容包括:使用R语言获取数据、数据分析中的数据处理与数据探索、生存分析、主成分分析、多维缩放、线性回归模型、逻辑回归模型、聚类模型、关联规则、随机森林、支持向量机、神经网络、文本挖掘、社交网络分析,以及关于R语言数据分析的两个延伸内容:H2O机器学习和R语言爬虫。

活动方式: 在本公众号下留言区留言,分享一下你学习R的经历或者其他感受,点赞数最高的2位小伙伴获得 《深入浅出R语言数据分析》 一书,免费包邮哦!截止时间2020年12月10日20点整。

0 人点赞