|机器学习| - R语言 - 随机森林 - 分类、回归、预测、筛选变量有史以来超详细解析（干货满满）

研究如何用R去实现随机森林也有三个月的时间了，从一开始的完全不理解，到现在的游刃有余，我似乎花了过多的时间，毕竟是初学者嘛。不知各位有没有发现，CSDN上随机森林的教程都说的有些模糊，好像在刻意回避着什么，而且很少有人说变量筛选的问题。所以，今日，我觉得有必要记录一下了。

随机森林 基于R

你即将从这里看到

1.如何使用随机森林做回归、分类、预测、交叉验证
2.如何评价模型效能（ROC、AUC）
3.如何筛选变量

在这里你不会看到

1.长篇大段没有注释的代码
2.故弄玄虚，简单问题复杂化
3.重复造轮子（本文以解决实际问题为主）

分类与回归

数据科学最基础的问题，分类与回归。
而随机森林做分类还是回归是取决于label是factor类型的还是numerical类型的

1.导包
一般需要四个包，如下，第一个是导入数据用的，第二个是随机森林包，第三个是筛选变量用的，第四个是绘制ROC曲线、计算AUC得分的包

library(xlsx)
library(randomForest)
library(varSelRF)
library(pROC)

2.数据读取、缺失值处理
读取数据，此数据是我做社会调查自己收集来的数据，各位可替换成鸢尾花数据集

d<-read.xlsx('socialrecode2.xlsx',1)#这种直接引用文件名的写法请确保文件在工作空间内
									#参数1代表使用第一行作为列名
d1<-na.omit(d)			#删去缺失值，随机森林分类器内不能传入缺失值
#如果手头没有待处理的数据，请使用鸢尾花数据集
d1<-data(iris)		#鸢尾花数据集是无缺失值的，故不需删去缺失值

如果你的数据量不够多，但是又含有许多缺失值，建议使用randomforest包自带的缺失值填充函数：na.roughfix(object, …)
按中位数或者模型进行填充
3.训练模型
在直觉上，似乎一片森林中决策树越多越好，但真的是这样吗？（此处有疑问不要紧，可参考文末原理部分）
为了测试出多少棵树时模型是最优的，我们先试探性的用800棵树试一下，即第二部分的代码。然后画出错误率与决策树的数量的关系曲线

#按index拆分训练集与测试集
index <- sample(2,nrow(iris),replace = TRUE,prob=c(0.7,0.3))
traindata <- iris[index==1,]
testdata <- iris[index==2,]

#试验性训练模型
set.seed(1234)#设置随机种子数，确保以后再执行代码时可以得到一样的结果
rf_ntree<- randomForest(Species ~ ., data=traindata,  #我这里把代码分成两行来写了
						ntree=800,important=TRUE,proximity=TRUE)
plot(rf_ntree)

我们可以看出，在树的数量大于200的时候就基本稳定了，在500左右会浮动一下，在600以后又是稳定的，而且错误率水平其实差不多，所以我们避开500左右的取值，在保证效能的情况下减少决策树的数量，减少运行时间。在这里我们就取ntree=200好了。

#正式训练
iris_rf <- randomForest(Species ~ ., data=traindata,
						ntree=200,important=TRUE,proximity=TRUE)

在这里需要注意的是，,important=TRUE意为保留指标重要性、proximity=TRUE意为保留近似矩阵，这两个参数分别在筛选变量、数据可视化的时候起很大作用，一定不可漏掉。
当然也有另一种写法，这种写法比较贴近我的习惯,谁是feature谁是label一目了然

iris_rf <- randomForest(x=traindata[,1:4], y=traindata[,5], 
						ntree=200,important=TRUE,proximity=TRUE)

然后我们查看模型具体情况

iris_rf

###输出结果如下（这是输出结果，不是代码）
Call:
 randomForest(formula = Species ~ ., data = traindata, ntree = 200,      important = TRUE, proximity = TRUE) 
               Type of random forest: classification
                     Number of trees: 200
No. of variables tried at each split: 2

        OOB estimate of  error rate: 5.66%
Confusion matrix:
           setosa versicolor virginica class.error
setosa         35          0         0  0.00000000
versicolor      0         32         3  0.08571429
virginica       0          3        33  0.08333333

4.预测
很简单，直接用predict函数就可以，第一个参数可以理解为“以什么样的方式进行预测”，第二个参数就是你要根据什么数据进行预测

iris_pred<-predict(iris_rf, newdata=testdata)
table(iris_pred,testdata$Species)

#输出结果如下
#iris_pred    setosa versicolor virginica
  #setosa         15          0         0
  #versicolor      0         14         1
  #virginica       0          1        13

混淆矩阵可以看出本来属于某个品种的鸢尾花被分到另一个品种的数量。
5.模型评价（AUC）

