logo

文章是学习《R语言实战》的成果整理。文中数据来源为威斯康星州乳腺癌数据集,本数据包含699个样本,11个变量。可在 UCI机器学习数据库中找到。

乳腺癌数据准备

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
loc = "https://archive.ics.uci.edu/ml/machine-learning-databases/"
ds = "breast-cancer-wisconsin/breast-cancer-wisconsin.data"
url = paste(loc, ds, sep="")
breast = read.table(url, sep=",", header=FALSE, na.strings="?")
names(breast) = c("ID", "clumpThickness", "sizeUniformity","shapeUniformity",
                  "maginalAdhesion","singleEpithelialCellSize", "bareNuclei",
                  "blandChromatin", "normalNucleoli", "mitosis", "class")
df = breast[-1]  # 去掉ID列
df$class = factor(df$class, levels=c(2,4),
                  labels=c("benign", "malignant"))

set.seed(1234)  # 设置随机数种子,方便重复性研究
train = sample(nrow(df), 0.7*nrow(df))  # 原数据的70%用来训练模型
df_train = df[train,]  # df_train为训练数据集
df_validate = df[-train,]  # df_validate为验证数据集

数据集中的变量

  • ID:样本ID号码
  • clumpThickness:肿块厚度
  • sizeUniformity:细胞大小的均匀性
  • shapeUniformity:细胞形状的均匀性
  • maginalAdhesion:边际附着力
  • singleEpithelialCellSize:单个上皮细胞大小
  • bareNuclei:裸核
  • blandChromatin:乏味染色体
  • normalNucleoli:正常核
  • mitosis:有丝分裂
  • class:类别。benign表示良性,malignant表示恶性。

数据说明

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
library(tidyverse)
df_train %>%
  mutate(value = 1) %>%
  group_by(class) %>%
  summarise(total = sum(value))->train_stat

df_validate %>%
  mutate(value = 1) %>%
  group_by(class) %>%
  summarise(total = sum(value))->validate_stat

第一个变量ID不纳入数据分析,最后一个变量(类别)即输出变量(编码为良性=2,恶性=4)。

对于每一个样本来说,另外九个变量是与判别恶性肿瘤相关的细胞特征,并且得到了记录。 这些细胞特征得分为1(最接近良性)至10(最接近病变)之间的整数。任一变量都不能单独作 为判别良性或恶性的标准,建模的目的是找到九个细胞特征的某种组合,从而实现对恶性肿瘤的 准确预测。Mangasarian和Wolberg在其1990年的文章中详细探讨了这个数据集。

下面给出R中数据准备流程。数据从UCI数据库中抽取,并随机分出训练集和验证集, 其中 训练集中包含489个样本单元 (占70%), 其中良性样本单元329个, 恶性样本单元160个; 验证集中包含210个样本单元 (占30%), 其中良性129个, 恶性81个。

逻辑回归(logistic regression)

逻辑回归(logistic regression)是广义线性模型的一种,可根据一组数值变量预测二元输出。 R中的基本函数 glm() 可用于拟合逻辑回归模型。 glm() 函数自动将预测变量中的分类变量编码为相应的虚拟变量。 威斯康星乳腺癌数据中的全部预测变量都是数值变量, 因此不必要对其编码。下面给出R中逻辑回归流程。

1
2
3
library(gmodels)
fit.logit = glm(class~., data=df_train, family=binomial())
summary(fit.logit)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
## 
## Call:
## glm(formula = class ~ ., family = binomial(), data = df_train)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.75813  -0.10602  -0.05679   0.01237   2.64317  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -10.42758    1.47602  -7.065 1.61e-12 ***
## clumpThickness             0.52434    0.15950   3.287  0.00101 ** 
## sizeUniformity            -0.04805    0.25706  -0.187  0.85171    
## shapeUniformity            0.42309    0.26775   1.580  0.11407    
## maginalAdhesion            0.29245    0.14690   1.991  0.04650 *  
## singleEpithelialCellSize   0.11053    0.17980   0.615  0.53871    
## bareNuclei                 0.33570    0.10715   3.133  0.00173 ** 
## blandChromatin             0.42353    0.20673   2.049  0.04049 *  
## normalNucleoli             0.28888    0.13995   2.064  0.03900 *  
## mitosis                    0.69057    0.39829   1.734  0.08295 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 612.063  on 482  degrees of freedom
## Residual deviance:  71.346  on 473  degrees of freedom
##   (6 observations deleted due to missingness)
## AIC: 91.346
## 
## Number of Fisher Scoring iterations: 8
1
2
3
4
5
prob = predict(fit.logit, df_validate, type="response")  # type="response"得出患肿瘤的概率
logit.pred = factor(prob > .5, levels=c(FALSE, TRUE),  # 概率大于0.5则为恶性,否则为良性
                     labels=c("benign", "malignant"))
