StatQuest学习笔记20——随机森林

前言——主要内容

这篇笔记是StatQuest系列视频的第53-55节，其中第53节讲的是随机森林，第54节讲的是缺失值的处理，第55节讲的是R与随机森林。

决策树的局限

随机森林（Random Forests）来源于决策树，因此如果不了解决策树的话，可以看上一篇笔记，决策树如下所示：

决策树很容易构建，也很好用，也很好题解。但事实上，决策树并非完美无缺，引用《The Elements of Statistical Learning》（又名《The Bible of Machine Learning》），中文译名为《统计学习基础：数据挖掘、推理与预测》的话讲，决策树的不精确准确（inaccuracy）决定了它无法作为一个理想的预测学习工具，如下所示：

换句话说，就是决策树的构建要利用大量的数据，但是当用这个构建好的决策树来对一批新的数据进行分类时，决策树并不能灵活地处理这些新数据，如下所示：

随机森林的优势

随机森林综合了决策树简洁的特征，同时又具备灵活性，因此随机森林在精确性方面也得到了极大的提高，如下所示：

构建随机森林

接着，我们来看一下随机森林的构建过程。

第一步：构建自举数据集

构建随机森林的第一步就是从原始数据中随机挑选数据，构建“自举”数据集（自举：bootstrap），如下所示：

下面是原始的数据集，这个数据集还是我们上篇笔记中的那个数据集，里面有4个变量，分别为胸痛（Chest Pain），良好血液循环（Good Blood Circ.），动脉阻塞（Blocked Arteries），体重（Weight），虽然这个数据集非常小，但是不影响使用，如下所示：

为了构建一个与原始数据集同样大小的自举数据集（bootstrapped dataset），我们仅需要随机地从原始数据中挑选样本即可，这个随机挑选的过程是有放回地挑选（这就是bootstrap），也就是我从原始数据集中随机挑一个样本，然后放回去，再随机挑一个样本，这样的话，相同有样本有可能被挑中好几次，如下所示：

我们先随机挑一个样本，这次挑中了第2个样本，如下所示：

把这个样本放到我们的自举数据集中，就成了自举数据集中的第1个样本，如下所示：

再随机挑一个样本，这次挑中了原始数据集中的第1个样本，把它放到我们的自举数据集中，就是自举数据集中的第2个样本，如下所示：

再从原始数据中随机挑一个，放到自举数据集中，如下所示：

再随机挑一个，放到自举数据集中，这里需要注意的是，我们挑了两次原始数据集中第4个样本，这样在自举数据集中就有了2个重复的样本，这个没关系，如下所示：

此时，我们就构建好了一个自举数据集，这个自举数据集的大小与原始数据集的大小相同，也就是说有同样数目的样本，如下所示：

第二步：构建树

构建随机森林的第二步就是利用我们构建的自举数据集来构建一个决策树，并且一次只能使用一个变量（或列）的随机亚集（subset）进行构建，在这个案例中，我们每一步只会使用2个变量（列），至于为什么使用2个变量，后文会讲解，如下所示：

因此，我们不会使用所有的4个变量来区分根节点，如下所示：

而是随机选择2个变量来区分根节点，如下所示：

在这个案例中，我们随机使用良好血液循环（Good Blood Circulation）和动脉阻塞（Blocked Arteries）这两个变量来区分根节点，如下所示：

就像我们这个案例的目的一样，我们假设良好血液循环（Good Blood Circulation）能够很好地区分样本，我们就用它作为根节点，如下所示：

由于我们使用过了良好血液循环（Good Blood Circulation）这个变量，因此这里把它变灰，以便于我们研究剩余的变量，如下所示：

现在我们看一下中间节点如何区分，如下所示：

与根节点一样，我们并不使用剩余的所有3个变量，而是随机选择2个变量用于区分，如下所示：

接着，我们就像以前构建决策树那样，构建树，不过我们一次只使用一个变量的随机亚集，如下所示：

最终的结果如下所示：

现在回顾一下我们构建这个树的过程：第一步，构建自举数据集；第二步，每一步随机选择一个变量亚集来构建节点，如下所示：

第三步：循环第1到到第2步

此时，我们就构建好了一个树，现在我们再回到第一步，重复操作，再生成一个自举数据集，还是一次使用变量的一个随机亚集来构建一个树（就是下图中第2个树），如下所示：

不断地重复这个过程，重复100次，这样我们会构建一系列的树，下图只是显示了其中的6种树，如下所示：

