2020-01-15 随机森林-原理及如何用R绘图

参考学习资料：菜鸟团生信专题

1. StatQuest生物统计学专题 - 决策树(1)
2. StatQuest生物统计学专题 - 决策树(2)
3. StatQuest生物统计学专题 - 随机森林(1) 构建与评价
4. StatQuest生物统计学专题 - 随机森林(2) R实例

对统计学原理不太熟练的时候听到一个名词随机森林，这是个什么鬼，这个森林里有啥，随机出现什么？这东西干啥的？啥也不知道啊，只好学习了。
根据搜到的结果终于知道了这个统计方法是怎么回事了。

随机森林就是使用随机的方法重新构建大量的决策树，根据多个决策树分类结果的“少数服从多数”的原则来完成对新数据的分类预测。

果然是森林，由很多树组成的，很多什么树呢，决策树，厉害了，这个树是有决策能力的的，能预测森林的走向，那么怎么实现的呢，决策树又是怎么来的呢，是来源于一堆数据，理清楚了，就是一堆数据通过某种算法变成了由决策树的形式展现出来，然后很多树对这个数据进行了深度分析，然后形成了一个随机森林可以根据已有的这堆数据对某个现象或某个疾病的进展做一个预后的推测，通常情况下预测精度还是蛮精准的。
这不就是赋予机器以思想的过程吗，估计深度学习就是这么来的。高级！
以上是我脑补的，具体过程还请看前面的参考资料。

随机森林的构建方法

构建一个决策树是很简单的，分为3步：

• 构建bootstrap数据集；
• 创建bootstrap数据集的决策树；
• 如此重复很多次，树多成林，就完成了随机森林的构建

那么重点来了

如何构建bootstrap数据集？

• bootstrap数据集就是以有放回的方式重新从原数据中抽取同样数量的记录构成新的数据集。
• 比如原始数据有4条记录，因此要随机抽取4条组成新的记录，假定我们随机抽取四次，结果为第2条、第1条、第4条、第4条。

什么是决策树？

决策树是一种有监督机器学习算法，所以需要一个已知数据用于计算，这种数据分为两部分：我们可以分别将其理解为“自变量”和“因变量”。

• 决策树的构建方法：
  • 首先从根节点开始，逐一尝试所有选项（自变量），并选择一个最佳的选项作为根节点；
  • 然后限定根节点，再逐一尝试剩余选项（自变量），并选择一个最佳的选项作为内部节点；

第一步——得到最佳的根节点
使用一种叫做Gini的方法，Gini用于衡量每一个决策树的impurity（不纯净度），Gini值越低越好，它的计算非常简单：
Gini impurity=1-(1-probability of "Yes")^2-(1-probability of "No")^2
第二步——限定根节点，同样方法找到最佳的内部节点
第三步——继续限定二级节点，完成最终的一侧分支
以此类推，直到完成全部的节点构建。（具体实例参考上述学习资料）

构建每一个bootstrap数据集的决策树

• 构建每一个bootstrap数据集的决策树，不同于决策树方法，此处要加一个调整，就是每次只选择特定数量的变量子集。
  如此重复很多次，就完成了随机森林的构建

• 将上述步骤重复很多次，比如是100次，那么这100个决策树就构成了随机森林。

• 随机森林的用法也很简单，当有新数据需要分类时，那么100个决策树都要纳入计算，根据100个分类结果，按照“少数服从多数”的原则进行定论，假如100个分类结果有80个Yes、20个No，那么最终分类结果就是”Yes“。

随机森林的模型评价方法

随机森林通过bootstrap和随机子集的方法可以产生大量的随机数据集，从而形成随机森林。而在bootstrap的过程中，由于是有放回的抽样，所以每一个boots数据集中绝大部分都会有重复样本，也就是说原始数据集中的一些样本不在bootstrap数据集中，这些样本称为Out-of-bag Dataset。

在每一次bootstrap过程中，Out-of-bag Dataset的样本数量大约为原始数据量的1/3（0.36左右）。

反过来说，假定一个随机森林有100个决策树，对于原始数据来说，每一个样本都有可能是某些树的Out-of-bag数据。

所以刚好就可以使用Out-of-bag Dataset数据集去评价随机森林的效果好坏。

具体如何评价呢？

方法就是将同Out-of-bag数据对应的决策树对Out-of-bag数据进行分类计算，看计算出来的分类结果和原始分类是否相符，计算不相符的Out-of-bag Dataset的比例，此比例就是随机森林的优劣程度评价。因此又叫做Out-of-bag Error。

最优子集数目的选择

由于我们已经知道了如何构建随机森林，也已经知道了随机森林的评价方法，所以寻找最优子集数目是一个非常简单的过程，多试试几个子集数目，寻找Out-of-bag Error最小的随机森林即可。

最优子集的选择是如下过程：

• 一般先从“样本数目的平方根”个数目开始，如过样本数目为4，则先定子集数目为2；

• 计算此时的随机森林和相应的Out-of-bag Error值；

• 然后再将“样本数目的平方根”周围的一些数作为子集数目，计算此时的随机森林和Out-of-bag Error；

• 根据Out-of-bag Error值，选择最佳的子集数目和相应的随机森林即可。

好，原理理清楚了，那就在R中进行随机森林的实践。

原始数据是使用年龄、性别、胸痛程度、血压等共14个参数来预测一个人是否罹患心脏病。原始数据位于http://archive.ics.uci.edu/ml/datasets/Heart+Disease。

用到的R包为ggplot2、randomForest。

library(ggplot2)
library(randomForest)

下载数据

使用read.csv函数读取网络数据。

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

数据清洗

修改列名