logit.perf = CrossTable(df_validate$class, logit.pred,
                    dnn=c("Actual", "Predicted"))
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  200 
## 
##  
##              | Predicted 
##       Actual |    benign | malignant | Row Total | 
## -------------|-----------|-----------|-----------|
##       benign |       118 |         2 |       120 | 
##              |    27.419 |    42.885 |           | 
##              |     0.983 |     0.017 |     0.600 | 
##              |     0.967 |     0.026 |           | 
##              |     0.590 |     0.010 |           | 
## -------------|-----------|-----------|-----------|
##    malignant |         4 |        76 |        80 | 
##              |    41.128 |    64.328 |           | 
##              |     0.050 |     0.950 |     0.400 | 
##              |     0.033 |     0.974 |           | 
##              |     0.020 |     0.380 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       122 |        78 |       200 | 
##              |     0.610 |     0.390 |           | 
## -------------|-----------|-----------|-----------|
## 
## 
1
logit.perf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## $t
##            y
## x           benign malignant
##   benign       118         2
##   malignant      4        76
## 
## $prop.row
##            y
## x               benign  malignant
##   benign    0.98333333 0.01666667
##   malignant 0.05000000 0.95000000
## 
## $prop.col
##            y
## x               benign  malignant
##   benign    0.96721311 0.02564103
##   malignant 0.03278689 0.97435897
## 
## $prop.tbl
##            y
## x           benign malignant
##   benign      0.59      0.01
##   malignant   0.02      0.38

首先,以类别为响应变量,其余变量为预测变量。基于 df_train 数据框中的数据构造逻 辑回归模型。接着给出了模型中的系数。

接着,采用基于 df_train 建立的模型来对 df_validate 数据集中的样本单元分类。 predict() 函数默认输出肿瘤为恶性的对数概率,指定参数 type=“response” 即可得到预测肿 瘤为恶性的概率。样本单元中,概率大于0.5的被分为恶性肿瘤类,概率小于等于0.5的被分为 良性肿瘤类。

最后给出预测与实际情况对比的交叉表(即混淆矩阵,confusion matrix)。模型正确判别 了118个类别为良性的患者和76个类别为恶性的患者。另外,df_validate 数据集中有10个样本 单元因包含缺失数据而无法判别。

在验证集上,正确分类的模型(即准确率,accuracy)为(76+118)/200=97%,同时要注意的是,模型中有三个预测变量(sizeUniformity、shapeUniformity和singleEpithelialCellSize)的系数未通过显著性检验(即p值大于0.1)。从预测的角度来说,我们一般不会将这些变量纳入最终模型。当这类不包含相关信息的变量特别多时,可以直接将其认定为模型中的噪声。

在这种情况下,可用逐步逻辑回归生成一个包含更少解释变量的模型,其目的是通过增加或 移除变量来得到一个更小的AIC值。具体到这一案例,可通过: $$logit.fit.reduced <- step(fit.logit)$$ 来得到一个精简的模型。这样,上面提到的三个变量就从最终模型中移除,这种精简后的模型在 验证集上的误差相对全变量模型更小。

k最邻近分类算法 (KNN)

K最近邻(k-Nearest Neighbor,KNN)分类算法,是一个理论上比较不成熟的方法,也是最简单的机器学习算法之一。该方法的思路是:如果一个样本在特征空间中的k个最相似(即特征空间中最邻近)的样本中的大多数属于某一个类别,则该样本也属于这个类别。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 建立数据标准化函数
normalize = function(x){
  (x-min(x))/(max(x)-min(x))
}

df1=na.omit(df)  # 缺失值占比很少,此处直接删除,否则无法计算和邻居的距离
df1_n = as.data.frame(lapply(df1[1:9], normalize))  # 对除class列的数据进行标准化处理
set.seed(1234)
train = sample(nrow(df1_n), 0.7*nrow(df1_n))
df1_train = df1_n[train,]
df1_validate = df1_n[-train,]
df1_train_labels = df1[train, 10]
df1_validate_labels = df1[-train, 10]

library(class)
library(gmodels)
knn.pred = knn(train = df1_train,
               test = df1_validate,
               cl = df1_train_labels,
               k=22)    # 参数取22是因为训练的样本有478个,开根后是22