通过利用一个自举样本，一次只利用变量的一个随机亚集，就会产生大量不同的树，这些树就构成了随机森林（一次生成一棵树，多次运算，就是多棵树，也就是森林了），它远比单一的决策树有效，如下所示：

随机森林的使用

此时，我们已经学会了构建随机森林的基本思想，那么随机森林有什么用处呢，如下所示：

我们看一个案例，例如，我们遇到了一个新病人，他的各种参数如下所示：

现在我们想知道，他是否得有心脏病，如下所示：

那么，我们把他的数据放到我们生成随机森林中，先放到第1棵树上，如下所示：

运算结果显示Yes，此时记录一下这个数据，如下所示：

现在把这个人的参数放到随机森林中的第2棵树上运行，如下所示：

计算结果是Yes，此时记录下这个数据，如下所示：

我们把这个人的数据再放到随机森林中的其他的树中运算，如下所示：

总之，就是不断地放到不同的树上运算，如下所示：

当运算结果后，我们记录下所有树的运算结果，我们会看到Yes或No哪个得票最多，如下所示：

在这个案例中，明显Yes得票最多，因此我们就能下结论，这个人得了心脏病，如下所示：

随机森林相关术语

自举数据并对它们进行汇总做出的决策称为袋装法（Bagging），如下所示：

随机森林的评估

我们已经了解了随机森林构建的基本思想，也知道了如何使用它，此时可能还有一个问题，这个模型的好坏程度到底如何，如下所示：

我们回到前面构建自举数据集的内容，我们在前文中说过，在一个自举数据集中，我们允许两条相同的记录（entry，记录可以理解为一个样本，也就是一行数据）存在，如下所示：

最终在自举数据集中，没有包含原始数据集中的第3条目录，如下所示：

通常来说，大概会有三分之一的原始数据集中的记录不会包含在自举数据集中，如下所示：

下图的右侧就是那条未包含在自举数据集中的记录，如下所示：

如果原始数据集非常大的话，那么就有不止一条记录不包含在自举数据集中，如下所示：

这些不包含在自举数据集上的记录被称为“出袋数据集”（out-of-bag dataset，未找到相应的中文译名，此处暂译为“出袋数据集”），这样命名主要是因为，如果把自举数据集看成一个袋子，这些数据并不在这个袋子中，如下所示：

由于这些出袋数据集并没有用于构建这些树，如下所示：

因此，我们可以把这个数据集当作是一个新的样本，把它放到随机森林中计算，看能不能得到想要的结果（没得心脏病），如下所示：

在这个案例中，这个样本最终的计算结果是No，也就是没得心脏病，如下所示：

然后我们把这些出袋数据集放到其他树中计算，如下所示：

有一个树计算的结果是Yes，这个明显是错误的，如下所示：

其他的树计算结果是No，如下所示：

对所有树的计算结果进行统计，我们发现，No的得票最多，那么我们就能下结论，这个样本的结果是没得心脏病，与样本本身的结果一致，如下所示：

我们对其他的出袋数据集继续进行计算（出袋数据据不止一个样本，有很多），这个样本的结果是Yes，与样本本身的结果一致，如下所示：

再计算另外的样本，样本本身的结果是No，但计算出来的结果却是Yes，并不一致，如下所示：

经过不断地计算，我们会得到很多计算结果，最终我们再计算一下，随机森林计算的样本结果与样本本身的结果一致的比例，我们就能得到随机森林的精确程度，其中随机森林计算的结果与样本本身的结果不一致的比例，我们称为出袋错误（Out-Of-Bag Error）（未找到相应的中文译名，此处暂译为“出袋错误”），如下所示：

至此，我们了解了这些内容：第一，构建随机森林；第二，使用随机森林；第三，评估随机森林的精确性。

既然我们知道了如何评估随机森林的精确性，那么此时我们再回到第一部分，了解一下随机森林是的构建原理，如下所示：

在构建随机森林时，我们最初在每一步中，使用的是随机选择的2个变量，如下所示：

现在我们就要比较一下，使用随机选择的2个变量与3个变量导致的出袋错误的差异，我们测试一些不同的参数，来研究一下如何选择一个最为精确的随机森林，如下所示：

换句话讲，我们做上述事件的顺序是：

第一，构建随机森林；

第二，评估随机森林的精确性。

接着，改变每一个步骤中的变量数目，不断地重复这个过程，得到一个最精确的随机森林，如下所示：

通常来说，我们最初使用的变量数目为总变量数目的平方根，例如在前面的案例中，总变量数目是4，它的平方根是2，那么我们一开始就先使用2个变量作为每次运算时变量的数目，然后在这个数目上下进行调整，如下所示：

