1. 簡(jiǎn)介

隨機(jī)森林是機(jī)器學(xué)習(xí)中的一種分類算法，這個(gè)算法在之前分享的TCGA生存模型的構(gòu)建以及模型預(yù)測(cè)和評(píng)估中使用過。在介紹隨機(jī)森林之前，非常有必要了解決策樹這種分類器。
決策樹是一種分類器，通過訓(xùn)練集構(gòu)建一顆決策樹，從而可以對(duì)新的數(shù)據(jù)預(yù)測(cè)其分類。一顆構(gòu)建好的決策樹如下：

可以看到這顆決策樹的目標(biāo)是將數(shù)據(jù)分成 "接受" 和 "拒絕" 兩類，分類的條件有樹中的節(jié)點(diǎn)來決定；

而隨機(jī)森林算法，可以看作好多顆決策樹構(gòu)成的分類器，首先通過有放回的抽樣從原始數(shù)據(jù)集中構(gòu)建多個(gè)子數(shù)據(jù)集，然后利用每個(gè)子數(shù)據(jù)集構(gòu)建一顆決策樹，最終的分類效果由多顆決策樹預(yù)測(cè)得到的眾數(shù)決定；

之所以叫做隨機(jī)森林，是因?yàn)閮蓚€(gè)核心觀點(diǎn)：
1）子數(shù)據(jù)集的構(gòu)建，通過隨機(jī)抽樣得到，所以有隨機(jī)這個(gè)關(guān)鍵詞
2）在這個(gè)分類器中，有多顆決策樹，所以有森林這個(gè)關(guān)鍵詞

隨機(jī)森林的算法可以用如下幾個(gè)步驟概括：
1）用有抽樣放回的方法（bootstrap）從樣本集中選取n個(gè)樣本作為一個(gè)訓(xùn)練集
2）用抽樣得到的樣本集生成一棵決策樹。在生成的每一個(gè)結(jié)點(diǎn)：
隨機(jī)不重復(fù)地選擇d個(gè)特征
3）利用這d個(gè)特征分別對(duì)樣本集進(jìn)行劃分，找到最佳的劃分特征（可用基尼系數(shù)、增益率或者信息增益判別）
4）重復(fù)步驟1到步驟2共k次，k即為隨機(jī)森林中決策樹的個(gè)數(shù)。
用訓(xùn)練得到的隨機(jī)森林對(duì)測(cè)試樣本進(jìn)行預(yù)測(cè)，并用票選法決定預(yù)測(cè)的結(jié)果。

楊凱, 侯艷, 李康. 隨機(jī)森林變量重要性評(píng)分及其研究進(jìn)展[J]. 2015.

2. 隨機(jī)森林之特征選擇

隨機(jī)森林具有一個(gè)重要特征：能夠計(jì)算單個(gè)特征變量的重要性。并且這一特征在很多方面能夠得到應(yīng)用。例如在銀行貸款業(yè)務(wù)中能否正確的評(píng)估一個(gè)企業(yè)的信用度，關(guān)系到是否能夠有效地回收貸款。但是信用評(píng)估模型的數(shù)據(jù)特征有很多，其中不乏有很多噪音，所以需要計(jì)算出每一個(gè)特征的重要性并對(duì)這些特征進(jìn)行一個(gè)排序，進(jìn)而可以從所有特征中選擇出重要性靠前的特征。

• 一：特征重要性

在隨機(jī)森林中某個(gè)特征X的重要性的計(jì)算方法如下：
a：對(duì)于隨機(jī)森林中的每一顆決策樹,使用相應(yīng)的OOB(袋外數(shù)據(jù))數(shù)據(jù)來計(jì)算它的袋外數(shù)據(jù)誤差,記為errOOB1.
b: 隨機(jī)地對(duì)袋外數(shù)據(jù)OOB所有樣本的特征X加入噪聲干擾(就可以隨機(jī)的改變樣本在特征X處的值),再次計(jì)算它的袋外數(shù)據(jù)誤差,記為errOOB2.
c：假設(shè)隨機(jī)森林中有Ntree棵樹,那么對(duì)于特征X的重要性=∑(errOOB2-errOOB1)/Ntree,之所以可以用這個(gè)表達(dá)式來作為相應(yīng)特征的重要性的度量值是因?yàn)椋喝艚o某個(gè)特征隨機(jī)加入噪聲之后,袋外的準(zhǔn)確率大幅度降低,則說明這個(gè)特征對(duì)于樣本的分類結(jié)果影響很大,也就是說它的重要程度比較高。