knn.perf = CrossTable(x=df1_validate_labels,y=knn.pred,
           dnn=c("Actual", "Predicted"))
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  205 
## 
##  
##              | Predicted 
##       Actual |    benign | malignant | Row Total | 
## -------------|-----------|-----------|-----------|
##       benign |       125 |         2 |       127 | 
##              |    22.030 |    40.695 |           | 
##              |     0.984 |     0.016 |     0.620 | 
##              |     0.940 |     0.028 |           | 
##              |     0.610 |     0.010 |           | 
## -------------|-----------|-----------|-----------|
##    malignant |         8 |        70 |        78 | 
##              |    35.870 |    66.259 |           | 
##              |     0.103 |     0.897 |     0.380 | 
##              |     0.060 |     0.972 |           | 
##              |     0.039 |     0.341 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       133 |        72 |       205 | 
##              |     0.649 |     0.351 |           | 
## -------------|-----------|-----------|-----------|
## 
## 
1
knn.perf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## $t
##            y
## x           benign malignant
##   benign       125         2
##   malignant      8        70
## 
## $prop.row
##            y
## x               benign  malignant
##   benign    0.98425197 0.01574803
##   malignant 0.10256410 0.89743590
## 
## $prop.col
##            y
## x               benign  malignant
##   benign    0.93984962 0.02777778
##   malignant 0.06015038 0.97222222
## 
## $prop.tbl
##            y
## x                benign   malignant
##   benign    0.609756098 0.009756098
##   malignant 0.039024390 0.341463415

经典决策树

经典决策树以一个二元输出变量(对应威斯康星州乳腺癌数据集中的良性/恶性)和一组预 测变量(对应九个细胞特征)为基础。具体算法如下。

1 选定一个最佳预测变量将全部样本单元分为两类,实现两类中的纯度最大化(即一类中 良性样本单元尽可能多,另一类中恶性样本单元尽可能多)。如果预测变量连续,则选定一个分 割点进行分类,使得两类纯度最大化;如果预测变量为分类变量(本例中未体现),则对各类别 进行合并再分类。

2 对每一个子类别继续执行步骤1。

3 重复步骤1~2,直到子类别中所含的样本单元数过少,或者没有分类法能将不纯度下 降到一个给定阈值以下。最终集中的子类别即终端节点(terminal node)。根据每一个终端节点中 样本单元的类别数众数来判别这一终端节点的所属类别。

4 对任一样本单元执行决策树,得到其终端节点,即可根据步骤3得到模型预测的所属类别。 不过,上述算法通常会得到一棵过大的树,从而出现过拟合现象。结果就是,对于训练集外 单元的分类性能较差。为解决这一问题,可采用10折交叉验证法选择预测误差最小的树。这一剪 枝后的树即可用于预测。

R中的rpart包支持rpart()函数构造决策树,prune()函数对决策树进行剪枝。下面给出 判别细胞为良性或恶性的决策树算法实现。

1
2
3
4
5
library(rpart)
set.seed(1234)
dtree = rpart(class ~ ., data=df_train, method="class",
              parms=list(split="information"))
dtree$cptable 
1
2
3
4
5
##         CP nsplit rel error  xerror       xstd
## 1 0.800000      0   1.00000 1.00000 0.06484605
## 2 0.046875      1   0.20000 0.30625 0.04150018
## 3 0.012500      3   0.10625 0.20625 0.03467089
## 4 0.010000      4   0.09375 0.18125 0.03264401

首先,rpart() 函数可用于生成决策树。 print(dtree) 和 summary(dtree) 可用于观测 所得模型,此时所得的树可能过大,需要剪枝。

rpart() 返回的 cptable 值中包括不同大小的树对应的预测误差,因此可用于辅助设定最终 的树的大小。其中,复杂度参数(cp)用于惩罚过大的树;树的大小即分支数(nsplit),有n 个分支的树将有n+1个终端节点; rel error 栏即训练集中各种树对应的误差;交叉验证误差 (xerror)即基于训练样本所得的10折交叉验证误差; xstd 栏为交叉验证误差的标准差。

1
plotcp(dtree)

1
dtree.pruned = prune(dtree, cp=.0125)

借助plotcp()函数可画出交叉验证误差与复杂度参数的关系图。对于所有 交叉验证误差在最小交叉验证误差一个标准差范围内的树,最小的树即最优的树。