在构建随机森林的过程中我们可能会遇到缺失值，缺失值有两种，我们先看第一种缺失值的处理。

第一种缺失值的处理

我们看一下我们的原始数据，如下所示：

原始数据一共是4个患者，其中第4名患者有2个数据缺失，用问号表示，如下所示：

随机森林中有2类缺失值，分别为：第一，用于构建随机森林的原始数据集中的数据缺失，就是上面的情况；第二，用于归类的，一个新的样本中的数据缺失情况，如下所示：

我们先讲第一种数据的缺失情况，如下所示：

我们要从上面的原始数据集中构建随机森林，但是不清楚第4个患者的动脉阻塞与体重数据，如下所示：

处理这种缺失数据的常规方式是观察其他数据，先做出一个可能不太好的初始猜测，然后逐步改进这个猜测，直到这个猜测变得比较好为止（希望如此），如下所示：

现在我们就做这样的猜测，我们先看血管阻塞这个变量，然后找到这个变量出现最多的数据就可做出猜测，如下所示：

我们发现，在这个变量中，No出现的频率最高（2个No，1个Yes），如下所示：

因此，我们就可以猜测，在这个变量中，第4个患者的值可能是No，如下所示：

猜测完血管阻塞这个变量中的缺失值后，我们再猜测一下体重这个变量在第4个患者中的缺失值，体重是一种连续型变量（数值），因此我们可以看一下这个变量中其他数值的中位数即可（在决策树中，这种缺失值的处理，要么使用中位数，要么使用中位数），如下所示：

在其他的患者中，体重这个变量的中位数是180，因此我们就可以把第4个患者的体重缺失值写为180，如下所示：

经过这样的处理，原始数据中缺失的数据就填充好了，如下所示：

填充好缺失值后，我们还要评估这些数据猜测得到底合理不合理，现在就需要看一下哪些样本与含有缺失值的样本比较类似，接着我们就会讲一下这种相似性如何计算，如下所示：

第一步：构建随机森林，如下所示：

第二步：将所有的原始数据放到随机森林中运行，如下所示：

我们先在随机森林中的第1棵树上运行，如下所示：

把所有患者的数据依次放进去，如下所示：

这4个患者的数据经过计算后，我们会发现，患者3和患者4最终位于同一片叶子节点上，如下所示：

这就说明，第3个患者与第4个患者相似，如下所示：

相似性矩阵

现在我们用一个相似性矩阵（Proximity Matrix）来研究相似的患者数据，如下所示：

相似性矩阵的每行代表一个样本，如下所示：

相似性矩阵的每列也代表一个样本，如下所示：

因为第3个样本与第4个样本最终位于同一个叶子节点上，如下所示：

我们把患者3和患者4交叉的地方填上1，如下所示：

此外，还有一个交叉点，也填上1，如下所示：

由于在随机森林中第1棵树上运行后，除了第3个患者和第4个患者的最终结果相同外，并没有其他的患者结果相同，因此我们在第1棵树上运行后，得到的相似性矩阵就如下所示：

现在我们把所有的数据放到随机森林的第2棵树上运行，如下所示：

运行过程如下所示：

最终我们发现，患者2，患者3和患者4位于同一个叶子节点上，如下所示：

我们再看一下，经过第1棵树运行后的相似性矩阵如下所示：

经过第2棵树的运算，我们把有相同结果的患者交叉处填上1，如果原来的交叉位置有1，就加1，如下所示：

从中我们可以看到，患者3和患者4还是位于同一个叶子节点上，因此就在原有的基础上再加1，如下所示：

患者2和患者3位于同一个节点，患者2和患者4位于同一个节点，如下所示：

现在我们把所有的数据放到第3棵树上运行，如下所示：

得到结果后，我们再更新一下相似性矩阵，如下所示：

其中我们发现，患者3和患者4还是位于同一个节点上，再加1，如下所示：

最终，我们运行完所有的树，得到一个最终的相似性矩阵，如下所示：

然后，我们把这个相似性矩阵中的所有值除以所有树的数目，在这个案例中，所有的树的数目是10，得到结果如下所示：

现在我们使用这个近似矩阵来评估一下我们对缺失值猜测的效果如何，如下所示：

对于动脉阻塞这个变量来说，我们使用这个近似值来计算Yes与No的加权频率，对于体重缺失值的处理，也是如此，如下所示：

