R语言实战(17)——分类
引言:不知道之前一起学习了R语言实战前16章的同学们是不是感觉还有些意犹未尽呢?没错,我们《R语言实战》学习栏目恢复更新啦!本期开始我们将陆续推出本书最后的第17-23章,欢迎小伙伴们继续和我们一起学习R语言的更多功能。
后台回复“R语言实战”即可获取二维码加入R语言实战学习讨论群。
本期我们将学习如何根据一组预测变量(或特征)来预测相对应的二分类结果(无信用风险/有信用风险,心脏病发作/心脏病未发作,是病毒邮件/不是病毒邮件)。
有监督机器学习领域中包含许多可用于分类的方法,如逻辑回归、决策树、随机森林、支持向量机、神经网络等。
安装所需的R包:
通过rpart、rpart.plot和party包来实现决策树模型及其可视化,通过randomForest包拟合随机森林,通过e1071包构造支持向量机,通过R中的基本函数glm()实现逻辑回归。
pkgs <- c("rpart", "rpart.plot", "party",
"randomForest", "e1071")
install.packages(pkgs, depend=TRUE)
17.1 数据准备
本章将用到UCI机器学习数据库(http://archive.ics.uci.edu/ml)中的威斯康星州乳腺癌数据,是一个由逗号分隔的txt文件。本数据集包含699个细针抽吸活检的样本单元,其中458个(65.5%)为良性样本单元,241个(34.5%)为恶性样本单元。
loc <- "http://archive.ics.uci.edu/ml/machine-learning-databases/"
ds <- "breast-cancer-wisconsin/breast-cancer-wisconsin.data"
url <- paste(loc, ds, sep="")
#在UCI机器学习数据库中找到该数据集
breast <- read.table(url, sep = ",", header = FALSE, na.strings = "?")
#表中未标明变量名,缺失数据用问号(?)表示
names(breast) <- c("ID", "clumpThickness", "sizeUniformity",
"shapeUniformity", "maginalAdhesion",
"singleEpithelialCellSize", "bareNuclei",
"blandChromatin", "normalNucleoli", "mitosis", "class")
数据集中包含的变量包括:ID、肿块厚度、细胞大小的均匀性、细胞形状的均匀性、边际附着力、单个上皮细胞大小、裸核、乏味染色体、正常核、有丝分裂、类别。
df <- breast[-1]
#第一个变量ID不纳入数据分析
df$class <- factor(df$class, levels=c(2,4),
labels=c("benign", "malignant"))
#变量(类别)即输出变量(编码为良性=2,恶性=4)
set.seed(1234)
train <- sample(nrow(df), 0.7*nrow(df))
df.train <- df[train,]
df.validate <- df[-train,]
#数据从UCI数据库中抽取,并随机分出训练集和验证集,其中训练集中包含489个样本单元(占70%)
table(df.train$class) #训练集
benign malignant
319 170
table(df.validate$class) #验证集
benign malignant
139 71
训练集将用于建立逻辑回归、决策树、条件推断树、随机森林、支持向量机等分类模型,测试集用于评估各个模型的有效性。
17.2 逻辑回归
使用glm()以类别为响应变量,其余变量为预测变量拟合逻辑回归模型:
fit.logit <- glm(class~., data=df.train, family=binomial())
summary(fit.logit)
Call:
glm(formula = class ~ ., family = binomial(), data = df.train)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.24605 -0.08012 -0.03110 0.00266 2.11576
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 **
blandChromatin 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
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 621.04 on 477 degrees of freedom
Residual deviance: 52.39 on 468 degrees of freedom
(11 observations deleted due to missingness)
AIC: 72.39
Number of Fisher Scoring iterations: 9
采用基于df.train建立的模型来对df.validate数据集中的样本单元分类:
prob <- predict(fit.logit, df.validate, type="response")
logit.pred <- factor(prob > .5, levels=c(FALSE, TRUE),
labels=c("benign", "malignant"))
给出预测与实际情况对比的交叉表(即混淆矩阵,confusion matrix):
logit.perf <- table(df.validate$class, logit.pred,
dnn=c("Actual", "Predicted"))
logit.perf
Predicted
Actual benign malignant
benign 129 6
malignant 1 69
可以看出,模型正确判别了129个类别为良性的患者和69个类别为恶性的患者。
df.validate数据集中有5个样本单元因包含缺失数据而无法判别。
在验证集上,正确分类的模型(即准确率,accuracy)为(69+129)/205=96.6%。
模型中有三个预测变量(sizeUniformity、shapeUniformity和singleEpithelialCellSize)的系数未通过显著性检验(即p值大于0.1)。
此时,可用逐步逻辑回归生成一个包含更少解释变量的模型:
logit.fit.reduced <- step(fit.logit)
Start: AIC=72.39
class ~ clumpThickness + sizeUniformity + shapeUniformity + maginalAdhesion + singleEpithelialCellSize + bareNuclei + blandChromatin + normalNucleoli + mitosis
Df Deviance AIC
- sizeUniformity 1 52.399 70.399
- shapeUniformity 1 52.686 70.686
- singleEpithelialCellSize 1 53.910 71.910
<none> 52.390 72.390
- normalNucleoli 1 57.465 75.465
- mitosis 1 57.966 75.966
- blandChromatin 1 62.856 80.856
- maginalAdhesion 1 63.044 81.044
- bareNuclei 1 64.982 82.982
- clumpThickness 1 68.323 86.323
Step: AIC=70.4
class ~ clumpThickness + shapeUniformity + maginalAdhesion + singleEpithelialCellSize + bareNuclei + blandChromatin + normalNucleoli + mitosis
Df Deviance AIC
- shapeUniformity 1 52.876 68.876
- singleEpithelialCellSize 1 53.918 69.918
<none> 52.399 70.399
- normalNucleoli 1 57.916 73.916
- mitosis 1 58.024 74.024
- blandChromatin 1 63.272 79.272
- maginalAdhesion 1 63.735 79.735
- bareNuclei 1 64.985 80.985
- clumpThickness 1 68.728 84.728
Step: AIC=68.88
class ~ clumpThickness + maginalAdhesion + singleEpithelialCellSize + bareNuclei + blandChromatin + normalNucleoli + mitosis
Df Deviance AIC
- singleEpithelialCellSize 1 54.154 68.154
<none> 52.876 68.876
- mitosis 1 59.402 73.402
- normalNucleoli 1 60.929 74.929
- blandChromatin 1 64.053 78.053
- maginalAdhesion 1 64.995 78.995
- bareNuclei 1 75.634 89.634
- clumpThickness 1 76.942 90.942
Step: AIC=68.15
class ~ clumpThickness + maginalAdhesion + bareNuclei + blandChromatin + normalNucleoli + mitosis
Df Deviance AIC
<none> 54.154 68.154
- mitosis 1 60.177 72.177
- normalNucleoli 1 60.937 72.937
- blandChromatin 1 64.056 76.056
- maginalAdhesion 1 65.022 77.022
- bareNuclei 1 76.417 88.417
- clumpThickness 1 77.027 89.027
通过增加或移除变量来得到一个更小的AIC值,上面提到的三个变量就从最终模型中移除,这种精简后的模型在验证集上的误差相对全变量模型更小。
17.3 决策树
决策树基本思想:对预测变量进行二元分离,从而构造一棵可用于预测新样本单元所属类别的树。
1. 经典决策树
经典决策树以一个二元输出变量(对应威斯康星州乳腺癌数据集中的良性/恶性)和一组预测变量(对应九个细胞特征)为基础。具体算法如下:
(1) 选定一个最佳预测变量将全部样本单元分为两类,实现两类中的纯度最大化(即一类中良性样本单元尽可能多,另一类中恶性样本单元尽可能多)。
如果预测变量为连续变量,则选定一个分割点进行分类,使得两类纯度最大化;
如果预测变量为分类变量,则对各类别进行合并再分类。
(2) 对每一个子类别继续执行步骤(1)。
(3) 重复步骤(1)~(2),直到子类别中所含的样本单元数过少,或者没有分类法能将不纯度下降到一个给定阈值以下。
最终集中的子类别即终端节点(terminal node)。
根据每一个终端节点中样本单元的类别数众数来判别这一终端节点的所属类别。
(4) 对任一样本单元执行决策树,得到其终端节点,即可根据步骤3得到模型预测的所属类别。
上述算法通常会得到一棵过大的树,从而出现过拟合现象。结果就是,对于训练集外单元的分类性能较差。可采用10折交叉验证法选择预测误差最小的树——剪枝。
R中的rpart包支持rpart()函数构造决策树,prune()函数对决策树进行剪枝。
使用rpart()函数创建分类决策树:
library(rpart)
set.seed(1234)
dtree <- rpart(class ~ .,data = df.train, method = "class",
parms = list(split="information"))
#生成决策树,此时所得的树可能过大,需要剪枝
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个终端节点;
rel error栏即训练集中各种树对应的误差;
交叉验证误差(xerror)即基于训练样本所得的10折交叉验证误差;
xstd栏为交叉验证误差的标准差。
对于所有交叉验证误差在最小交叉验证误差一个标准差范围内的树,最小的树即最优的树。
plotcp(dtree)
#画出交叉验证误差与复杂度参数的关系图
图17-1 复杂度参数与交叉验证误差
虚线是基于一个标准差准则得到的上限(0.12+1×0.0264=0.15)。从图像来看,应选择虚线下最左侧cp值对应的树。
最小的交叉验证误差为0.12,标准误差为0.0264,则最优的树为交叉验证误差在0.12±0.0264(0.09和0.15)之间的树。
由cptable表可知,五个终端节点(即四次分割)的树满足要求(交叉验证误差为0.123 52);根据图17-1也可以选得最优树,即四次分割(五个节点)对应的树。
prune()函数根据复杂度参数剪掉最不重要的枝:
dtree.pruned <- prune(dtree, cp=.01)
从cptable中可以看到,三次分割对应的复杂度参数为0.01,prune(dtree, cp=0.0125)可得到一个理想大小的树。
prp()函数可用于画出最终的决策树:
library(rpart.plot)
prp(dtree.pruned,type = 2,extra = 104,
fallen.leaves = TRUE, main="Decision Tree")
图17-2 用剪枝后的传统决策树预测癌症状态
当观测点到达终端节点时,分类结束。每一个节点处都有对应类别的概率以及样本单元的占比。
type=2可画出每个节点下分割的标签
extra=104可画出每一类的概率以及每个节点处的样本占比
fallen.leaves=TRUE可在图的底端显示终端节点
predict()函数对验证集中的观测点分类:
dtree.pred <- predict(dtree.pruned, df.validate, type="class")
dtree.perf <- table(df.validate$class, dtree.pred,
dnn = c("Actual", "Predicted"))
dtree.perf #实际类别与预测类别的交叉表
Predicted
Actual benign malignant
benign 129 10
malignant 2 69
验证集中的准确率达到了94%。与逻辑回归不同的是,验证集中的210个样本单元都可由最终树来分类。
2. 条件推断树
条件推断树(conditional inference tree)与传统决策树类似,但变量和分割的选取是基于显著性检验的,而不是纯净度或同质性一类的度量。显著性检验是置换检验。算法如下:
(1) 对输出变量与每个预测变量间的关系计算p值。
(2) 选取p值最小的变量。
(3) 在因变量与被选中的变量间尝试所有可能的二元分割(通过排列检验),并选取最显著的分割。
(4) 将数据集分成两群,并对每个子群重复上述步骤。
(5) 重复直至所有分割都不显著或已到达最小节点为止。
条件推断树可由party包中的ctree()函数获得。
使用ctree()函数创建条件推断树:
library(party)
fit.ctree <- ctree(class~., data = df.train)
plot(fit.ctree, main="Conditional Inference Tree")
ctree.pred <- predict(fit.ctree, df.validate, type="response")
ctree.perf <- table(df.validate$class, ctree.pred,
dnn=c("Actual", "Predicted"))
ctree.perf
Predicted
Actual benign malignant
benign 131 8
malignant 4 67
图17-3 乳腺癌数据的条件推断树
对于条件推断树来说,剪枝不是必需的。每个节点中的阴影区域代表这个节点对应的恶性肿瘤比例。
17.4 随机森林
随机森林(random forest)是一种组成式的有监督学习方法,可以同时生成多个预测模型,并将模型的结果汇总以提升分类准确率。
随机森林的算法涉及对样本单元和变量进行抽样,从而生成大量决策树。对每个样本单元来说,所有决策树依次对其进行分类。所有决策树预测类别中的众数类别即为随机森林所预测的这一样本单元的类别。
假设训练集中共有N个样本单元,M个变量,则随机森林算法如下:
(1) 从训练集中随机有放回地抽取N个样本单元,生成大量决策树。
(2) 在每一个节点随机抽取m(3) 完整生成所有决策树,无需剪枝(最小节点为1)。
(4) 终端节点的所属类别由节点对应的众数类别决定。
(5) 对于新的观测点,用所有的树对其进行分类,其类别由多数决定原则生成。
生成树时没有用到的样本点所对应的类别可由生成的树估计,与其真实类别比较即可得到袋外预测(out-of-bag,OOB)误差。无法获得验证集时,这是随机森林的一大优势。
randomForest包中的randomForest()函数可用于生成随机森林。函数默认生成500棵树,并且默认在每个节点处抽取sqrt(M)个变量,最小节点为1。
生成随机森林:
library(randomForest)
set.seed(1234)
fit.forest <- randomForest(class~., data=df.train, na.action=na.roughfix, importance=TRUE)
fit.forest
Call:
randomForest(formula = class ~ ., data = df.train, importance = TRUE, na.action = na.roughfix)
Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 3
OOB estimate of error rate: 2.66%
Confusion matrix:
benign malignant class.error
benign 313 6 0.01880878
malignant 7 163 0.04117647
randomForest()函数从训练集中有放回地随机抽取489个观测点,在每棵树的每个节点随机抽取3个变量,从而生成了500棵传统决策树。
na.action=na.roughfix参数可将数值变量中的缺失值替换成对应列的中位数,类别变量中的缺失值替换成对应列的众数类(若有多个众数则随机选一个)。
设置information=TRUE参数,可得到变量重要性。
importance()函数输出变量重要性:
importance(fit.forest, type=2)
MeanDecreaseGini
clumpThickness 11.382432
sizeUniformity 63.037519
shapeUniformity 49.027649
maginalAdhesion 4.275345
singleEpithelialCellSize 14.504981
bareNuclei 42.736364
blandChromatin 22.484755
normalNucleoli 11.375285
mitosis 1.755382
type=2,得到变量相对重要性为分割该变量时节点不纯度(异质性)的下降总量对所有树取平均。
节点不纯度由Gini系数定义。
本例中,sizeUniformity是最重要的变量,mitosis是最不重要的变量。
通过随机森林算法对验证集中的样本单元进行分类:
forest.pred <- predict(fit.forest, df.validate)
forest.perf <- table(df.validate$class, forest.pred,
dnn=c("Actual", "Predicted"))
forest.perf
Predicted
Actual benign malignant
benign 128 7
malignant 2 68
分类时剔除验证集中有缺失值的单元。对验证集的预测准确率高达95.6%。
随机森林的优点:
分类准确率通常更高
可处理大规模问题(即多样本单元、多变量)
可处理训练集中有大量缺失值的数据
可应对变量远多于样本单元的数据
可计算袋外预测误差(OOB error)、度量变量重要性
随机森林的缺点:
分类方法(此例中相当于500棵决策树)较难理解和表达。
17.5 支持向量机
支持向量机(SVM)是一类可用于分类和回归的有监督机器学习模型。
SVM旨在在多维空间中找到一个能将全部样本单元分成两类的最优平面,这一平面应使两类中距离最近的点的间距(margin)尽可能大,在间距边界上的点被称为支持向量(support vector,它们决定间距),分割的超平面位于间距的中间。
对于一个N维空间(即N个变量)来说,最优超平面(即线性决策面,linear decision surface)为N–1维。
图17-4 线性可分的二分类问题
圆圈和三角形分别代表两个不同类别,间距即两根虚线间的距离。虚线上的点(实心的圆圈和三角形)即支持向量。最优超平面即间距中的黑色实线。
若数据是非线性的,SVM通过核函数将数据投影到高维,使其在高维线性可分。
图17-5 线性不可分的特征
可以想象对图17-5的数据投影,从而将圆圈从纸上分离出来,使其位于三角形上方的平面。一种方法是将二维数据投影到三维空间:
这样就可以用一张硬纸片将三角形与圆圈分开(一个二维平面变成了一个三维空间)。
SVM可以通过R中kernlab包的ksvm()函数和e1071包中的svm()函数实现。
通过svm()函数对威斯康星州乳腺癌数据建立SVM模型:
library(e1071)
set.seed(1234)
fit.svm <- svm(class~., data=df.train)
fit.svm
Call:
svm(formula = class ~ ., data = df.train)
Parameters:
SVM-Type: C-classification
SVM-Kernel: radial
cost: 1
Number of Support Vectors: 74
svm.pred <- predict(fit.svm, na.omit(df.validate))
svm.perf <- table(na.omit(df.validate)$class,
svm.pred, dnn=c("Actual", "Predicted"))
svm.perf
Predicted
Actual benign malignant
benign 126 9
malignant 0 70
由于方差较大的预测变量通常对SVM的生成影响更大,svm()函数默认在生成模型前对每个变量标准化,使其均值为0、标准差为1。
从结果来看,SVM的预测准确率为95.6%。与随机森林算法不同的是,SVM在预测新样本单元时不允许有缺失值出现。
选择调和参数
svm()函数默认通过径向基函数(Radial Basis Function,RBF)将样本单元投射到高维空间。
RBF核是一种非线性投影,可以应对类别标签与预测变量间的非线性关系。
用带RBF核的SVM拟合样本时,可能影响最终结果的两个参数:
gamma,核函数的参数,控制分割超平面的形状。
gamma越大,通常导致支持向量越多。也可将gamma看作控制训练样本“到达范围”的参数,即gamma越大意味着训练样本到达范围越广,而越小则意味着到达范围越窄。gamma必须大于0。cost,成本参数代表犯错的成本。
一个较大的成本意味着模型对误差的惩罚更大,从而将生成一个更复杂的分类边界,对应的训练集中的误差也会更小,可能存在过拟合问题,即对新样本单元的预测误差可能很大。较小的成本意味着分类边界更平滑,但可能会导致欠拟合。成本参数也恒为正。
svm()函数默认设置gamma为预测变量个数的倒数,成本参数为1。
可以通过tune.svm()对每个参数设置一个候选范围,tune.svm()函数对每一个参数组合生成一个SVM模型,并输出在每一个参数组合上的表现。
对不同的gamma和成本拟合一个带RBF核的SVM模型:
set.seed(1234)
tuned <- tune.svm(class~., data=df.train,
gamma=10^(-6:1),
cost=10^(-10:10))
tuned
Parameter tuning of ‘svm’:
- sampling method: 10-fold cross validation
- best parameters:
gamma cost
0.01 1
- best performance: 0.02740302
尝试八个不同的gamma(从0.000 001到10)以及21个成本参数(从0.01到10^10),共拟合了168(8×21)个模型,并比较了其结果。训练集中10折交叉验证误差最小的模型所对应的参数为gamm=0.01,成本参数为1。
对全部训练样本拟合出新的SVM模型,对验证集中的样本单元进行预测:
fit.svm <- svm(class~., data=df.train, gamma=.01, cost=1)
svm.pred <- predict(fit.svm, na.omit(df.validate))
svm.perf <- table(na.omit(df.validate)$class,
svm.pred, dnn=c("Actual", "Predicted"))
svm.perf
Predicted
Actual benign malignant
benign 128 7
malignant 0 70
在本例中,调和后的模型轻微减少了错分个数(从7减少到6)。一般来说,为SVM模型选取调和参数通常可以得到更好的结果。
SVM的优点:
适用面比较广
可以应用于变量数远多于样本单元数的问题
SVM的缺点:
分类准则比较难以理解和表述
对大量样本建模时不如随机森林
17.6 选择预测效果最好的解
表17-1 预测准确性度量
统计量 | 解释 |
---|---|
敏感度sensitivity | 正类的样本单元被成功预测的概率,也叫正例覆盖率(true positive)或召回率(recall) |
特异性specificity | 负类的样本单元被成功预测的概率,也叫负例覆盖率(true negative) |
正例命中率positive predictive power | 被预测为正类的样本单元中,预测正确的样本单元占比,也叫精确度(precision) |
负例命中率negative predictive power | 被预测为负类的样本单元中,预测正确的样本单元占比 |
准确率accuracy | 被正确分类的样本单元所占比重,也叫ACC |
评估二分类准确性:
performance <- function(table, n=2){
if(!all(dim(table) == c(2,2)))
stop("Must be a 2 x 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),
"\nPositive Predictive Value = ", round(ppp, n),
"\nNegative Predictive Value = ", round(npp, n),
"\nAccuracy = ", round(hitrate, n), "\n", sep="") #输出结果
cat(result)
}
给定真值(行)和预测值(列),performance()函数可给出这五个准确性度量的值。
函数首先提取出负类中正确的个数(良性组织被判别为良性)、负类中错分的个数(恶性组织被判为良性)、正类中错分的个数(良性组织被判为恶性)、正类中正确的个数(恶性组织被判为恶性)。
计算敏感度、特性性、正例命中率、负例命中率和准确率。
显示规范后的结果。
乳腺癌数据分类器的性能:
performance(logit.perf)#逻辑回归
Sensitivity = 0.99
Specificity = 0.96
Positive Predictive Value = 0.92
Negative Predictive Value = 0.99
Accuracy = 0.97
performance(dtree.perf)#传统决策树
Sensitivity = 0.97
Specificity = 0.93
Positive Predictive Value = 0.87
Negative Predictive Value = 0.98
Accuracy = 0.94
performance(ctree.perf)#条件推断树
Sensitivity = 0.94
Specificity = 0.94
Positive Predictive Value = 0.89
Negative Predictive Value = 0.97
Accuracy = 0.94
performance(forest.perf)#随机森林
Sensitivity = 0.97
Specificity = 0.95
Positive Predictive Value = 0.91
Negative Predictive Value = 0.98
Accuracy = 0.96
performance(svm.perf)#支持向量机
Sensitivity = 1
Specificity = 0.95
Positive Predictive Value = 0.91
Negative Predictive Value = 1
Accuracy = 0.97
17.7 小结
本章介绍了一系列用于二分类的机器学习方法,包括逻辑回归分类方法、传统决策树、条件推断树、集成性的随机森林以及支持向量机。下一期我们将学习如何处理缺失数据。
后台回复“R语言实战”即可获取二维码加入R语言实战学习讨论群。