复杂度参数与交叉验证误差。虚线是基于一个标准差准则得到的上限 (0.18+1×0.0326=0.21)。从图像来看,应选择虚线下最左侧 cp 值对应的树

本例中,最小的交叉验证误差为0.18,标准误差为0.0326,则最优的树为交叉验证误差在 0.18±0.0326(0.15和0.21)之间的树。由代码清单17-3的 cptable 表可知,四个终端节点(即三 次分割)的树满足要求(交叉验证误差为0.206 25);根据图也可以选得最优树,即三次分割 (四个节点)对应的树。

在完整树的基础上,prune()函数根据复杂度参数剪掉最不重要的枝,从而将树的大小控制 在理想范围内。cptable 中可以看到,三次分割对应的复杂度参数为0.0125, 从而 prune(dtree, cp=0.0125) 可得到一个理想大小的树。

1
2
3
library(rpart.plot)
prp(dtree.pruned, type = 2, extra = 104,
    fallen.leaves = TRUE, main="Decision Tree")

1
2
3
dtree.pred = predict(dtree.pruned, df_validate, type="class")
dtree.perf = CrossTable(df_validate$class, dtree.pred,
                    dnn=c("Actual", "Predicted"))
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  210 
## 
##  
##              | Predicted 
##       Actual |    benign | malignant | Row Total | 
## -------------|-----------|-----------|-----------|
##       benign |       122 |         7 |       129 | 
##              |    27.573 |    39.756 |           | 
##              |     0.946 |     0.054 |     0.614 | 
##              |     0.984 |     0.081 |           | 
##              |     0.581 |     0.033 |           | 
## -------------|-----------|-----------|-----------|
##    malignant |         2 |        79 |        81 | 
##              |    43.912 |    63.315 |           | 
##              |     0.025 |     0.975 |     0.386 | 
##              |     0.016 |     0.919 |           | 
##              |     0.010 |     0.376 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       124 |        86 |       210 | 
##              |     0.590 |     0.410 |           | 
## -------------|-----------|-----------|-----------|
## 
## 
1
dtree.perf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## $t
##            y
## x           benign malignant
##   benign       122         7
##   malignant      2        79
## 
## $prop.row
##            y
## x               benign  malignant
##   benign    0.94573643 0.05426357
##   malignant 0.02469136 0.97530864
## 
## $prop.col
##            y
## x               benign  malignant
##   benign    0.98387097 0.08139535
##   malignant 0.01612903 0.91860465
## 
## $prop.tbl
##            y
## x               benign  malignant
##   benign    0.58095238 0.03333333
##   malignant 0.00952381 0.37619048

rpart.plot 包中的 prp() 函数可用于画出最终的决策树。prp() 函数中有许多可供选择的参数 (详见?prp),如 type=2 可画出每个节点下分割的标签,extra=104可画出每一类的概率以及每个节点处的样本占比,fallen.leaves=TRUE可在图的底端显示终端节点。对观测点分类时,从树的顶端开始,若满足条件则从左枝往下,否则从右枝往下,重复这个过程直到碰到一个终端节点为止。该终端节点即为这一观测点的所属类别。

最后,predict() 函数用来对验证集中的观测点分类。整体来看,验证集中的准确率达到了96%。与逻辑回归不同的是,验证集中的210个样本单元都可由最终树来分类。值得注意的是,对于水平数很多或缺失值很多的预测变量,决策树可能会有偏。

条件推断树

在介绍随机森林之前,我们先介绍传统决策树的一种重要变体,即条件推断树(conditional inference tree)。条件推断树与传统决策树类似,但变量和分割的选取是基于显著性检验的,而不 是纯净度或同质性一类的度量。显著性检验是置换检验。

条件推断树的算法如下。

1 对输出变量与每个预测变量间的关系计算p值。

2 选取p值最小的变量。

3 在因变量与被选中的变量间尝试所有可能的二元分割(通过排列检验),并选取最显著的 分割。

4 将数据集分成两群,并对每个子群重复上述步骤。

5 重复直至所有分割都不显著或已到达最小节点为止。

条件推断树可由 party 包中的 ctree() 函数获得。

1
2
3
library(party)
fit.ctree = ctree(class~., data=df_train)
plot(fit.ctree, main="Conditional Inference Tree")

