随机森林R语言回归学习笔记和一个失败的试验记录

随机森林R语言回归学习笔记

随机森林简介

• 随机森林是一种包含很多决策树（Decision Trees）的集成分类器（Ensemble Classifier）。
• 它输出的类是单个树的类输出的模式(Breiman 2001)。
• 可以处理小n大p问题，高阶相互作用，相关的预测变量等。
• 随机森林可以进行分类或回归分析，得到变量的重要性、模型的误差等指标，并可以进行预测。

随机森林算法发展

• Breiman 2001发表随机森林算法
• 变量的比较（Archer and Kirnes，2008; Groemping，2009）
• 变量间的交互作用( Winham et al．， 2012)

随机森林原理

• 同其他模型一样，随机森林可以解释若干自变量（、、…、）对因变量的作用。
• 如果因变量有个观测值，有个自变量与之相关；
• 在构建分类树的时候，随机森林会随机的在原始数据中重新选择个观测值，其中有的观测值被选择多次，有的没有被选到，这是Bootstrap重新抽样的方法。
• 随机森林随机的从个自变量选择部分变量进行分类树节点的确定，每次构建的分类树可能都不一样；
• 一般情况下，随机森林随机的生成几百个至几千个分类树，然后选择重复程度最高的树作为最终结果。

随机森林优势

1. 对于多数数据集，可以提供较高的分类精度；
2. 可以处理大量的输入变量；
3. 在分类的时候可以估计变量的重要性；
4. 在森林建立的过程中可以生成泛化误差的内部无偏估计；
5. 在大量数据缺失的情况下仍能保持较高的准确性；
6. 它提供了一种检测变量相互作用的实验方法；
7. 它可以平衡类别间数据量不均衡的数据集的误差；
8. 它可以计算案例间的邻近度，对于聚类、边界探测、可视化数据很有用；
9. 使用上述方法，可以扩展到未标记数据，从而实现无监督聚类、异常值检测和数据视图
10. 学习速度快

随机森林劣势

1. 对于某些数据集，随机森林容易过度拟合。这在多噪声的分类/回归任务中更为明显。
2. 随机森林不能把握大量不相关的特征以及减少熵的决策树的集合。
3. 选择随机决策边界比熵减少决策边界更有效，从而使更大的集成更可行。虽然起初这似乎是一个优势，但它具有将计算从训练时间转移到评估时间的效果，这对于大多数应用程序来说实际上是一个缺点。

随机森林的应用

#随机森林回归
rf = randomForest(NDVI~., data = ProvWide1, importance=T, ntree=1000)
varImpPlot(rf)   #绘制自变量排序图

权重排序图

• %IncMSE：percentage of increase of mean square error（Increase in MSE(%)）
• 通过对每一个预测变量随机赋值，如果该预测变量更为重要，那么其值被随机替换后模型预测的误差会增大。
• 和MeanDecreaseAccuracy 是一个概念的
• 该值越大表示该变量的重要性越大
• IncNodePurity: Increase in Node Purity
• 通过残差平方和来度量，代表了每个变量对分类树每个节点上观测值的异质性的影响，从而比较变量的重要性
• 和Mean Decrease Gini一个概念
• 该值越大表示该变量的重要性越大

记录一个失败的试验

最近几天尝试使用NDVI年度数据和气候因子做了一个随机森林的试验。结果不理想，看网上有人说随机森林适合多观测值、多变量的情况，而我这只有20年的数据也就是20个观测值，5个变量，所以可能并不适用。但是这个代码应该还是有一些参考意义的，发出来供大家参考。

贴一个评估论坛上面的说法，供参考

要点总结

• 使用rfPermute包进行了随机森林的计算，得到了%IncMSE、IncNodePurity及对应的p值

%IncMSE、IncNodePurity及对应的p值

library(tidyverse)
library(readxl)
#读取NDVI分省数据
NDVIProv = read.csv("./RFcsv/NDVIprov.csv", header = T)
NDVIProvLong = gather(NDVIProv, key = Year, value = Value, 3:22)
NDVIProvLong$Year = substr(NDVIProvLong$Year, 2, 5)  #去除年份中的“X”
NDVIProvLong$Year = as.numeric(NDVIProvLong$Year)
NDVIProvLong = NDVIProvLong[,2:4]
names(NDVIProvLong) = c("Code", "Year", "NDVI")

#读取风速时间序列
WindSpead = read.csv(file = "./RFcsv/ClimateFactor/1951～2020年中国各省份平均风速年度数据.csv", header=T, encoding="UTF-8")
names(WindSpead) = c("Year", "Code", "Prov", "WindSpead")
#读取气温时间序列
TMP = read_excel(path = "./RFcsv/ClimateFactor/1951～2020年中国各省份平均气温年度数据.xlsx", 
                 sheet = 1, col_names = T)
names(TMP) = c("Year", "Code", "Prov", "TMP")
#读取相对湿度Relative humidity
RH = read_excel(path = "./RFcsv/ClimateFactor/1951～2020年中国各省份相对湿度年度数据.xlsx", 
                sheet = 1, col_names = T)
names(RH) = c("Year", "Code", "Prov", "RH")
#读取日照时数Sunshine Hours
SH = read_excel(path = "./RFcsv/ClimateFactor/1960～2020年中国各省份日照时数年度数据.xlsx", 
                sheet = 1, col_names = T)
names(SH) = c("Year", "Code", "Prov", "SH")
#读取降水量
Pre = read.csv(file = "./RFcsv/ClimateFactor/2000年～2020年中国各省份降水量年度数据.csv", header=T, encoding="UTF-8")
names(Pre) = c("Year", "Code", "Prov", "Pre")

