R语言多元 logistic回归 r语言多分类logistic回归分析

转载

夜行者3号 2023-09-30 11:06:53

文章标签 R语言多元 logistic回归 logit回归模型 R 回归虚拟变量na 决策树 文章分类 R语言后端开发

数据分析师常需要基于一组预测变量预测一个二分类问题，如根据个人信息和财务历史记录预测其是否会还贷等。
有监督学习领域中包含许多可用于分类的方法，如逻辑回归，决策树，随机森林，支持向量机，神经网络等。
本文主要介绍用R实现几种分类模型，采用相同的数据集，因此可以直接比较各个方法的结果，对于模型数学原理不做详细讨论。

主要内容

数据准备
逻辑回归
决策树分类
随机森林集成分类
支持向量机模型
评价分类准确性

1 数据准备

本文数据来源于威斯康星州的乳腺癌数据集(http://archive.ics.uci.edu/ml)，包含699个细针抽吸活检的样本单元，11个变量，最后一个变量(class)即为类别，其中458个为良性，241个恶性。
首先，从数据库UCI数据库抽取数据，并随机分成训练集和验证集，训练集70%，验证集30%。

##安装必备程序包
pkgs"rpart","rpart.plot","party","randomForest","e1071")
install.packages(pkgs,depend=TRUE)
loc"http://archive.ics.uci.edu/ml/machine-learning-da
tabases/"
ds"breast-cancer-wisconsin/breast-cancer-wisconsin.data"
url" ")
##缺失值用?表示
breast",",header=FALSE,na.strings="?")
a"ID","clumpThickness","sizeUniformity","shapeUniformity","maginalAdhesion", "singleEpithelialCellSize","bareNuclei","blandChromation","normalNucleoli","mitosis","class")
names(breast)df1]
df$class2,4),labels=c("benign","malignant"))
set.seed(1234)
train0.7*nrow(df))
df.traindf.validatetable(df.train$class)
table(df.validate$class)

##安装必备程序包
pkgs"rpart","rpart.plot","party","randomForest","e1071")
install.packages(pkgs,depend=TRUE)
loc"http://archive.ics.uci.edu/ml/machine-learning-da
tabases/"
ds"breast-cancer-wisconsin/breast-cancer-wisconsin.data"
url" ")
##缺失值用?表示
breast",",header=FALSE,na.strings="?")
a"ID","clumpThickness","sizeUniformity","shapeUniformity","maginalAdhesion", "singleEpithelialCellSize","bareNuclei","blandChromation","normalNucleoli","mitosis","class")
names(breast)df1]
df$class2,4),labels=c("benign","malignant"))
set.seed(1234)
train0.7*nrow(df))
df.traindf.validatetable(df.train$class)
table(df.validate$class)

显示前五条数据：

>  df[1:5,]
  clumpThickness sizeUniformity shapeUniformity maginalAdhesion singleEpithelialCellSize bareNuclei
1              5              1               1               1                        2          1
2              5              4               4               5                        7         10
3              3              1               1               1                        2          2
4              6              8               8               1                        3          4
5              4              1               1               3                        2          1
  blandChromation normalNucleoli mitosis  class
1               3              1       1 benign
2               3              2       1 benign
3               3              1       1 benign
4               3              7       1 benign
5               3              1       1 benign


>  df[1:5,]
  clumpThickness sizeUniformity shapeUniformity maginalAdhesion singleEpithelialCellSize bareNuclei
1              5              1               1               1                        2          1
2              5              4               4               5                        7         10
3              3              1               1               1                        2          2
4              6              8               8               1                        3          4
5              4              1               1               3                        2          1
  blandChromation normalNucleoli mitosis  class
1               3              1       1 benign
2               3              2       1 benign
3               3              1       1 benign
4               3              7       1 benign
5               3              1       1 benign

2 逻辑回归

逻辑回归是广义线性模型的一种，R中glm()可用于拟合逻辑回归问题，glm()函数自动将分类变量编码为相应的虚拟变量，本文中数据全部为数值变量，因此不需要编码。

fit.logitsummary(fit.logit)
probtype="response")
fit.logitsummary(fit.logit)
#predict()函数默认输出肿瘤为恶性的对数概率，指定参数type="response"
#即可得到预测肿瘤为恶性的概率,概率大于0.5被分为恶性肿瘤，小于0.5良性
probtype="response")
#生成混淆矩阵logit.pred,评估模型预测准确性
logit.pred.5,level=c(FALSE,TRUE),labels=c("begin","malignant"))
logit.perf$class,logit.pred,dnn=c("Actual","Predicted"))
logit.perf