• 二：特征選擇

在論文 Variable Selection using Random Forests中詳細(xì)的論述了基于隨機(jī)森林的特征選擇方法，這里我們進(jìn)行一些回顧。

首先特征選擇的目標(biāo)有兩個(gè)：
1）找到與應(yīng)變量高度相關(guān)的特征變量。
2）選擇出數(shù)目較少的特征變量并且能夠充分的預(yù)測(cè)應(yīng)變量的結(jié)果。

其次一般特征選擇的步驟為：
1）初步估計(jì)和排序
a: 對(duì)隨機(jī)森林中的特征變量按照VI（Variable Importance）降序排序。
b: 確定刪除比例,從當(dāng)前的特征變量中剔除相應(yīng)比例不重要的指標(biāo)，從而得到一個(gè)新的特征集。
c: 用新的特征集建立新的隨機(jī)森林,并計(jì)算特征集中每個(gè)特征的VI,并排序。
d: 重復(fù)以上步驟,直到剩下m個(gè)特征。
2）根據(jù)1中得到的每個(gè)特征集和它們建立起來的隨機(jī)森林,計(jì)算對(duì)應(yīng)的袋外誤差率(OOB err),將袋外誤差率最低的特征集作為最后選定的特征集。

來源：隨機(jī)森林之特征選擇

3. randomForest 包的使用

randomForest 包提供了利用隨機(jī)森林算法解決分類和回歸問題的功能；我們這里只關(guān)注隨機(jī)森林算法在分類問題中的應(yīng)用

安裝randomForest包

install.packages("randomForest")
library(randomForest)

載入iris演示數(shù)據(jù)集
iris 數(shù)據(jù)集共有150行，5列，其中第5列Species為分類變量，共有3種分類情況。
這個(gè)數(shù)據(jù)集可以看做有150個(gè)樣本，將Species這個(gè)分類變量作為因變量，使用其他4個(gè)指標(biāo)來預(yù)測(cè)花屬于哪個(gè)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

接下來調(diào)用randomForest 函數(shù)進(jìn)行分類
調(diào)用該函數(shù)時(shí)，通過一個(gè)表達(dá)式指定分類變量 Species 和對(duì)應(yīng)的數(shù)據(jù)集data 就可以了。Species ~ .的意思是iris數(shù)據(jù)集中以Species為被解釋變量，其他列為解釋變量。后面的importance 和 proximity 是計(jì)算每個(gè)變量的重要性和樣本之間的距離

set.seed(71)
iris.rf <- randomForest(Species ~ ., data=iris, importance=TRUE,
                             proximity=TRUE)

分類器構(gòu)建完畢之后，首先看一下這個(gè)分類器的準(zhǔn)確性

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 的結(jié)果中，OOB estimate of error rate 表明了分類器的錯(cuò)誤率為4.67%， Confusion matrix 表明了每個(gè)分類的詳細(xì)的分類情況；

對(duì)于setosa 這個(gè)group而言，基于隨機(jī)森林算法的分類器，有50個(gè)樣本分類到了setosa 這個(gè)group, 而且這50個(gè)樣本和iris 中屬于setosa 這個(gè)group的樣本完全一致，所以對(duì)于setosa 這個(gè)group而言，分類器的錯(cuò)誤率為0；

對(duì)于versicolor 這個(gè)group而言，基于隨機(jī)森林算法的分類器，有47個(gè)樣本分類到了versicolor 這個(gè)group, 3個(gè)樣本分類到了virginica 這個(gè)group，有3個(gè)樣本分類錯(cuò)誤，在iris 中屬于versicolor 這個(gè)group的樣本有50個(gè)，所以對(duì)于versicolor 這個(gè)group而言，分類器的錯(cuò)誤率為3/50 = 0.06 ；

對(duì)于virginica 這個(gè)group而言，基于隨機(jī)森林算法的分類器，有3個(gè)樣本分類到了versicolor 這個(gè)group, 47個(gè)樣本分類到了virginica 這個(gè)group，有3個(gè)樣本分類錯(cuò)誤，在iris 中屬于virginica 這個(gè)group的樣本有50個(gè)，所以對(duì)于virginica這個(gè)group而言，分類器的錯(cuò)誤率為3/50 = 0.06 ；

