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%负例错误率）。从癌症诊断的角度来说，特异性（即 成功鉴别恶性样本的概率）这一指标格外重要。

  1. 随机森林官网：https://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm [return]