#NDVI和气候因子合并
NDVIProvLong %>% left_join(Pre, by=c("Year", "Code")) %>%
  left_join(WindSpead, by=c("Year", "Code")) %>%
  left_join(TMP, by=c("Year", "Code")) %>%
  left_join(RH, by=c("Year", "Code")) %>%
  left_join(SH, by=c("Year", "Code")) -> NDVICF
NDVICF= NDVICF[,c(1:3,5,7,9,11,13)] 
write.csv(NDVICF, "./RFcsv/NDVICF.csv")

把数据整理成了这个样子，用做随机森林的输入：

NDVICF

后面我先用了一个省的数据进行了试验，然后通过for循环计算多个省的随机森林结果。

#随机森林
library(rfPermute)
CodeUnique = unique(NDVICF$Code)
NDVICF%>% filter(Code==650000) -> NDVICF110000
NDVICF110000 = NDVICF110000[, 3:8]
rf_110000 = rfPermute(NDVI~., NDVICF110000, ntree=80)
plotImportance(rf_110000, plot.type = "bar", ranks = F, imp.type = "IncNodePurity")
rf_110000$rf

improtance_110000 = as.data.frame(importance(rf_110000))

#批量进行随机森林计算，计算每个省的随机森林
for (i in 1:length(CodeUnique)) {
  NDVICF%>%filter(Code == CodeUnique[i]) -> NDVICF110000
  NDVICF110000 = NDVICF110000[, 3:8]
  rf_110000 = rfPermute(NDVI~., NDVICF110000, ntree=80)
  improtance_110000 = as.data.frame(importance(rf_110000))
  write.csv(improtance_110000, paste0("./RFresult/", CodeUnique[i], ".csv"))
}

#读取随机森林计算结果
csvlist = list.files(path = "./RFresult/", pattern = ".csv$")
csvdir = paste0("./RFresult/", csvlist)
csvdata = lapply(csvdir, read.csv)
names(csvdata) = CodeUnique
csvdata2 = as.data.frame(csvdata)

csv1 = read.csv(csvdir[1], header = T)[,1:2]
Prov = rep(CodeUnique[1], 5)
csv1 = cbind(Prov, csv1)
csvdata = csv1

for (i in 2:length(CodeUnique)) {
  csv2 = read.csv(csvdir[i], header = T)[,1:2]
  Prov = rep(CodeUnique[i], 5)
  csv2 = cbind(Prov, csv2)
  csvdata = rbind(csvdata, csv2)
}

names(csvdata) = c("Code", "ClimFactor", "Importance")

#提取省名和代码
ProvCode = unique(RH[,2:3])
csvdata%>%left_join(ProvCode, by="Code") -> csvdata2
write.csv(csvdata2, "./RF/RFimportance.csv")
csvdata2 = read.csv("./RF/RFimportance.csv", header = T)

#以%IncMSE度量重要性
ggplot(data = csvdata2, aes(x=ClimFactor, y=Importance))+
  geom_bar(stat = 'identity')+
  facet_wrap(~Prov, scales = "free_y", ncol = 4)+
  theme_bw()

以%IncMSE度量重要性，这个结果感觉很奇怪，有的是正的，有的是负的，如果有熟悉随机森林的大佬还请留言赐教

上图没有区分通过显著性检验和未通过显著性检验的随机森林结果，下面我换了一种重要性的度量方式，并且加入了显著性检验。

#以IncNodePurity来度量重要性
csvlist = list.files(path = "./RFresult/", pattern = ".csv$")
csvdir = paste0("./RFresult/", csvlist)
csvdata = lapply(csvdir, read.csv)
names(csvdata) = CodeUnique
csvdata2 = as.data.frame(csvdata)

csv1 = read.csv(csvdir[1], header = T)[,c(1, 4:5)]
Prov = rep(CodeUnique[1], 5)
csv1 = cbind(Prov, csv1)
csvdata = csv1

for (i in 2:length(CodeUnique)) {
  csv2 = read.csv(csvdir[i], header = T)[,c(1, 4:5)]
  Prov = rep(CodeUnique[i], 5)
  csv2 = cbind(Prov, csv2)
  csvdata = rbind(csvdata, csv2)
}

#判断显著性，IncNodePurity.pval<0.05则为显著，TRUE
csvdata%>%
  mutate(sigTF = if_else(IncNodePurity.pval<0.05, T, F)) ->csvdata2
csvdata2 = csvdata2[, c(1:3,5)]
#改名
names(csvdata2) = c("Code", "ClimFactor", "Importance", "sig")
#挂接省名
csvdata2%>%left_join(ProvCode, by="Code") -> csvdata2
#重要性绘图
p1 = ggplot(data = csvdata2, aes(x=ClimFactor, y=Importance))+
  geom_bar(stat = 'identity', aes(fill=sig))+
  facet_wrap(~Prov, scales = "free_y", ncol = 4)+
  theme_bw()
p1

write.csv(csvdata2, "./RF/RFimportanceIncNodePurity.csv")

以IncNodePurity来度量重要性并判断显著性，IncNodePurity.pval<0.05则为显著，TRUE（绿色）

看了半天，我实在是没法解释这个图，可能随机森林不适合这种情况吧。

参考文献

1. 李欣海. 随机森林模型在分类与回归分析中的应用[J]. 应用昆虫学报, 2013, 50(04): 1190–1197.
2. Breiman L. Random Forests[J]. Machine Learning, 2001, 45(1): 5–32.
3. http://www.corina.cc/article/28/
4. https://bbs.pinggu.org/thread-3100555-1-1.html