其中动脉阻塞这个变量中，除了缺失值外，在剩下的数据中，Yes的频率是1/3，No的频率是2/3，如下所示：

那么对Yes的加权频率计算如下所示：

患者2中Yes对应的近似值为0.1，如下所示：

然后除以患者4中所有近似值的和，如下所示：

Yes的权重是0.1，如下所示：

那么Yes的加权频率就是0.03，如下所示：

再来计算No的权加频率，如下所示：

它的公式为：

No的权重是0.9，如下所示：

最终计算结果是0.6，我们发现，No的加权频率更高，也就是说，我们可以在这个缺失值处填上No，如下所示：

对于体重缺失值的处理，我们也是采用近似值来计算它的加权平均数，如下所示：

加权平均数的的公式如下所示：

第一个样本的加权平均体重是0.1，用125乘以它即可，如下所示：

接着计算患者2的权加平均数，如下所示：

最终计算出加权平均数，如下所示：

此时，把这个结果（也就是198.5）填到缺失值处，如下所示：

现在我们再重新审视一下我们的猜测，重新再运行一次数据，然后构建一个随机森林，在所有的树中计算数据，重新计算近似值，再次计算缺失值，这个过程要重复6或7次，直接覆盖所有的缺失值（直到我们再次计算时，这个数据不再变化为止），如下所示：

相似矩阵的拓展

这里我们先看一下相似矩阵的一些很有意思的东西，如下所示：

我们再看一下在除以10（随机森林的所有树的数目）之间前的相似矩阵，如下所示：

如果样本3和样本4最终都落在了所有的10个树的同一个叶子节点上，那么上述的数字8就是10，如下所示：

我们再除以10（所有树的数目），我们就看到在这个相似矩阵中最大的数值就是1，如下所示：

在这个近似矩阵中，1意味着这两个样本是尽可能接近的，如下所示：

此时，我们用1减去这些近似值，就得到了某个距离，那么样本3和样本4的近似值是1，1减去1，就是0，也就是说这两个样本之间没有距离，如下所示：

样本1和样本4之间的的估计值是0.1，用1减去这个数值，得到0.9，那么这就表示，这两个样本的距离很大，也就表示它们之间离得远，如下所示：

我们把1减去这个近似值矩阵得到的新矩阵称为距离矩阵（distance matrix），画出热图就是下面的这个样子：

我们也可以绘制出MDS图，如下所示：

我们对相似矩阵的这些拓展计算很有意思，因为，无论你的数据是什么形式（有序数据，多个选择，连续型数据），如果我们使用这些数据构建一个树，我们就能画出热图或者是MDS图，用于研究不同样本的之间的相似性，如下所示：

第二种缺失值的处理

此时，再回到缺失值的处理上来，前面我们只是介绍了一种缺失值的处理，也就是用于构建随机森林的原始数据中的缺失值处理，现在我们介绍一下另外一种缺失值的处理，也就是我们想用于分类的新的样本中的缺失值处理，如下所示：

试想，我们已经通过现有的数据构建了一个随机森林，此时我们要将下面的的这个患者进行分类，不过他的一个动脉阻塞参数是缺失的，如下所示：

但我们还是想要知道这个患者是否患有心脏病，如下所示：

但我们不知道动脉阻塞这个参数，因此我们需要对这个参数进行一个猜测，其过程就是把这个患者的其他参数放到随机森林中运行，如下所示：

第一步：我们先猜测这个患者得了心脏病与不得心脏病这两种情况，做出两份数据，如下所示：

然后我们使用先前提到的迭代方法来对这动脉阻塞这个参数进行猜测，如下所示：

猜测结果一个是Yes，一个是No，如下所示：

然后把这两个数据放到随机森林中运算，如下所示：

然后我们观察这两个数据在随机森林中的运行结果，看一下哪个结果数目最多，如下所示：

在第一个数据中，所有的3个树都标记了其为Yes，如下所示：

在第二个数据中，只有一个树正确地标记了No，如下所示：

因此第一个数据是正确的数据，因为所有的树中它被正确计算的次数最多，如下所示：

R语言与随机森林

加载包

在使用R语言构建随机森林的时候，需要使用如下三个包：

1
2
3

library(ggplot2)
library(cowplot)
library(randomForest)

数据集的处理

原始数据使用UCI机器学习的数据，这里的数据很多，我们只用其中的心脏病方面的一个数据集，如下所示：

代码如下所示：

1
2
3

url <- "http://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.cleveland.data"
data <- read.csv(url, header=FALSE)
head(data)

数据如下所示：