fit.logitsummary(fit.logit)
probtype="response")
fit.logitsummary(fit.logit)
#predict()函数默认输出肿瘤为恶性的对数概率，指定参数type="response"
#即可得到预测肿瘤为恶性的概率,概率大于0.5被分为恶性肿瘤，小于0.5良性
probtype="response")
#生成混淆矩阵logit.pred,评估模型预测准确性
logit.pred.5,level=c(FALSE,TRUE),labels=c("begin","malignant"))
logit.perf$class,logit.pred,dnn=c("Actual","Predicted"))
logit.perf

回归结果以及预测效果如下：

Coefficients:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -12.4412     2.0547  -6.055  1.4e-09 ***
clumpThickness             0.7407     0.2262   3.275  0.00106 ** 
sizeUniformity            -0.0320     0.3399  -0.094  0.92500    
shapeUniformity            0.2073     0.3715   0.558  0.57680    
maginalAdhesion            0.5194     0.1708   3.041  0.00236 ** 
singleEpithelialCellSize  -0.3217     0.2613  -1.231  0.21831    
bareNuclei                 0.5851     0.1881   3.111  0.00187 ** 
blandChromation            0.8599     0.2923   2.942  0.00326 ** 
normalNucleoli             0.4036     0.1828   2.208  0.02725 *  
mitosis                    0.8923     0.3552   2.512  0.01200 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> logit.perf
           Predicted
Actual      begin malignant
  benign      129         6
  malignant     1        69

Coefficients:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -12.4412     2.0547  -6.055  1.4e-09 ***
clumpThickness             0.7407     0.2262   3.275  0.00106 ** 
sizeUniformity            -0.0320     0.3399  -0.094  0.92500    
shapeUniformity            0.2073     0.3715   0.558  0.57680    
maginalAdhesion            0.5194     0.1708   3.041  0.00236 ** 
singleEpithelialCellSize  -0.3217     0.2613  -1.231  0.21831    
bareNuclei                 0.5851     0.1881   3.111  0.00187 ** 
blandChromation            0.8599     0.2923   2.942  0.00326 ** 
normalNucleoli             0.4036     0.1828   2.208  0.02725 *  
mitosis                    0.8923     0.3552   2.512  0.01200 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> logit.perf
           Predicted
Actual      begin malignant
  benign      129         6
  malignant     1        69

从混淆矩阵可以看出，模型正确识别了129个类别为良性的患者和69个类别为恶性的患者，且df.validate数据集有4个样本单元含有缺失值无法判断。
在验证集上，准确率为(129+69)/206=96.1%。

3 决策树

3.1 经典决策树

经典决策树以一个二元变量和一组预测变量为基础建立模型，具体算法如下：
(1)选定一个最佳预测变量将全部样本单元分为两类，实现两类中的纯度最大化(即一类中良性样本尽可能多，另一类中恶性尽可能多，纯度主要用信息熵和基尼系数表示)。如果预测变量连续，则选定一个分割点进行分类，使得两类纯度最大化；如果预测变量为分类变量，则对各类别进行合并再分类。
(2)对每一个子类别继续执行步骤(1)
(3)重复(1)~(2),直到子类别中所含的样本单元过少，或者没有分类能将不纯度下降到一个给定阈值以下。
(4)对任一样本执行决策树，得到其终端节点(terminal node)，即可根据模型预测其所属类别。
其中熵定义为：
基尼系数：
当接近0或1时，D和G都接近于0，纯度最大；接近0.5时，纯度最小，信息最混乱。
注：上述算法可能会得到一颗过大的树，可采用10折交叉验证选择预测误差最小的树，即剪枝。
R中rpart包的rpart()函数用于构造决策树，prune()函数用于剪枝。

##经典决策树
library(rpart)
set.seed(1234)#设置随机种子
dtree"class",parms=list(split="information"))#创建决策树
dtree$cptable#返回不同大小树对应的预测误差
dtree.pruned#剪枝
library(rpart.plot)
prp(dtree.pruned,type=2,extra=104,fallen.leaves=TRUE,main="Decision Tree")画出树
dtree.predtype="class")
dtree.perf$class,dtree.pred,dnn=c("Actual","Predicted"))
dtree.perf#得出预测矩阵