roc<-multiclass.roc (as.ordered(testdata$Species) ,as.ordered(iris_pred$iris_pred))
roc
###输出结果为
Call:
multiclass.roc.default(response = as.ordered(testdata$Species),     predictor = as.ordered(iris_pred$iris_pred))

Data: as.ordered(iris_pred$iris_pred) with 3 levels of as.ordered(testdata$Species): setosa, versicolor, virginica.
Multi-class area under the curve: 0.977

Multi-class area under the curve: 0.977，AUC得分为0.977，相当不错。
那roc曲线怎么绘制？pROC只能画二分类问题的曲线，假如你使用其他的数据，并且是二分类问题，那么可使用如下代码

roc<-roc(as.ordered(testdata$Species) ,as.ordered(iris_pred$iris_pred))
plot(roc,print.auc=T, auc.polygon=T, grid=c(0.1, 0.2), grid.col=c("green","red"), 
     max.auc.polygon=T, auc.polygon.col="skyblue",print.thres=T)

交叉验证

# 交叉验证
##先清除NA的数据样本，交叉验证不允许有NA
airquality=na.omit(airquality)
myairquality=cbind(airquality[1:6],matrix(runif(96*nrow(airquality)),nrow(airquality),96))
# 交叉验证中添加随机数的训练集、分组、交叉验证的次数
result=rfcv(myairquality,airquality$Ozone,cv.fold = 5)
# 绘制错误率曲线，观察错误率与使用Markers数量的变化
with(result,plot(n.var,error.cv,log="x",type="o",lwd=2))


##使用replicate多次交叉验证
result=replicate(5,rfcv(myairquality,airquality$Ozone),simplify = FALSE)
error.cv=sapply(result,"[[","error.cv")
matplot(result[[1]]$n.var,cbind(rowMeans(error.cv),error.cv),type = "l",lwd = c(2,rep(1,ncol(error.cv))),col = 1,lty = 1,log="x",xlab = "Number of variables",ylab="CV Error")

多次交叉验证曲线重合，说明结果可信

变量筛选

我用我自己的社会调查数据举例。我的错误率曲线画出来是这样的（代码和上文交叉验证相同）

这说明特征数量不宜超过20左右的某个数字，否则错误率会上升。
那好，我们输出重要性前20个变量，就是我们模型所需要的最佳变量。
注意，这里不是用鸢尾花数据集做的，建议换数据来测试一下。参数可以根据自己的数据调整。

varImpPlot(rf,n.var=20,cex=2)

当然，输出鸢尾花数据集的变量重要性也是可以的

varImpPlot(iris_rf)

换成你的数据，就取前二十个就好了。
那有没有另一种方法筛选变量？
有的。
那就是varSelRF包
1.varSelRF()函数

varSelRF(xdata, Class, c.sd = 1, mtryFactor = 1, ntree = 5000,
         ntreeIterat = 2000, vars.drop.num = NULL, vars.drop.frac = 0.2,
         whole.range = TRUE, recompute.var.imp = FALSE, verbose = FALSE,
         returnFirstForest = TRUE, fitted.rf = NULL, keep.forest = FALSE)

等于单独做一次随机森林。

2.randomVarImpsRF()

randomVarImpsRF(xdata, Class, forest, numrandom = 100,
                whichImp = "impsUnscaled", usingCluster = TRUE,
                TheCluster = NULL)

用现成的随机森林筛选变量，第三个参数放置随机森林对象

3.变量重要性可视化方法
有两种方法
这是基于varSelRF()

rf.vs1 <- varSelRF(x, cl, ntree = 200, ntreeIterat = 100,
                   vars.drop.frac = 0.2)
rf.vs1
plot(rf.vs1)     

这是基于randomVarImpsRF()，可以设置显不显示变量名

randomImportances<-randomVarImpsRF(xdata, Class, forest, numrandom = 100,
                whichImp = "impsUnscaled", usingCluster = TRUE,
                TheCluster = NULL)
randomVarImpsRFplot(randomImportances, forest,
                    whichImp = "impsUnscaled", nvars = NULL,
                    show.var.names = FALSE, vars.highlight = NULL,
                    main = NULL, screeRandom = TRUE,
                    lwdBlack = 1.5,
                    lwdRed = 2,
                    lwdLightblue = 1,
                    cexPoint = 1,
                    overlayTrue = FALSE,
                    xlab = NULL,
                    ylab = NULL)

以我自己社会调查的数据为例，如下，我用的是第二个工具，以现成的随机森林判断指标重要性

数据可视化

随机森林包提供了一个很方便的函数接口
MDSplot(iris_rf,iris$Species,palette = rep((2:6),2),cex=1)
第一个参数是你用的模型名称，后面一个参数意为确定一种特征，根据这种特征标色，在这里当然是根据鸢尾花的种类标色
palette是调色板，要注意rep（）函数的用法。 rep((2:6),2)意思是向量2:6循环2次，而2:6每一个元素都是一种颜色，即2、3、4、5、6分别表示一种颜色。如果颜色不够用就循环使用，循环两次。cex表示点的大小。