> head(data)
  V1 V2 V3  V4  V5 V6 V7  V8 V9 V10 V11 V12 V13 V14
1 63  1  1 145 233  1  2 150  0 2.3   3 0.0 6.0   0
2 67  1  4 160 286  0  2 108  1 1.5   2 3.0 3.0   2
3 67  1  4 120 229  0  2 129  1 2.6   2 2.0 7.0   1
4 37  1  3 130 250  0  0 187  0 3.5   3 0.0 3.0   0
5 41  0  2 130 204  0  2 172  0 1.4   1 0.0 3.0   0
6 56  1  2 120 236  0  0 178  0 0.8   1 0.0 3.0   0

从数据中我们可以发现，第一行对于每个变量的命名并不是特别的规范，因此，根据网站提供的信息，我们把这些变量名改一下，如下所示：

colnames(data) <- c(
  "age",
  "sex",
  "cp",
  "trestbps",
  "chol",
  "fbs",
  "restecg",
  "thalach",
  "exang",
  "oldpeak",
  "slope",
  "ca",
  "thal",
  "hd"
)
head(data)

更改后的结果如下所示：

> head(data)
  age sex cp trestbps chol fbs restecg thalach exang oldpeak slope  ca thal hd
1  63   1  1      145  233   1       2     150     0     2.3     3 0.0  6.0  0
2  67   1  4      160  286   0       2     108     1     1.5     2 3.0  3.0  2
3  67   1  4      120  229   0       2     129     1     2.6     2 2.0  7.0  1
4  37   1  3      130  250   0       0     187     0     3.5     3 0.0  3.0  0
5  41   0  2      130  204   0       2     172     0     1.4     1 0.0  3.0  0
6  56   1  2      120  236   0       0     178     0     0.8     1 0.0  3.0  0

查看一下数据结构，如下所示：

> str(data)
'data.frame':	303 obs. of  14 variables:
 $ age     : num  63 67 67 37 41 56 62 57 63 53 ...
 $ sex     : num  1 1 1 1 0 1 0 0 1 1 ...
 $ cp      : num  1 4 4 3 2 2 4 4 4 4 ...
 $ trestbps: num  145 160 120 130 130 120 140 120 130 140 ...
 $ chol    : num  233 286 229 250 204 236 268 354 254 203 ...
 $ fbs     : num  1 0 0 0 0 0 0 0 0 1 ...
 $ restecg : num  2 2 2 0 2 0 2 0 2 2 ...
 $ thalach : num  150 108 129 187 172 178 160 163 147 155 ...
 $ exang   : num  0 1 1 0 0 0 0 1 0 1 ...
 $ oldpeak : num  2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
 $ slope   : num  3 2 2 3 1 1 3 1 2 3 ...
 $ ca      : Factor w/ 5 levels "?","0.0","1.0",..: 2 5 4 2 2 2 4 2 3 2 ...
 $ thal    : Factor w/ 4 levels "?","3.0","6.0",..: 3 2 4 2 2 2 2 2 4 4 ...
 $ hd      : int  0 2 1 0 0 0 3 0 2 1 ...

现在解释一下各个变量：

sex是性别，其中0表示女性，1表示男性，在分析时，需要把sex这个变量根据1或0更改为M或F。

cp是chest pain的缩写，其数值是1-3，这个范围表示不同类型的胸痛，4表示没有胸痛。

ca与thal是因子数据类型，我们看到，其中有的数值是?，这个是不符合规范的，因此后面我们会把它改为NA。

此外，还要把其它的变量更改为因子类型，如下所示：

data[data == "?"] <- NA
data[data$sex == 0,]$sex <- "F"
data[data$sex == 1,]$sex <- "M"
data$sex <- as.factor(data$sex)
data$cp <- as.factor(data$cp)
data$fbs <- as.factor(data$fbs)
data$restecg <- as.factor(data$restecg)
data$exang <- as.factor(data$exang)
data$slope <- as.factor(data$slope)
data$ca <- as.integer(data$ca)
data$ca <- as.factor(data$ca)
data$thal <- as.factor(data$thal)
data$hd <- ifelse(test=data$hd == 0, yes = "Healthy", no="Unhealthy")
data$hd <- as.factor(data$hd)

查看一下数据类型，如下所示：