1
2
3
ctree.pred = predict(fit.ctree, df_validate, type="response")
ctree.perf = CrossTable(df_validate$class, ctree.pred,
                        dnn=c("Actual", "Predicted"))
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  210 
## 
##  
##              | Predicted 
##       Actual |    benign | malignant | Row Total | 
## -------------|-----------|-----------|-----------|
##       benign |       122 |         7 |       129 | 
##              |    26.624 |    39.153 |           | 
##              |     0.946 |     0.054 |     0.614 | 
##              |     0.976 |     0.082 |           | 
##              |     0.581 |     0.033 |           | 
## -------------|-----------|-----------|-----------|
##    malignant |         3 |        78 |        81 | 
##              |    42.401 |    62.354 |           | 
##              |     0.037 |     0.963 |     0.386 | 
##              |     0.024 |     0.918 |           | 
##              |     0.014 |     0.371 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       125 |        85 |       210 | 
##              |     0.595 |     0.405 |           | 
## -------------|-----------|-----------|-----------|
## 
## 
1
ctree.perf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## $t
##            y
## x           benign malignant
##   benign       122         7
##   malignant      3        78
## 
## $prop.row
##            y
## x               benign  malignant
##   benign    0.94573643 0.05426357
##   malignant 0.03703704 0.96296296
## 
## $prop.col
##            y
## x               benign  malignant
##   benign    0.97600000 0.08235294
##   malignant 0.02400000 0.91764706
## 
## $prop.tbl
##            y
## x               benign  malignant
##   benign    0.58095238 0.03333333
##   malignant 0.01428571 0.37142857

值得注意的是,对于条件推断树来说,剪枝不是必需的,其生成过程相对更自动化一些。另 外, party 包也提供了许多图像参数。上图展示了一棵条件推断树,每个节点中的阴影区域代 表这个节点对应的恶性肿瘤比例。

随机森林(random forest)

随机森林(random forest)是一种组成式的有监督学习方法。在随机森林中,我们同时 生成多个预测模型,并将模型的结果汇总以提升分类准确率。Leo Breiman和Adele Cutler在 https://mng.bz/7Nul1上有关于随机森林的详尽介绍。

随机森林的算法涉及对样本单元和变量进行抽样,从而生成大量决策树。对每个样本单元来 说,所有决策树依次对其进行分类。所有决策树预测类别中的众数类别即为随机森林所预测的这 一样本单元的类别。

假设训练集中共有N个样本单元,M个变量,则随机森林算法如下。

  • 从训练集中随机有放回地抽取N个样本单元,生成大量决策树。

  • 在每一个节点随机抽取m<M个变量,将其作为分割该节点的候选变量。每一个节点处的 变量数应一致。

  • 完整生成所有决策树,无需剪枝(最小节点为1)。

  • 终端节点的所属类别由节点对应的众数类别决定。

  • 对于新的观测点,用所有的树对其进行分类,其类别由多数决定原则生成。

生成树时没有用到的样本点所对应的类别可由生成的树估计,与其真实类别比较即可得到袋 外预测(out-of-bag,OOB)误差。无法获得验证集时,这是随机森林的一大优势。随机森林算 法可计算变量的相对重要程度,这将在下文中介绍。

randomForest 包中的 randomForest() 函数可用于生成随机森林。函数默认生成500棵树, 并且默认在每个节点处抽取 sqrt(M) 个变量,最小节点为1。

1
2
3
4
5
6
library(randomForest)
set.seed(1234)
fit.forest = randomForest(class~., data = df_train,
                          na.action = na.roughfix,
                          importance = TRUE)
print(fit.forest)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
## 
## 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: 3.48%
## Confusion matrix:
##           benign malignant class.error
## benign       320         9  0.02735562
## malignant      8       152  0.05000000

首先,randomForest() 函数从训练集中有放回地随机抽取489个观测点,在每棵树的每个 节点随机抽取3个变量,从而生成了500棵传统决策树。 na.action=na.roughfix 参数可将数 值变量中的缺失值替换成对应列的中位数,类别变量中的缺失值替换成对应列的众数类(若有多 个众数则随机选一个)。

1
importance(fit.forest, type=2)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
##                          MeanDecreaseGini
## clumpThickness                  11.847884
## sizeUniformity                  63.362023
## shapeUniformity                 44.921067
## maginalAdhesion                  5.442715
## singleEpithelialCellSize        12.725937
## bareNuclei                      35.489054
## blandChromatin                  18.472646
## normalNucleoli                  20.934481
## mitosis                          1.832045

随机森林可度量变量重要性,通过设置 importance = TRUE 参数得到,并通过 importance() 函数输出。由 type=2 参数得到的变量相对重要性就是分割该变量时节点不纯度(异质性)的 下降总量对所有树取平均。节点不纯度由Gini系数定义。本例中,sizeUniformity 是最重要的 变量,mitosis 是最不重要的变量。