##经典决策树
library(rpart)
set.seed(1234)#设置随机种子
dtree"class",parms=list(split="information"))#创建决策树
dtree$cptable#返回不同大小树对应的预测误差
dtree.pruned#剪枝
library(rpart.plot)
prp(dtree.pruned,type=2,extra=104,fallen.leaves=TRUE,main="Decision Tree")画出树
dtree.predtype="class")
dtree.perf$class,dtree.pred,dnn=c("Actual","Predicted"))
dtree.perf#得出预测矩阵

首先查看cptable内容:

> dtree$cptable
    CP         nsplit  rel error    xerror  xstd
1 0.81764706      0 1.00000000 1.0000000 0.06194645
2 0.04117647      1 0.18235294 0.1823529 0.03169642
3 0.01764706      3 0.10000000 0.1588235 0.02970979
4 0.01000000      4 0.08235294 0.1235294 0.02637116

> dtree$cptable
    CP         nsplit  rel error    xerror  xstd
1 0.81764706      0 1.00000000 1.0000000 0.06194645
2 0.04117647      1 0.18235294 0.1823529 0.03169642
3 0.01764706      3 0.10000000 0.1588235 0.02970979
4 0.01000000      4 0.08235294 0.1235294 0.02637116

cptable包含不同大小的树对应的预测误差，可以用来辅助设定最终的树的大小。复杂度参数(cp)用于惩罚过大的树，树的大小即分支数(nsplit),有n个分支的树将有n+1个终端节点；real error为训练集中各种树对应的误差，xerror是交叉验证误差。
查看dtree.perf预测矩阵:

> dtree.perf
           Predicted
Actual      benign malignant
  benign       129        10
  malignant      2        69

> dtree.perf
           Predicted
Actual      benign malignant
  benign       129        10
  malignant      2        69

准确率(129+69)/210=94.2%，与逻辑回归不同的是，验证集中的210个样本都可以由最终决策树分类，但是对于水平数很多或缺失值很多的预测变量，决策树可能有偏。

下面是最终生成的决策树：

R语言多元 logistic回归 r语言多分类logistic回归分析_R语言多元 logistic回归

上图中，每一个节点处都有对应类别的概率以及样本单元的占比，比如最顶端的根节点，100%的样本单元中，有85%的单元良性，35%恶性。从树的顶端开始，若满足条件sizeUniformity<3,则从左枝往下，否则右枝往下重复这个过程直到碰到一个终端节点为止，该终端节点即为这一观测点的所属类别。

3.2 条件决策树

条件决策树与经典决策树类似，但变量与分割的选取是基于显著性检验的，而不是纯净度或同质性一类的度量，用到的显著性检验是置换检验。条件推断树算法如下：
(1)对输出变量与每个预测变量间的关系计算p值；
(2)选取p值最小的变量；
(3)在因变量与被选变量间尝试所有可能的二元分割(通过排列检验)，并选取最显著的分割；
(4)将数据集分成两群，并对每个子群重复上述步骤；
(5)重复直至所有分割都不显著或已达最小节点为止。
R中party包的ctree()函数可用于生成条件决策树：

##条件决策树####
library(party)
fit.ctree#plot()画出决策树，参数type="extended"表示展示条件推断树
plot(fit.ctree,type="extended",main="Conditional Inference Tree")
ctree.predtype="response")
ctree.pref$class,ctree.pred,dnn=c("Actual","Predicted"))
ctree.pref
> ctree.pref
           Predicted
Actual      benign malignant
  benign       131         8
  malignant      4        67

##条件决策树####
library(party)
fit.ctree#plot()画出决策树，参数type="extended"表示展示条件推断树
plot(fit.ctree,type="extended",main="Conditional Inference Tree")
ctree.predtype="response")
ctree.pref$class,ctree.pred,dnn=c("Actual","Predicted"))
ctree.pref
> ctree.pref
           Predicted
Actual      benign malignant
  benign       131         8
  malignant      4        67

上述条件决策树准确率为(131+67)/210=94.2%，下图展示了一棵条件推断树，每个节点中的阴影区域为这个节点处对应的恶性肿瘤比例。

R语言多元 logistic回归 r语言多分类logistic回归分析_logit回归模型_02

4 随机森林

随机森林是一种组成式的有监督学习方法，在随机森林中，同时生成多个预测模型，并将模型的结果汇总以提升准确率。
随机森林算法涉及对样本单元和变量进行抽样，从而生成大量决策树，对每个样本单元来说，所有决策树依次对其进行分类，所有决策树预测类别中的众数即为随机森林所预测的这一样本单元的类别。
对于N个样本单元，M个变量，其算法如下：