> str(data)
'data.frame':	303 obs. of  14 variables:
 $ age     : num  63 67 67 37 41 56 62 57 63 53 ...
 $ sex     : Factor w/ 2 levels "F","M": 2 2 2 2 1 2 1 1 2 2 ...
 $ cp      : Factor w/ 4 levels "1","2","3","4": 1 4 4 3 2 2 4 4 4 4 ...
 $ trestbps: num  145 160 120 130 130 120 140 120 130 140 ...
 $ chol    : num  233 286 229 250 204 236 268 354 254 203 ...
 $ fbs     : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 1 1 1 2 ...
 $ restecg : Factor w/ 3 levels "0","1","2": 3 3 3 1 3 1 3 1 3 3 ...
 $ thalach : num  150 108 129 187 172 178 160 163 147 155 ...
 $ exang   : Factor w/ 2 levels "0","1": 1 2 2 1 1 1 1 2 1 2 ...
 $ oldpeak : num  2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
 $ slope   : Factor w/ 3 levels "1","2","3": 3 2 2 3 1 1 3 1 2 3 ...
 $ ca      : Factor w/ 4 levels "2","3","4","5": 1 4 3 1 1 1 3 1 2 1 ...
 $ thal    : Factor w/ 4 levels "?","3.0","6.0",..: 3 2 4 2 2 2 2 2 4 4 ...
 $ hd      : Factor w/ 2 levels "Healthy","Unhealthy": 1 2 2 1 1 1 2 1 2 2 ...

补充缺失值

由于我们要生成随机森林，此时我们设一个种子，用于生成随机数，这样我们可以重复我们的计算结果，在构建随机森林之间，由于我们的数据集中的缺失值，此时要对缺失值进行补充，补充缺失值使用的是rfImpute()函数，如下所示：

1 2	set.seed(42) data.imputed <- rfImpute(hd ~.,data=data,iter=6)

hd~.是第1个参数，它表示，我们使用hd（heart disease）这一列作为最终的预测，也就是说拿一个患者的数据，预测患者是否得有心脏病。

data = data这个表示，我们使用data这个数据集来构建随机森林。

iter=6用来指定迭代多少次来评估缺失值，从理论上讲，这个值是4到6就足够了，如下所示：

运行结果如下所示：

> data.imputed <- rfImpute(hd ~.,data=data,iter=6)
ntree      OOB      1      2
  300:  17.49% 12.80% 23.02%
ntree      OOB      1      2
  300:  16.83% 14.02% 20.14%
ntree      OOB      1      2
  300:  17.82% 13.41% 23.02%
ntree      OOB      1      2
  300:  17.49% 14.02% 21.58%
ntree      OOB      1      2
  300:  17.16% 12.80% 22.30%
ntree      OOB      1      2
  300:  18.15% 14.63% 22.30%

经过每次迭代，rflmpute()就会在第2列输出出袋错误率（OOB，即Out-of-Bag），如果估计值在不断地改善，那么这个值就会不断地变小，但在这个结果中，我们发现这个值并没有不断地变小，那么我们就可以认为，我们利用此函数得到的估计值可以满足计算，如下所示：

构建随机森林

接着，我们使用randomForest()函数来构建随机森林，如下所示：

1	model <- randomForest(hd ~., data=data.imputed, proximity= TRUE)

在这些参数中，第1个参数不用解释，跟前面的一样，第2个参数data.imputed是填充了估计值后的数据集，第3个参数是proximity = TRUE，它表示要返回相似矩阵，最终我们可以使用这个数据来进行聚类。

现在我们看一下结果，如下所示：

> model
Call:
 randomForest(formula = hd ~ ., data = data.imputed, proximity = TRUE) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 3
        OOB estimate of  error rate: 16.83%
Confusion matrix:
          Healthy Unhealthy class.error
Healthy       142        22   0.1341463
Unhealthy      29       110   0.2086331

现在说明一下这个结果，先看第一部分，如下所示：

1 2	Call: randomForest(formula = hd ~ ., data = data.imputed, proximity = TRUE)

这一部分是构建随机森林的一些参数。

再看第二部分，如下所示：

1	Type of random forest: classification

这一部分的意思是讲，我们构建随机森林的类型是分类（classification），如果我们使用随机森林来预测体重或身高这种连续型变量，那么这个就是回归（regression），如果我们忽略随机森林要预测的某个变量（我的理解就是，即不预测分类，也就是说不预测心脏病，也不预测回归，也就是说不预测体重或身高），那么这里就是无监督（unsupervised），如下所示：

再看下一部分，如下所示：

1	Number of trees: 500

这里是说，随机森林中有多少棵树，默认的数字是500，随后我们会看一下500棵树能否做出一个很好的分类。

再看下一部分，如下所示：

1	No. of variables tried at each split: 3