1
2
3
forest.pred = predict(fit.forest, df_validate)
forest.perf = CrossTable(df_validate$class, forest.pred,
                    dnn=c("Actual", "Predicted"))
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  200 
## 
##  
##              | Predicted 
##       Actual |    benign | malignant | Row Total | 
## -------------|-----------|-----------|-----------|
##       benign |       117 |         3 |       120 | 
##              |    30.147 |    43.383 |           | 
##              |     0.975 |     0.025 |     0.600 | 
##              |     0.992 |     0.037 |           | 
##              |     0.585 |     0.015 |           | 
## -------------|-----------|-----------|-----------|
##    malignant |         1 |        79 |        80 | 
##              |    45.221 |    65.074 |           | 
##              |     0.012 |     0.988 |     0.400 | 
##              |     0.008 |     0.963 |           | 
##              |     0.005 |     0.395 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       118 |        82 |       200 | 
##              |     0.590 |     0.410 |           | 
## -------------|-----------|-----------|-----------|
## 
## 
1
forest.perf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## $t
##            y
## x           benign malignant
##   benign       117         3
##   malignant      1        79
## 
## $prop.row
##            y
## x           benign malignant
##   benign    0.9750    0.0250
##   malignant 0.0125    0.9875
## 
## $prop.col
##            y
## x                benign   malignant
##   benign    0.991525424 0.036585366
##   malignant 0.008474576 0.963414634
## 
## $prop.tbl
##            y
## x           benign malignant
##   benign     0.585     0.015
##   malignant  0.005     0.395

最后,再通过随机森林算法对验证集中的样本单元进行分类,并计算预测准确率。分类时 剔除验证集中有缺失值的单元。总体来看,对验证集的预测准确率高达98%。

randomForest 包根据传统决策树生成随机森林,而 party 包中的 cforest() 函数则可基于 条件推断树生成随机森林。当预测变量间高度相关时,基于条件推断树的随机森林可能效果更好。

相较于其他分类方法,随机森林的分类准确率通常更高。另外,随机森林算法可处理大规模 问题(即多样本单元、多变量),可处理训练集中有大量缺失值的数据,也可应对变量远多于样本 单元的数据。可计算袋外预测误差(OOB error)、度量变量重要性也是随机森林的两个明显优势。

随机森林的一个明显缺点是分类方法(此例中相当于500棵决策树)较难理解和表达。另外, 我们需要存储整个随机森林以对新样本单元分类。

支持向量机(SVM)

支持向量机(SVM)是一类可用于分类和回归的有监督机器学习模型。其流行归功于两个方 面:一方面,他们可输出较准确的预测结果;另一方面,模型基于较优雅的数学理论。

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

对于一个N维空间(即N个变量)来说,最优超平面(即线性决策面,linear decision surface) 为N–1维。当变量数为2时,曲面是一条直线;当变量数为3时,曲面是一个平面;当变量数为10 时,曲面就是一个九维的超平面。当然,这并不是太好想象。

SVM可以通过R中 kernlab 包的 ksvm() 函数和 e1071 包中的 svm() 函数实现。ksvm()功能 更强大,但svm()相对更简单。

1
2
3
4
library(e1071)
set.seed(1234)
fit.svm = svm(class~., data=df_train)
fit.svm
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
## 
## Call:
## svm(formula = class ~ ., data = df_train)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.1111111 
## 
## Number of Support Vectors:  76
1
2
3
svm.pred = predict(fit.svm, na.omit(df_validate))
svm.perf = CrossTable(na.omit(df_validate)$class,
                 svm.pred, dnn=c("Actual", "Predicted"))
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  200 
## 
##  
##              | Predicted 
##       Actual |    benign | malignant | Row Total | 
## -------------|-----------|-----------|-----------|
##       benign |       116 |         4 |       120 | 
##              |    27.859 |    40.929 |           | 
##              |     0.967 |     0.033 |     0.600 | 
##              |     0.975 |     0.049 |           | 
##              |     0.580 |     0.020 |           | 
## -------------|-----------|-----------|-----------|
##    malignant |         3 |        77 |        80 | 
##              |    41.789 |    61.394 |           | 
##              |     0.037 |     0.963 |     0.400 | 
##              |     0.025 |     0.951 |           | 
##              |     0.015 |     0.385 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       119 |        81 |       200 | 
##              |     0.595 |     0.405 |           | 
## -------------|-----------|-----------|-----------|
## 
## 
1
svm.perf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## $t
##            y
## x           benign malignant
##   benign       116         4
##   malignant      3        77
## 
## $prop.row
##            y
## x               benign  malignant
##   benign    0.96666667 0.03333333
##   malignant 0.03750000 0.96250000
## 
## $prop.col
##            y
## x               benign  malignant
##   benign    0.97478992 0.04938272
##   malignant 0.02521008 0.95061728
## 
## $prop.tbl
##            y
## x           benign malignant
##   benign     0.580     0.020
##   malignant  0.015     0.385

