1. 简介
随机森林是机器学习中的一种分类算法，这个算法在之前分享的TCGA生存模型的构建以及模型预测和评估中使用过。在介绍随机森林之前，非常有必要了解决策树这种分类器。
决策树是一种分类器，通过训练集构建一颗决策树，从而可以对新的数据预测其分类。一颗构建好的决策树如下：
可以看到这颗决策树的目标是将数据分成 "接受" 和 "拒绝" 两类，分类的条件有树中的节点来决定；
而随机森林算法，可以看作好多颗决策树构成的分类器，首先通过有放回的抽样从原始数据集中构建多个子数据集，然后利用每个子数据集构建一颗决策树，最终的分类效果由多颗决策树预测得到的众数决定；
之所以叫做随机森林，是因为两个核心观点：
1）子数据集的构建，通过随机抽样得到，所以有随机这个关键词
2）在这个分类器中，有多颗决策树，所以有森林这个关键词
随机森林的算法可以用如下几个步骤概括：
1）用有抽样放回的方法（bootstrap）从样本集中选取n个样本作为一个训练集
2）用抽样得到的样本集生成一棵决策树。在生成的每一个结点：
随机不重复地选择d个特征
3）利用这d个特征分别对样本集进行划分，找到最佳的划分特征（可用基尼系数、增益率或者信息增益判别）
4）重复步骤1到步骤2共k次，k即为随机森林中决策树的个数。
用训练得到的随机森林对测试样本进行预测，并用票选法决定预测的结果。
2. 随机森林之特征选择
随机森林具有一个重要特征：能够计算单个特征变量的重要性。并且这一特征在很多方面能够得到应用。例如在银行贷款业务中能否正确的评估一个企业的信用度，关系到是否能够有效地回收贷款。但是信用评估模型的数据特征有很多，其中不乏有很多噪音，所以需要计算出每一个特征的重要性并对这些特征进行一个排序，进而可以从所有特征中选择出重要性靠前的特征。
- 一：特征重要性
在随机森林中某个特征X的重要性的计算方法如下：
a：对于随机森林中的每一颗决策树,使用相应的OOB(袋外数据)数据来计算它的袋外数据误差,记为errOOB1.
b: 随机地对袋外数据OOB所有样本的特征X加入噪声干扰(就可以随机的改变样本在特征X处的值),再次计算它的袋外数据误差,记为errOOB2.
c：假设随机森林中有Ntree棵树,那么对于特征X的重要性=∑(errOOB2-errOOB1)/Ntree,之所以可以用这个表达式来作为相应特征的重要性的度量值是因为：若给某个特征随机加入噪声之后,袋外的准确率大幅度降低,则说明这个特征对于样本的分类结果影响很大,也就是说它的重要程度比较高。
- 二：特征选择
在论文 Variable Selection using Random Forests中详细的论述了基于随机森林的特征选择方法，这里我们进行一些回顾。
首先特征选择的目标有两个：
1）找到与应变量高度相关的特征变量。
2）选择出数目较少的特征变量并且能够充分的预测应变量的结果。
其次一般特征选择的步骤为：
1）初步估计和排序
a: 对随机森林中的特征变量按照VI（Variable Importance）降序排序。
b: 确定删除比例,从当前的特征变量中剔除相应比例不重要的指标，从而得到一个新的特征集。
c: 用新的特征集建立新的随机森林,并计算特征集中每个特征的VI,并排序。
d: 重复以上步骤,直到剩下m个特征。
2）根据1中得到的每个特征集和它们建立起来的随机森林,计算对应的袋外误差率(OOB err),将袋外误差率最低的特征集作为最后选定的特征集。
来源：随机森林之特征选择
3. randomForest 包的使用
randomForest 包提供了利用随机森林算法解决分类和回归问题的功能；我们这里只关注随机森林算法在分类问题中的应用
安装randomForest包
install.packages("randomForest")
library(randomForest)
载入iris演示数据集
iris 数据集共有150行，5列，其中第5列Species为分类变量，共有3种分类情况。
这个数据集可以看做有150个样本，将Species这个分类变量作为因变量，使用其他4个指标来预测花属于哪个Species。
data(iris)
str(iris)
# 'data.frame': 150 obs. of 5 variables:
# $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
# $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
# $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
# $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
# $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
head(iris)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species
# 1 5.1 3.5 1.4 0.2 setosa
# 2 4.9 3.0 1.4 0.2 setosa
# 3 4.7 3.2 1.3 0.2 setosa
# 4 4.6 3.1 1.5 0.2 setosa
# 5 5.0 3.6 1.4 0.2 setosa
# 6 5.4 3.9 1.7 0.4 setosa
接下来调用randomForest 函数进行分类
调用该函数时，通过一个表达式指定分类变量 Species 和对应的数据集data 就可以了。Species ~ .的意思是iris数据集中以Species为被解释变量，其他列为解释变量。后面的importance 和 proximity 是计算每个变量的重要性和样本之间的距离
set.seed(71)
iris.rf <- randomForest(Species ~ ., data=iris, importance=TRUE,
proximity=TRUE)
分类器构建完毕之后，首先看一下这个分类器的准确性
print(iris.rf)
# Call:
# randomForest(formula = Species ~ ., data = iris, importance = TRUE, proximity = TRUE)
# Type of random forest: classification
# Number of trees: 500
# No. of variables tried at each split: 2
# OOB estimate of error rate: 4.67%
# Confusion matrix:
# setosa versicolor virginica class.error
# setosa 50 0 0 0.00
# versicolor 0 47 3 0.06
# virginica 0 4 46 0.08
print 的结果中，OOB estimate of error rate 表明了分类器的错误率为4.67%， Confusion matrix 表明了每个分类的详细的分类情况；
对于setosa 这个group而言，基于随机森林算法的分类器，有50个样本分类到了setosa 这个group, 而且这50个样本和iris 中属于setosa 这个group的样本完全一致，所以对于setosa 这个group而言，分类器的错误率为0；
对于versicolor 这个group而言，基于随机森林算法的分类器，有47个样本分类到了versicolor 这个group, 3个样本分类到了virginica 这个group，有3个样本分类错误，在iris 中属于versicolor 这个group的样本有50个，所以对于versicolor 这个group而言，分类器的错误率为3/50 = 0.06 ；
对于virginica 这个group而言，基于随机森林算法的分类器，有3个样本分类到了versicolor 这个group, 47个样本分类到了virginica 这个group，有3个样本分类错误，在iris 中属于virginica 这个group的样本有50个，所以对于virginica这个group而言，分类器的错误率为3/50 = 0.06 ；
然后看一下样本之间的距离
iris.mds <- cmdscale(1 - iris.rf$proximity, eig=TRUE)
通过调用cmdscale 函数进行样本之间的距离，proximity 是样本之间的相似度矩阵，所以用1减去之后得到样本的类似距离矩阵的一个矩阵
iris.mds 的结果如下
str(iris.mds)
# List of 5
# $ points: num [1:150, 1:2] -0.565 -0.565 -0.566 -0.566 -0.565 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:150] "1" "2" "3" "4" ...
# .. ..$ : NULL
# $ eig : num [1:150] 24 21.15 2.3 1.34 1.18 ...
# $ x : NULL
# $ ac : num 0
# $ GOF : num [1:2] 0.737 0.799
head(iris.mds$points)
# [,1] [,2]
# 1 -0.5649559 -0.04385526
# 2 -0.5651790 -0.04419585
# 3 -0.5656595 -0.04355633
# 4 -0.5657237 -0.04388197
# 5 -0.5651469 -0.04390029
# 6 -0.5655295 -0.04356653
在iris.mds 中points可以看做每个样本映射到2维空间中的坐标，每一维空间是一个分类特征，但是不是最原始的4个特征，而是由4个特征衍生得到的新的分类特征，根据这个坐标，可以画一张散点图，得到每个样本基于两个分类变量的分组情况
plot(iris.mds$points, col = rep(c("red", "blue", "green"), each = 50))
图中不同分类的样本用不同的颜色标注，可以看到基于两个新的分类特征，样本的分组效果还是很好的，不同组的样本明显区分开来。
最后，在看一下4个特征，每个特征的重要性
iris.rf$importance
# setosa versicolor virginica MeanDecreaseAccuracy MeanDecreaseGini
# Sepal.Length 0.02837234 0.0169917984 0.043822619 0.029560910 9.371080
# Sepal.Width 0.01019719 0.0007153689 0.008415053 0.006240133 2.446698
# Petal.Length 0.32043066 0.2938837690 0.284731709 0.297205369 42.134304
# Petal.Width 0.34286020 0.3205494688 0.281664051 0.312822368 45.284763
之前调用randomForest 函数时，通过指定importance = TRUE 来计算每个特征的importance , 在 iris.rf$importance 矩阵中，有两个值是需要重点关注的MeanDecreaseAccuracy 和 MeanDecreaseGini
我们还可以利用
varImpPlot(iris.rf, main = "Top 30 - variable importance")
图中和坐标为importance 结果中的MeanDecreaseAccuracy 和 MeanDecreaseGini 指标的值，纵坐标为对应的每个分类特征，该函数默认画top30个特征，由于这个数据集只有4个分类特征，所以4个都出现了。
⚠️⚠️⚠️使用randomForest对转录组的差异基因进行进一步的筛选
比如转录组测序中最后筛选出12个差异基因，想要对他们进行重要度的排序，需要把expr表格整理成如下格式：前面12列是12个基因的表达量，必须是数值型。最后一列是分组，必须是factor型。然后就可以建模
rf <- randomForest(group ~ ., data=expr, importance=TRUE, proximity=TRUE) varImpPlot(rf, main = "Variable importance")
4. 划分训练、测试集
前面的演示中，150个样本全部作为了训练集，在这里演示一下划分了训练、测试集的模型
4.1 载入包和演示数据集
library(randomForest)
library(caret)
library(pROC)
data(iris)
4.2 划分训练、测试集（一般来说比例是8:2，或者7:3）
⚠️训练集和测试集的变量的数目必须是一致的，否则会报错（这里都是5）。
# 80%是训练集，20%是测试集。
trainlist <- createDataPartition(iris$Species,p=0.8,list = FALSE)
trainset <- iris[trainlist,]
testset <- iris[-trainlist,]
dim(trainset)
# [1] 120 5
dim(testset)
# [1] 30 5
4.3 使用训练集建模
set.seed(6666)
rf.train <- randomForest(as.factor(Species) ~ ., data=trainset, importance=TRUE,
proximity=TRUE,na.action = na.pass)
rf.train
# Call:
# randomForest(formula = as.factor(Species) ~ ., data = trainset, importance = TRUE, proximity = TRUE, na.action = na.pass)
# Type of random forest: classification
# Number of trees: 500
# No. of variables tried at each split: 2
# OOB estimate of error rate: 6.67%
# Confusion matrix:
# setosa versicolor virginica class.error
# setosa 40 0 0 0.000
# versicolor 0 37 3 0.075
# virginica 0 5 35 0.125
plot(rf.train,main = 'randomforest origin') #绘图查看训练效果
4.4 预测测试集
set.seed(6666)
rf.test <- predict(rf.train,newdata = testset,type = 'class')
# 11 19 20 21 31 35 37 44
# setosa setosa setosa setosa setosa setosa setosa setosa
# 47 49 60 63 65 66 69 75
# setosa setosa versicolor versicolor versicolor versicolor versicolor versicolor
# 76 79 93 97 105 106 108 109
# versicolor versicolor versicolor versicolor virginica virginica virginica virginica
# 112 113 116 118 121 142
# virginica virginica virginica virginica virginica virginica
# Levels: setosa versicolor virginica
统计预测结果
rf.cf <- caret::confusionMatrix(as.factor(rf.test),as.factor(testset$Species))
4.5 ROC, AUC
首先要注意的是，ROC需要的是概率，而不是预测的种类。所以要再预测一测，把type改成prob
rf.test2 <- predict(rf.train,newdata = testset,type = 'prob')
head(rf.test2)
# setosa versicolor virginica
# 11 1.00 0.00 0
# 19 0.98 0.02 0
# 20 1.00 0.00 0
# 21 1.00 0.00 0
# 31 1.00 0.00 0
# 35 1.00 0.00 0
绘制roc（因为我们有3个class，所以使用的是multiclass.roc，多种类的ROC）
roc.rf <- multiclass.roc(testset$Species,rf.test2)
roc.rf
# Call:
# multiclass.roc.default(response = testset$Species, predictor = rf.test2)
# Data: multivariate predictor rf.test2 with 3 levels of testset$Species: setosa, versicolor, virginica.
# Multi-class area under the curve: 1
可以看到，AUC=1。