这一段文字是说，在每个内部节点上，使用的变量（也就是列）的数目，用于分类的随机森林默认使用变量数目的平方根进行计算，例如在这个案例中，变量的数目是13（不包括hd这个变量），那么13的平方根就是3.6，取整数（这里不是四舍五入），就是3，至于这个3是不是最好的数值，我们随后会进行调整，如果是用于回归的随机森林，那么这个数字就是变量的数目除以3，如下所示：

再看下一部分，如下所示：

1	OOB estimate of error rate: 16.83%

这一部分是出袋错误率（OOB，即Out-of-Bag），这表明，这个随机森林能够将83.5%的OOB样本正确地归类，如下所示：

最后一部分是混淆矩阵（confusion matrix），如下所示（可能与教程中的有些出入，但整体相差不大）：

Confusion matrix:
          Healthy Unhealthy class.error
Healthy       142        22   0.1341463
Unhealthy      29       110   0.2086331

混淆矩阵是总结分类模型预测结果的情形分析表，它以矩阵形式将数据集中的记录按照真实的类别与分类模型作出的分类判断两个标准进行汇总，那么，在上面的这个矩阵中，我们可以看到：

有142名健康的病人被正确地标记为“健康”；
有29名健康的病人被错误地标记为“健康”；
有22名健康的病人被错误地标记为“不健康”；
有110名不健康的病人被正确地标记为“不健康”。

绘制错误率图

此时，我们看一下，500棵树能否做出一个良好的分类，我们先画出错误率，绘图使用的是ggplot2包，因此要构建一个数据框用于ggplot2的识别，如下所示：

oob.error.data <- data.frame(
  Trees=rep(1:nrow(model$err.rate),times=3),
  Type=rep(c("OOB","Healthy","Unhealthy"),each=nrow(model$err.rate)),
  Error=c(model$err.rate[,"OOB"],
          model$err.rate[,"Healthy"],
          model$err.rate[,"Unhealthy"]))

这段代码有点复杂，现在解释一下：

这段代码的大部分都是基于model这个数据中的err.rate来构建的，我们查看下这个数据，如下所示：

> head(model$err.rate)
           OOB   Healthy Unhealthy
[1,] 0.2571429 0.2765957 0.2413793
[2,] 0.2470588 0.2068966 0.2891566
[3,] 0.2783019 0.2758621 0.2812500
[4,] 0.2704918 0.2647059 0.2777778
[5,] 0.2592593 0.2600000 0.2583333
[6,] 0.2597865 0.2467532 0.2755906
> tail(model$err.rate)
             OOB   Healthy Unhealthy
[495,] 0.1683168 0.1341463 0.2086331
[496,] 0.1683168 0.1341463 0.2086331
[497,] 0.1683168 0.1341463 0.2086331
[498,] 0.1683168 0.1341463 0.2086331
[499,] 0.1683168 0.1341463 0.2086331
[500,] 0.1683168 0.1341463 0.2086331

其中第1列是OOB错误率，第2列是健康患者错误率（就是说有多大的可能性会把健康患者错误分类），第3列是不健康患者错误率（就是说有多大的可能性会把不健康患者错误分类）。

每1行表示的是，在构建随机森林的不同阶段的的错误率，如下所示：

其中，第1行表示构建了第1棵树后的错误率，如下所示：

第2行表示，构建了前2棵树后的错误率，如下所示：

那么最后1行则表示，已经构建了随机森林中500棵树后的错误率，如下所示：

因此，使用上述代码，我们最终生成的数据框则是下面的这个样本，如下所示：

> head(oob.error.data)
  Trees Type     Error
1     1  OOB 0.2571429
2     2  OOB 0.2470588
3     3  OOB 0.2783019
4     4  OOB 0.2704918
5     5  OOB 0.2592593
6     6  OOB 0.2597865

第1列数字则表示的是随机森林中树的数目，如下所示：

第2列则是错误类型，如下所示：

再看一下其他的：

最后1列是错误的实际数值，如下所示：

调用ggplot函数，如下所示：

1 2	ggplot(data=oob.error.data,aes(x=Trees,y=Error)) + geom_line(aes(color=Type))

绘图结果如下所示：

其中蓝线是对不健康患者进行归类时的错误率，绿线是总的OOB错误率，红线是对健康患者进行归类时的错误率。通常我们发现，随着随机森林中的树越来越多，错误率是逐步下降的，如下所示

各项参数调整

添加树的数目

我们试想一下，如果我们增加更多的树，错误率会不会更低，如下所示：