由于方差较大的预测变量通常对SVM的生成影响更大,svm() 函数默认在生成模型前对每个 变量标准化,使其均值为0、标准差为1。从结果来看,SVM的预测准确率不错,但不如随机森林方法, 与随机森林算法不同的是,SVM在预测新样本单元时不允许有缺失值出现。

svm() 函数默认通过径向基函数(Radial Basis Function,RBF)将样本单元投射到高维空间。 一般来说RBF核是一个比较好的选择,因为它是一种非线性投影,可以应对类别标签与预测变量 间的非线性关系。

在用带RBF核的SVM拟合样本时,两个参数可能影响最终结果:gamma和成本(cost)。gamma是核函数的参数,控制分割超平面的形状。gamma越大,通常导致支持向量越多。我们也可将gamma看作控制训练样本“到达范围”的参数,即gamma越大意味着训练样本到达范围越广,而越小则意味着到达范围越窄。gamma必须大于0。

成本参数代表犯错的成本。一个较大的成本意味着模型对误差的惩罚更大,从而将生成一个 更复杂的分类边界,对应的训练集中的误差也会更小,但也意味着可能存在过拟合问题,即对新 样本单元的预测误差可能很大。较小的成本意味着分类边界更平滑,但可能会导致欠拟合。与 gamma一样,成本参数也恒为正。

svm() 函数默认设置gamma为预测变量个数的倒数,成本参数为1。不过gamma与成本参数的 不同组合可能生成更有效的模型。在建模时,我们可以尝试变动参数值建立不同的模型,但利用 格点搜索法可能更有效。可以通过 tune.svm() 对每个参数设置一个候选范围, tune.svm() 函 数对每一个参数组合生成一个SVM模型,并输出在每一个参数组合上的表现。

  • 带RBF核的SVM模型
1
2
3
4
5
set.seed(1234)
tuned = tune.svm(class~., data=df_train,
                  gamma=10^(-3:1),
                  cost=10^(-3:1))
print(tuned)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  gamma cost
##   0.01    1
## 
## - best performance: 0.02904092

首先,对不同的gamma和成本拟合一个带RBF核的SVM模型。我们一共将尝试八个不同的gamma以成本参数。总体来说,我们共拟合了多个模型,并比较了其结果。训练集中10折交叉验证误差最小的模型所对应的参数为 gamm=0.1,成本参数为1。

1
2
3
4
fit.svm = svm(class~., data=df_train, gamma=.01, cost=1)
svm.pred = predict(fit.svm, na.omit(df_validate))
svm.perf = CrossTable(na.omit(df_validate)$class, svm.pred,
                 dnn=c("Actual", "Predicted"))
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  200 
## 
##  
##              | Predicted 
##       Actual |    benign | malignant | Row Total | 
## -------------|-----------|-----------|-----------|
##       benign |       117 |         3 |       120 | 
##              |    28.125 |    42.188 |           | 
##              |     0.975 |     0.025 |     0.600 | 
##              |     0.975 |     0.037 |           | 
##              |     0.585 |     0.015 |           | 
## -------------|-----------|-----------|-----------|
##    malignant |         3 |        77 |        80 | 
##              |    42.188 |    63.281 |           | 
##              |     0.037 |     0.963 |     0.400 | 
##              |     0.025 |     0.963 |           | 
##              |     0.015 |     0.385 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |       120 |        80 |       200 | 
##              |     0.600 |     0.400 |           | 
## -------------|-----------|-----------|-----------|
## 
## 
1
svm.perf
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## $t
##            y
## x           benign malignant
##   benign       117         3
##   malignant      3        77
## 
## $prop.row
##            y
## x           benign malignant
##   benign    0.9750    0.0250
##   malignant 0.0375    0.9625
## 
## $prop.col
##            y
## x           benign malignant
##   benign    0.9750    0.0375
##   malignant 0.0250    0.9625
## 
## $prop.tbl
##            y
## x           benign malignant
##   benign     0.585     0.015
##   malignant  0.015     0.385