然后看一下樣本之間的距離

iris.mds <- cmdscale(1 - iris.rf$proximity, eig=TRUE)

通過調(diào)用cmdscale 函數(shù)進(jìn)行樣本之間的距離，proximity 是樣本之間的相似度矩陣，所以用1減去之后得到樣本的類似距離矩陣的一個(gè)矩陣
iris.mds 的結(jié)果如下

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可以看做每個(gè)樣本映射到2維空間中的坐標(biāo)，每一維空間是一個(gè)分類特征，但是不是最原始的4個(gè)特征，而是由4個(gè)特征衍生得到的新的分類特征，根據(jù)這個(gè)坐標(biāo)，可以畫一張散點(diǎn)圖，得到每個(gè)樣本基于兩個(gè)分類變量的分組情況

plot(iris.mds$points, col = rep(c("red", "blue", "green"), each = 50))

圖中不同分類的樣本用不同的顏色標(biāo)注，可以看到基于兩個(gè)新的分類特征，樣本的分組效果還是很好的，不同組的樣本明顯區(qū)分開來。

最后，在看一下4個(gè)特征，每個(gè)特征的重要性

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

之前調(diào)用randomForest 函數(shù)時(shí)，通過指定importance = TRUE 來計(jì)算每個(gè)特征的importance , 在 iris.rf$importance 矩陣中，有兩個(gè)值是需要重點(diǎn)關(guān)注的MeanDecreaseAccuracy 和 MeanDecreaseGini
我們還可以利用

varImpPlot(iris.rf, main = "Top 30 - variable importance")

圖中和坐標(biāo)為importance 結(jié)果中的MeanDecreaseAccuracy 和 MeanDecreaseGini 指標(biāo)的值，縱坐標(biāo)為對(duì)應(yīng)的每個(gè)分類特征，該函數(shù)默認(rèn)畫top30個(gè)特征，由于這個(gè)數(shù)據(jù)集只有4個(gè)分類特征，所以4個(gè)都出現(xiàn)了。

??????使用randomForest對(duì)轉(zhuǎn)錄組的差異基因進(jìn)行進(jìn)一步的篩選
比如轉(zhuǎn)錄組測(cè)序中最后篩選出12個(gè)差異基因，想要對(duì)他們進(jìn)行重要度的排序，需要把expr表格整理成如下格式：

前面12列是12個(gè)基因的表達(dá)量，必須是數(shù)值型。最后一列是分組，必須是factor型。

然后就可以建模

rf <- randomForest(group ~ ., data=expr, importance=TRUE, proximity=TRUE)
varImpPlot(rf, main = "Variable importance")

4. 劃分訓(xùn)練、測(cè)試集

前面的演示中，150個(gè)樣本全部作為了訓(xùn)練集，在這里演示一下劃分了訓(xùn)練、測(cè)試集的模型

4.1 載入包和演示數(shù)據(jù)集

library(randomForest)
library(caret)
library(pROC)
data(iris)

4.2 劃分訓(xùn)練、測(cè)試集（一般來說比例是8:2，或者7:3）

??訓(xùn)練集和測(cè)試集的變量的數(shù)目必須是一致的，否則會(huì)報(bào)錯(cuò)（這里都是5）。

# 80%是訓(xùn)練集，20%是測(cè)試集。
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 使用訓(xùn)練集建模

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') #繪圖查看訓(xùn)練效果

黑色的線是均值，可以看到50之后，幾條彩色的線都可以分的比較開，說明50棵樹以上都可以分的比較好。（樹的數(shù)量默認(rèn)是500）

4.4 預(yù)測(cè)測(cè)試集

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

統(tǒng)計(jì)預(yù)測(cè)結(jié)果

rf.cf <- caret::confusionMatrix(as.factor(rf.test),as.factor(testset$Species))

從混淆矩陣中可以看出來，對(duì)這30個(gè)樣本的預(yù)測(cè)是全對(duì)的，說明模型構(gòu)建的很不錯(cuò)。Accuracy和Kappa值越接近于1，預(yù)測(cè)效果越好。

4.5 ROC, AUC

首先要注意的是，ROC需要的是概率，而不是預(yù)測(cè)的種類。所以要再預(yù)測(cè)一測(cè)，把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（因?yàn)槲覀冇?個(gè)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。