# 取回的数据没有列名
head(data)
# 手动修改列名
colnames(data) <- c(
  "age",  "sex",  "cp",   "trestbps", 
  "chol",   "fbs",    "restecg", 
  "thalach",   "exang",    "oldpeak", 
  "slope",  "ca",   "thal",  "hd" 
) # 共14个
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

修改数据格式

从官网可以知道14个数据变量的含义，他们分别是：

  age
  sex: 0 = female, 1 = male
  cp: chest pain
          # 1 = typical angina,
          # 2 = atypical angina,
          # 3 = non-anginal pain,
          # 4 = asymptomatic
  trestbps: resting blood pressure (in mm Hg)
  chol: serum cholestoral in mg/dl
  fbs: fasting blood sugar greater than 120 mg/dl, 1 = TRUE, 0 = FALSE
  restecg: resting electrocardiographic results
          # 1 = normal
          # 2 = having ST-T wave abnormality
          # 3 = showing probable or definite left ventricular hypertrophy
  thalach: maximum heart rate achieved
  exang: exercise induced angina, 1 = yes, 0 = no
  oldpeak: ST depression induced by exercise relative to rest
  slope: the slope of the peak exercise ST segment
          # 1 = upsloping
          # 2 = flat
          # 3 = downsloping
  ca: number of major vessels (0-3) colored by fluoroscopy
  thal： this is short of thalium heart scan
          # 3 = normal (no cold spots)
          # 6 = fixed defect (cold spots during rest and exercise)
          # 7 = reversible defect (when cold spots only appear during exercise)
  hd： (the predicted attribute) - diagnosis of heart disease
          # 0 if less than or equal to 50% diameter narrowing
          # 1 if greater than 50% diameter narrowing

然后我们使用str(data)看一下此时的数据结构可以发现不规范的变量有，

> 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代表男；
• cp、fbs、restecg、exang、slope：应该是因子变量；
• ca：应该使用“NA”表示缺失值，而不是“？”；
• thal：和ca情况一样；
• hd：因子变量，0是健康，1是心脏病；

按照上述情况，进行以下修改：

#先将“？”更换为“NA”
data[data=="?"]<-NA

#将sex修改为正确的因子变量，并处理好性别
data$sex <- factor(data$sex,levels = c(0,1), labels = c("F","M"))
#将cp、fbs、restecg、exang、slope修改为因子变量
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)

#由于ca刚开始使用的’？‘代表的缺失值，因此ca里面的变量都是字符类型的因子变量
#因此需要先转换为数值型，再进行因子变量转换
data$ca <- as.integer(data$ca)
data$ca <- as.factor(data$ca)
#同理，将thal做同样处理
data$thal <- as.integer(data$thal) 
data$thal <- as.factor(data$thal)

#将hd修改为因子变量，0为健康，1为不健康
data$hd <- ifelse(data$hd == 0, "Healthy", "Unhealthy")
data$hd <- as.factor(data$hd) 

转换数据格式，这里注意as.的用法，以及函数ifelse的用法

构建随机森林

填补缺失值

rfImpute函数用于填补缺失值，随机森林的缺失值填补是根据相似度进行填补的一种迭代算法。

简单来说，缺失值会首先根据数据类别不同进行填充，也就是数值数据填充当前变量的中位数，分类数据填充当前变量的众数。然后构建随机森林，并根据随机森林来决定所有的记录之间的相似度，最后根据相似度和当前变量的其他数据对缺失值进行填补。

• 一般上述过程会迭代4-6次，4-6次被认为是比较合理的迭代次数。
• R中ifImpute函数默认创建300棵决策树的随机森林。
• 最佳子集数目根据数据类别不同进行设定，数值数据为总变量数除以3，分类数据为总变量数的平方根。

# rfImpute填补缺失值
# hd~. 其中hd为相应变量，.为其他所有变量，其意义为使用其他所有变量预测
hddata.imputed <- rfImpute(hd ~ ., data = data, iter=6)

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

结果会输出每次迭代后的OOB值，越低越好。

构建随机森林

# 构建随机森林
# 默认创建500棵决策树的随机森林。
# 最佳子集数目根据数据类别不同进行设定，数值数据为总变量数除以3，分类数据为总变量数的平方根。

model <- randomForest(hd ~ ., data=data.imputed, proximity=TRUE)
# proximity参数不是必须的，加上后，则会输出proximity矩阵，此矩阵可用于热图或MDS（PCoA）

随机森林评价

决策树的数量

默认是创建500棵决策树，此时的OOB（out of bag）值可以用于评价随机森林的模型如何。

我们可以看看此时从第1棵树到第500棵决策树时，OOB的变化趋势：

# model下err.rate中是OOB的数据，它有三列，分别是总OOB、健康人的OOB以及不健康人的OOB
# 创建数据框用于ggplot绘图
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"]))
# 绘图
ggplot(oob.error.data, aes(x=Trees, y=Error)) + geom_line(aes(color=Type))

可以看出，大概从150以后的OOB的值趋于稳定了，默认的500是非常稳健的数值了。

Error
最佳子集数目
默认子集数目为总变量数的平方跟，也就是13的平方根，约为3.6，所以默认的子集数目为3。

我们可以改变不同的子集数目以确认最佳子集数目是多少，比如可以看一下子集数目分别为1-10时的结果：

> oob.values <- vector(length=10)
> for (i in 1:10){
+   temp.model <- randomForest(hd~.,data = data.imputed,mtry=i)
+   oob.values[i] <- temp.model$err.rate[nrow(temp.model$err.rate),1]
+ }
> oob.values
 [1] 0.1683168 0.1749175 0.1650165 0.1815182 0.1782178 0.1881188 0.1815182
 [8] 0.1947195 0.1980198 0.1947195

可以发现最低的OOB确实是子集数目为3。