基于这一参数值组合,我们对全部训练样本拟合出新的SVM模型,然后用这一模型对验证 集中的样本单元进行预测,并给出错分个数。在本例中,调和后的模型微减少了错分个数 (从7减少到6)。一般来说,为SVM模型选取调和参数通常可以得到更好的结果。

如前所述,由于SVM适用面比较广,它目前是很流行的一种模型。SVM也可以应用于变量 数远多于样本单元数的问题,而这类问题在生物医药行业很常见,因为在DNA微序列的基因 表示中,变量数通常比可用样本量的个数高1~2个量级。

与随机森林类似,SVM的一大缺点是分类准则比较难以理解和表述。SVM从本质上来说是 一个黑盒子。另外,SVM在对大量样本建模时不如随机森林,但只要建立了一个成功的模型, 在对新样本分类时就没有问题了。

预测效果评价

预测准确性度量

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
library(knitr)
library(kableExtra)
stat = c("敏感度","特异性","正例命中率","正例命中率","准确率")
explain = c(
  "正类的样本单元被成功预测的概率,也叫正例覆盖率(true positive)或召回率(recall)",
  "负类的样本单元被成功预测的概率,也叫负例覆盖率(true negative)",
  "被预测为正类的样本单元中,预测正确的样本单元占比,也叫精确度(precision)",
  "被预测为负类的样本单元中,预测正确的样本单元占比",
  "被正确分类的样本单元所占比重,也叫ACC")
data.frame(stat, explain) %>%
  kable() %>%
  kable_styling("striped", full_width = F) 
stat explain
敏感度 正类的样本单元被成功预测的概率,也叫正例覆盖率(true positive)或召回率(recall)
特异性 负类的样本单元被成功预测的概率,也叫负例覆盖率(true negative)
正例命中率 被预测为正类的样本单元中,预测正确的样本单元占比,也叫精确度(precision)
正例命中率 被预测为负类的样本单元中,预测正确的样本单元占比
准确率 被正确分类的样本单元所占比重,也叫ACC

编写函数计算统计量

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
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)
}

给出计算结果

  • 逻辑回归模型
1
performance(logit.perf$`t`)
1
2
3
4
5
## Sensitivity = 0.95
## Specificity = 0.98
## Positive Predictive Value = 0.97
## Negative Predictive Value = 0.97
## Accuracy = 0.97
  • k近邻模型
1
performance(knn.perf$`t`)
1
2
3
4
5
## Sensitivity = 0.9
## Specificity = 0.98
## Positive Predictive Value = 0.97
## Negative Predictive Value = 0.94
## Accuracy = 0.95
  • 经典决策树
1
performance(dtree.perf$`t`)
1
2
3
4
5
## Sensitivity = 0.98
## Specificity = 0.95
## Positive Predictive Value = 0.92
## Negative Predictive Value = 0.98
## Accuracy = 0.96
  • 条件推断数
1
performance(ctree.perf$`t`)
1
2
3
4
5
## Sensitivity = 0.96
## Specificity = 0.95
## Positive Predictive Value = 0.92
## Negative Predictive Value = 0.98
## Accuracy = 0.95
  • 随机森林
1
performance(forest.perf$`t`)
1
2
3
4
5
## Sensitivity = 0.99
## Specificity = 0.98
## Positive Predictive Value = 0.96
## Negative Predictive Value = 0.99
## Accuracy = 0.98
  • 支持向量机
1
performance(svm.perf$`t`)
1
2
3
4
5
## Sensitivity = 0.96
## Specificity = 0.98
## Positive Predictive Value = 0.96
## Negative Predictive Value = 0.98
## Accuracy = 0.97

在这个案例中,这些分类器(逻辑回归、K近邻、传统决策树、条件推断树、随机森林和支持向量机) 都表现得相当不错。不过在现实中并不总是这样。

在这个案例中,随机森林的表现相对更好。不过各个分类器的差距较小,因此随机森林的优 势可能具有一定的偶然性。随机森林成功鉴别了99%的恶性样本和98%的良性样本,总体来说预 测准确率高达99%。96%被判为恶性组织的样本单元确实是恶性的(即4%正例错误率),99%被 判为良性组织的样本单元确实是良性的(即1%负例错误率)。从癌症诊断的角度来说,特异性(即 成功鉴别恶性样本的概率)这一指标格外重要。