此时，我们将树的数目设为1000，如下所示：

1	model <- randomForest(hd ~., data=data.imputed, ntree = 1000, proximity= TRUE)

这里我们需要注意一下，前面我们并没有使用ntree=1000这个参数，这是因为randomForest默认的ntree=500，因此这个参数如果使用默认值，就可以不添加，如果要改变这个数字，就需要更改，现在看一下结果：

Call:
 randomForest(formula = hd ~ ., data = data.imputed, ntree = 1000,      proximity = TRUE) 
               Type of random forest: classification
                     Number of trees: 1000
No. of variables tried at each split: 3
        OOB estimate of  error rate: 16.5%
Confusion matrix:
          Healthy Unhealthy class.error
Healthy       142        22   0.1341463
Unhealthy      28       111   0.2014388

看一下OOB estimate of error rate: 16.5%这一行，它的OOB值是16.5%，而当ntree=500时，这个数值是16.83%，基本上两者是相同的，我们再看一下ntree=500时的混淆矩阵，如下所示：

Confusion matrix:
          Healthy Unhealthy class.error
Healthy       142        22   0.1341463
Unhealthy      29       110   0.2086331

这两者的数值也差不多，再看一下绘图结果，如下所示：

从这个绘图结果来看，当树的数目超过500后，错误率基本上就稳定了，因此我们添加再多的树，意义不大，如下所示：

更改节点变量数目

现在我们再来看一下构建原来内部节点时使用的变量数目，原来使用的是3，如下所示：

现在我们改一下这个数字，此时我们需要一个空向量，它所含元素的数目是10，如下所示：

1	oob.values <- vector(length = 10)

然后再使用一个循环，用于测试构建每个内部节点的不同变量数目，如下所示：

代码为：

for(i in 1:10){
  temp.model <- randomForest(hd ~ .,data=data.imputed, mtry=i,ntree=1000)
  oob.values[i] <- temp.model$err.rate[nrow(temp.model$err.rate),1]
}

查看一下输出结果：

1
2
3

oob.values
[1] 0.1650165 0.1716172 0.1650165 0.1815182 0.1848185 0.1815182 0.1848185 0.1914191
[9] 0.1881188 0.1881188

其中第3个数值对应的就是mtry=3，它是构建随机森林的默认数字，从中我们可以发现，这个数字基本上是最小的（0.1716172），也就是说它有最低的OOB错误率，因此，默认数值已经是最好的数值了，不过，如果我们不知道最优的数值，那么就要尝试一下其他的数值，如下所示：

绘制随机森林的MDS图

最后，我们绘制一下随机森林的MDS图，首先利用dist()函数从1- 相似矩中构建一个距离矩阵，如下所示：

1	distance.matrix <- dist(1-model$proximity)

然后使用缩放函数cmdscale，再计算距离矩阵的x轴与y轴的占总变异的百分比，如下所示：

1 2	mds.stuff <- cmdscale(distance.matrix,eig=TRUE,x.ret=TRUE) mds.var.per <- round(mds.stuff$eig/sum(mds.stuff$eig)*100,1)

构建ggplot绘图的数据集并绘图，如下所示：

mds.values <- mds.stuff$points
mds.data <- data.frame(Sample = rownames(mds.values),
                       X=mds.values[,1],
                       Y=mds.values[,2],
                       Status=data.imputed$hd)
ggplot(data=mds.data, aes(x=X, y=Y, label=Sample)) +
  geom_text(aes(color=Status)) +
  theme_bw()+
  xlab(paste("MDS1 - ", mds.var.per[1], "%", sep="")) +
  ylab(paste("MDS2 - ", mds.var.per[2], "%", sep="")) +
  ggtitle("MDS plot using(1 - Random Forest Proximities)")

最终绘图结果如下所示：

从图上我们可以得到如下信息：

不健康的患者位于左侧，健康的患者位于右侧；
但在左侧，有几个患者被错误归类了（绿色中的红色文字）；
我们再看一下x轴，可以发现，MDS1占了总变异的47.7%，y轴上是MDS2，它占据了总异的14.2%；
整个图形呈一个V型（教程中是倒V型，我不太理解为什么我绘图出来是反的），它的开口是朝上x轴的相反方向的，这表明，在V型的两臂上，两个样本相离得越远，他们之间的差异就越大；
如果我们拿到一个新的患者的数据，在不知道他是否得了心脏病的情况下，如果经计算，他被归类到V型的左臂，那么我们就可以得到一个肯定的结论，他得了心脏病，如下所示：