从训练集中随机有放回地抽取N个单元，生成大量决策树；
在每一个节点随机抽取m
完整生成所有决策树，无需剪枝(最小节点为1);
终端节点的所属类别由节点对应的众数类别决定；
对于新的观测点，用所有的树对其进行分类，其类别由多数决定原则生成。注：生成树没有用到的样本点所对应的类别可由生成的树估计，与真实类别比较即可得到预测误差(也称为袋外预测，out-of-bag，OOB)，这一特性在无法获得验证集时，是一大优势。
R中randomForest包中的randomForest()函数可用于生成随机森林，函数默认生成500棵树，每个节点处抽取个变量，最小节点为1。

####随机森林###
library("randomForest")
set.seed(1234)
fit.forest#na.action=na.roughfix：可将变量中缺失值换成对应的中位数，importance=TRUE度量变量重要性
fit.forest
importance(fit.forest,type=2)#输出变量重要性
##type=2:得到的变量相对重要性，就是分割该变量时节点不纯度的下降总量对所有树取平均
#节点不纯度由Gini系数确定
forest.predforest.pref$class,forest.pred,dnn=c("Actual","Predicted"))
forest.pref

plot(fit.forest)
> forest.pref
           Predicted
Actual      benign malignant
  benign       128         7
  malignant      2        68


####随机森林###
library("randomForest")
set.seed(1234)
fit.forest#na.action=na.roughfix：可将变量中缺失值换成对应的中位数，importance=TRUE度量变量重要性
fit.forest
importance(fit.forest,type=2)#输出变量重要性
##type=2:得到的变量相对重要性，就是分割该变量时节点不纯度的下降总量对所有树取平均
#节点不纯度由Gini系数确定
forest.predforest.pref$class,forest.pred,dnn=c("Actual","Predicted"))
forest.pref

plot(fit.forest)
> forest.pref
           Predicted
Actual      benign malignant
  benign       128         7
  malignant      2        68

通过对验证集的样本单元进行分类，得到准确率为(128+68)/205=95.6%，结合下面残差图，在500棵树时，误差已经足够小。

R语言多元 logistic回归 r语言多分类logistic回归分析_R语言多元 logistic回归_03

5 支持向量机

支持向量机(SVM)是一类用于分类和回归的有监督学习模型。SVM有两个优势，可输出比较准确的模型和基于较优雅的数学理论。

SVM旨在在多维空间中找到一个能将全部样本单元分成两类的最优平面，这一平面应使两类中距离最近的点的间距(margin)尽可能大，在间距边界上的点称为支持向量(support vector,它们决定间距)，分割的超平面位于间距的中间。

对于一个N维空间(即N个变量)来说，最优超平面为N-1维。当变量数为2时，曲面是一条直线；变量数为3时，曲面是一个平面。

下图是一个SVM分二维到三维分类展示：

R语言多元 logistic回归 r语言多分类logistic回归分析_logit回归模型_04

SVM的数学解释比较复杂，这里不做详细讨论。

R中e1071包中的svm()函数可实现SVM分类:

#####SVM###
library(e1071)
set.seed(1234)
fit.svmfit.svm
svm.predsvm.pref$class,svm.pred,dnn=c("Actual","Predicted"))
> svm.pref
           Predicted
Actual      benign malignant
  benign       126         9
  malignant     0         70       


#####SVM###
library(e1071)
set.seed(1234)
fit.svmfit.svm
svm.predsvm.pref$class,svm.pred,dnn=c("Actual","Predicted"))
> svm.pref
           Predicted
Actual      benign malignant
  benign       126         9
  malignant     0         70

svm()函数默认将变量标准化，由上述结果，svm()的预测准确率不错，为(126+70)/205=95.6%，与随机森林不同的是，SVM在预测新样本单元时不允许有缺失值出现。

6 评价分类的准确性

以上我们介绍了几种有监督机器学习对细针抽活检细胞进行分类，但如何从中选出最准确的方法呢？最常用的是前文提到的准确率(accuary)，即分类器是否正确地划分样本单元，不过这一指标有其不足之处，比如：
假设我们现在需要判别一个人是否患有精神分裂症。精神分裂症是一种极少见的的生理障碍，人群中的患病率约为1%。如果一种分类方法将全部人都判为未患病，则这一分类器的准确率达到99%，但它会把所有患精神病的人都识别为健康人。从这个角度看，它显然不是一个好的分类器。因此，在准确率之外，我们还应该思考以下问题：