MDSplot(iris_rf,iris$Species)

可以看到除了少数不同种类被分到一起，效果还是不错的。

利用随机森林聚类（无监督学习）

set.seed(1234)
rf2=randomForest(iris[,-5])

MDSplot(iris.urf,iris$Species)

可以看出模型中是没有监督项label的，数据可视化后如下。比在监督学习下集中度差了一些。

随机森林简单原理

2020.5.21更新
没想到一拖就是一年哈，久等了各位。
最近考研效率低下，不想学习。。。所以我决定，这是我这篇文章最后一次更新啦~

让我们从种一棵决策树开始

前面提到，随机森林的“随机”二字体现在两个方面，一是变量的随机，二是样本的随机。
我知道这么说很抽象，让我们从最简单开始入手
现在，我们只用一棵决策树(CART)来跑数据，那么其工作原理是什么样的呢？
假设现在对你进行调查，调查的题目是“朋友约你，你会不会出去玩”，目的是要收集朋友每一次约你的情况，假设总共约你约了1000次，我们设置了以下几个变量：降雨、雾霾、气温、活动范围（室内还是室外）。
也就是说，每一次邀约，都是一类样本，而你接不接受邀约，取决于：降雨、雾霾、气温、活动范围

这是一棵决策树决策的全过程，看懂了吗？每个节点里，都有对应的一类样本。比如说图中“雾霾”那个黑色节点，存储了这1000次邀约当中，当时天气没有下雨的那一部分，所以要进一步通过“雾霾”这一变量，进行进一步的细化分类。再比如图中气温在18-30度的那个“出门”节点，它存储的是1000个样本中，当时天气不下雨、没有雾霾、且气温在18-30度的样本。
那么问题来了，决策树自上到下，为什么呈降雨、雾霾、气温、活动范围这一变量排列顺序进行分类？就不能先看活动范围，再看降雨、雾霾，用这种顺序进行分类？
其实是在决策树算法中，我们会先计算熵、或者基尼系数、或者节点纯度，计算公式参考其他文献吧，我就不写了。这些指标都是衡量某个节点当中样本的混乱度的。节点当中样本label越单一，混乱度就越低，节点就越纯。
再用上图举例子，譬如最上面的根节点，包含了所有的可能性“出门、或者“不出门”，所以一定是最混乱的、熵的值是最高的。但是对于决策树第二层“不出门”这个节点，其节点当中全部是“不出门”这一结果的样本，所以其混乱度最低，熵值最小。
那么变量在决策树中的先后顺序，就是为了能使整体决策树的熵，可以最快的下降。成熟的决策树算法（CRAT）会把最优顺序筛选出来，然后使用最优顺序运行。
变量对熵的减小，叫信息增益，可以看作变量的重要性。这也是MeanDecreaseGini的含义——因为有了这个变量，模型的混乱度下降了多少。

随机森林是一片种满了决策树的森林~

好了，我们知道了决策树是怎么工作的，那么随机森林，其实就是N多决策树通过投票原则捆绑在一起的“强分类器”。这里的“强”是相对于单棵决策树来说的。
换句话说，其实这种观点很朴素，有一部分人认为，一颗决策树就那么厉害了，那我多种几棵树，岂不是更牛？
但问题是，我把同样的一棵树，复制1000000次，结果也不会改变啊，我种植那么多一样的树，有什么意义咩？
这就是“随机森林”当中“随机”的来源了，如果我在总样本中随机取样，取总样本的一部分进行训练，那结果势必要发生改变；如果我在所有的变量当中随机取一部分作为进行分类的变量（比如“降雨、雾霾、气温、活动范围”这四个变量我只取其中三个），即使是样本不变，结果也一定会发生变化。
每一棵树，都会对抽取的那一部分样本进行训练，我们把每一棵树判断的样本属于哪一种label的结果，记录下来，取平均值，或者“投票”。比如有100棵树，其中80棵在随机抽样的时候抽到了样本1,这80棵树里有65棵都判断样本1为positive，那么整片森林投票结果就是positive，这也是随机森林这个分类器的分类结果。
衡量变量重要性的MeanDecreaseGini，自然也是各个树在一起取平均值咯.
归根结底，随机，求的就是 变数，有了变化，事情才有了转机。
当然，这只是通俗理解。这么做会不会使分类器性能提升，并不是一个确定的答案。主要的考虑还是避免了单棵决策树的剪枝问题，避免过拟合。

结语

很感慨，终于了却了一桩心事。
我是小白这一事实短期内不会改变。
你们有没有特别不想学习的时候啊，，，感觉就像咸鱼一样，颓废好久了。。。我好怕。。。
大家祝我学习效率高涨哈！