患有精神病分裂症的人中有多大比例成功鉴别？
未患病的人中有多大比例被成功鉴别？
如果一个人被鉴别为精神分裂症患者，这个判别有多大概率是准确的呢？
如果一个人被鉴别为未患病，这个判别又有多大概率是准确的呢？上述问题涉及一个分类器的敏感度(sensitivity)，特异性(sensitivity),正例命中率(positive predicitive power),负例命中率(negative predicitive)，具体解释如下表：然后我们定义一个函数计算这几个统计量：

performancefunction(table,n=2){
  if(!all(dim(table)==c(2,2)))
    stop("Must be a 2*2 table")
  tn=table[1,1]
  fp=table[1,2]
  fn=table[2,1]
  tp=table[2,2]
  sensitivity=tp/(tp+fn)
  specificity=tn/(tn+fp)
  ppp=tp/(tp+fp)
  npp=tn/(tn+fn)
  hitrate = (tp+tn)/(tp+tn+fp+fn)
  result = paste("Sensitivity=",round(sensitivity,n),
                 "\nSpecificity=",round(specificity,n),
                 "\nNegative Predictive Value=",round(ppp,n),
                 "\nPositive Predictive Value=",round(npp,n),
                 "\nAccuracy=",round(hitrate,n),"\n",sep=""
                 )
  cat(result)
}

performancefunction(table,n=2){
  if(!all(dim(table)==c(2,2)))
    stop("Must be a 2*2 table")
  tn=table[1,1]
  fp=table[1,2]
  fn=table[2,1]
  tp=table[2,2]
  sensitivity=tp/(tp+fn)
  specificity=tn/(tn+fp)
  ppp=tp/(tp+fp)
  npp=tn/(tn+fn)
  hitrate = (tp+tn)/(tp+tn+fp+fn)
  result = paste("Sensitivity=",round(sensitivity,n),
                 "\nSpecificity=",round(specificity,n),
                 "\nNegative Predictive Value=",round(ppp,n),
                 "\nPositive Predictive Value=",round(npp,n),
                 "\nAccuracy=",round(hitrate,n),"\n",sep=""
                 )
  cat(result)
}

用performance()函数比较上述五种分类器的效果：

> performance(logit.perf)
Sensitivity=0.99
Specificity=0.96
Negative Predictive Value=0.92
Positive Predictive Value=0.99
Accuracy=0.97
> performance(dtree.perf)
Sensitivity=0.97
Specificity=0.93
Negative Predictive Value=0.87
Positive Predictive Value=0.98
Accuracy=0.94
> performance(ctree.pref)
Sensitivity=0.94
Specificity=0.94
Negative Predictive Value=0.89
Positive Predictive Value=0.97
Accuracy=0.94
> performance(forest.pref)
Sensitivity=0.97
Specificity=0.95
Negative Predictive Value=0.91
Positive Predictive Value=0.98
Accuracy=0.96
> performance(svm.pref)
Sensitivity=1
Specificity=0.93
Negative Predictive Value=0.89
Positive Predictive Value=1
Accuracy=0.96

> performance(logit.perf)
Sensitivity=0.99
Specificity=0.96
Negative Predictive Value=0.92
Positive Predictive Value=0.99
Accuracy=0.97
> performance(dtree.perf)
Sensitivity=0.97
Specificity=0.93
Negative Predictive Value=0.87
Positive Predictive Value=0.98
Accuracy=0.94
> performance(ctree.pref)
Sensitivity=0.94
Specificity=0.94
Negative Predictive Value=0.89
Positive Predictive Value=0.97
Accuracy=0.94
> performance(forest.pref)
Sensitivity=0.97
Specificity=0.95
Negative Predictive Value=0.91
Positive Predictive Value=0.98
Accuracy=0.96
> performance(svm.pref)
Sensitivity=1
Specificity=0.93
Negative Predictive Value=0.89
Positive Predictive Value=1
Accuracy=0.96

在这个案例中，这些分类都表现得相当不错，但是实际问题中并不总是这样。随机森林相对更好，成功鉴别了97%的恶性样本和95%的良性样本，总体来说预测准确率达到96%。91%被判为恶性组织的样本单元确实是恶性的(即9%正例错误率)，98%被判为良性组织的样本单元确实是良性的(即2%负例错误率)。虽然SVM的精确度也很高，但从特异性和敏感性来看，对于癌症诊断，特异性(成功鉴别恶性样本的概率)更重要。因此，评估分类模型的优劣，需要结合实际问题的研究目的，选择合适的评价指标。

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。