随机森林调参用于训练最佳模型（数据工程）

#5.1随机森林挖掘影响SOC的关键指标#载入工作目录setwd("D:/data")#加载数据所在路径或者通过rstudio左上角工具栏中的seesion→set Working Directory手动载入#载入待使用的所有包library(readxl)#读取excel表格library(doParallel)#多线程运行library(randomForest)#随机森林模型library(foreach)#用于简化并行计算和循环操作library(rfPermute)#计算变量重要性library(psych)#计算相关性library(reshape2)#将宽数据框转换成长数据框library(tidyverse)#整理数据格式library(patchwork)#拼图library(export)#导出为ppt格式library(ggplot2)#绘图#读取并初步处理数据data.raw <- read_xlsx("土壤数据.xlsx",col_names=T,sheet=2)data.raw <- na.omit(data.raw)#去除存在缺失值的样本#将原始数据命名为data.raw，下图为示例数据简要展示，建议正式作图不要在数据中出现中文，可能导出时出现乱码问题，关于此问题后面也给出了解决方案

#逐步建立最优随机森林模型## 将mtry试验数据以列表形式储存，目前还主要是模仿修改，所以看起来有些繁琐，实际逻辑就是将ntree固定为一个较大的值，然后mtry取值为1：特征变量数目，借助for循环探寻##最佳参数，因为其结果文件会包括1：ntree的结果值，所以对于ntree固定为一个较大值即可；另外考虑到可能存在欠拟合导致的高解释率问题，所以提取不同mtry下，解释率前五#的对应参数name <- c("besttest")names = c()for (i in name){  names = c(names,paste(i,seq(1,8,1),sep = "_"))}length(names) #8个组合err = lapply(strsplit(names,"_"),as.numeric) #指定筛选函数，以rsqe为评判标准，rsq值越大，模型解释率越高detect_para = function(err){  set.seed(12345); #保证结果可重复  rf = randomForest( y = data$Soc, x = data[,-c(1,3)],importance = TRUE,#第一列为分组变量，第三列为因变量，因此去除，具体使用时要根据数据情况调整x和y                     mtry = err[[2]],                     ntree = 10000 # 计算error rate时使用,需~23s。  )  res1 = data.frame(mtry =err[[2]] ,                   ntree = order(rf$rsq, decreasing = TRUE)[1],                   rsq = rf$rsq[order(rf$rsq, decreasing = TRUE)[1]]# 解释方差  )   res2 = data.frame(mtry =err[[2]] ,                    ntree = order(rf$rsq, decreasing = TRUE)[2],                    rsq = rf$rsq[order(rf$rsq, decreasing = TRUE)[2]]# 解释方差  )   res3 = data.frame(mtry =err[[2]] ,                    ntree = order(rf$rsq, decreasing = TRUE)[3],                    rsq = rf$rsq[order(rf$rsq, decreasing = TRUE)[3]]# 解释方差  )   res4 = data.frame(mtry =err[[2]] ,                    ntree = order(rf$rsq, decreasing = TRUE)[4],                    rsq = rf$rsq[order(rf$rsq, decreasing = TRUE)[4]]# 解释方差  )   res5 = data.frame(mtry =err[[2]] ,                    ntree = order(rf$rsq, decreasing = TRUE)[5],                    rsq = rf$rsq[order(rf$rsq, decreasing = TRUE)[5]]# 解释方差  )   res <- cbind(res1,res2,res3,res4,res5)  return(res)}#筛选函数建立完毕，接下来针对不同分组展开探究即可
#全部样本最佳随机森林模型data <- data.raw# parallel多核运行函数--加快运行速度cores = detectCores()-8#确定并行运算cpu数，依据电脑配置cl = makeCluster(cores);clusterExport(cl, c("data"));clusterEvalQ(cl, library(randomForest));system.time(  {    results <- parSapply(cl, err, detect_para); # 8 tasks  } )stopCluster(cl)

#这里手动合并，不知道怎么调整结果文件格式，将其调整为三行n列格式，方便查看结果result <-as.data.frame(t(cbind(results[1:3,],results[4:6,],         results[7:9,],results[10:12,],results[13:15,])))#这里是因为list格式文件没办法排序和转置，所以将其转换为数值文件result <- lapply(result, as.numeric)result <- data.frame(result[1],result[2],result[3])result[order(result[,3],decreasing = TRUE)[1:5],]

#去除欠拟合模型，即决策树数量过少时，会存在样本预测值为NA，因此全部样本随机森林模型最优参数为：7,79,rsq=0.8689571#最佳模型建立及不同特征重要性和显著性提取set.seed(12345)rf_all <- randomForest( y = data$Soc, x = data[,-c(1,3)],importance = TRUE,ntree =79,mtry = 7);rf_all

rf_all$rsq#此结果即为当mtry为7时，ntree从1至79时分别对应的模型解释率，这也是为什么不特意对ntree进行调参，因为其结果文件会包括ntree从1至设定值的模型解释率，仅需##将ntree设为较大值即可，我们前面设置为10000

set.seed(12345)all_RF<- rfPermute(y = data$Soc, x = data[,-c(1,3)],importance = TRUE,ntree =79,mtry = 7,nperm=99)importance(all_RF)[,1:2]#最优参数下的变量重要性及显著性

#草地样本最佳随机森林模型data <- data.raw[data.raw$土地利用类型=="草地",]# parallel多核运行函数--加快运行速度cores = detectCores()-8;cl = makeCluster(cores);clusterExport(cl, c("data"));clusterEvalQ(cl, library(randomForest));system.time(  {    results <- parSapply(cl, err, detect_para); # 8 tasks  } )stopCluster(cl)result <-as.data.frame(t(cbind(results[1:3,],results[4:6,],                               results[7:9,],results[10:12,],results[13:15,])))result <- lapply(result, as.numeric)result <- data.frame(result[1],result[2],result[3])result[order(result[,3],decreasing = TRUE)[1:5],]#去除欠拟合模型，最佳4,8,0.5302694#最佳模型建立及不同特征重要性和显著性提取set.seed(12345)rf_grass <- randomForest( y = data$Soc, x = data[,-c(1,3)],importance = TRUE,ntree =8,mtry = 4);rf_grassset.seed(12345)grass_RF<- rfPermute(y = data$Soc, x = data[,-c(1,3)],importance = TRUE,ntree =8,mtry = 4,nperm=99)importance(grass_RF)[,1:2]
#耕地样本最佳随机森林模型data <- data.raw[data.raw$土地利用类型=="耕地",]# parallel多核运行函数--加快运行速度cores = detectCores()-8;cl = makeCluster(cores);clusterExport(cl, c("data"));clusterEvalQ(cl, library(randomForest));system.time(  {    results <- parSapply(cl, err, detect_para); # 8 tasks  } )stopCluster(cl)result <-as.data.frame(t(cbind(results[1:3,],results[4:6,],                               results[7:9,],results[10:12,],results[13:15,])))result <- lapply(result, as.numeric)result <- data.frame(result[1],result[2],result[3])result[order(result[,3],decreasing = TRUE)[1:5],]#去除欠拟合模型，最佳8,81,0.8172151#最佳模型建立及不同特征重要性和显著性提取set.seed(12345)rf_land <- randomForest( y = data$Soc, x = data[,-c(1,3)],importance = TRUE,ntree =81,mtry = 8);rf_landset.seed(12345)land_RF<- rfPermute(y = data$Soc, x = data[,-c(1,3)],importance = TRUE,ntree =81,mtry = 8,nperm=99)importance(land_RF)[,1:2]
#林地样本最佳随机森林模型参数调整data <- data.raw[data.raw$土地利用类型=="林地",]# parallel多核运行函数--加快运行速度cores = detectCores()-8;cl = makeCluster(cores);clusterExport(cl, c("data"));clusterEvalQ(cl, library(randomForest));system.time(  {results <- parSapply(cl, err, detect_para)} )stopCluster(cl)result <-as.data.frame(t(cbind(results[1:3,],results[4:6,],                               results[7:9,],results[10:12,],results[13:15,])))result <- lapply(result, as.numeric)result <- data.frame(result[1],result[2],result[3])result[order(result[,3],decreasing = TRUE)[1:5],]#去除欠拟合模型，最佳6,6,0.01731922#最佳模型建立及不同特征重要性和显著性提取set.seed(12345)rf_forest <- randomForest( y = data$Soc, x = data[,-c(1,3)],importance = TRUE,ntree =6,mtry = 6);rf_forestset.seed(12345)forest_RF<- rfPermute(y = data$Soc, x = data[,-c(1,3)],importance = TRUE,ntree =6,mtry = 6,nperm=99)importance(forest_RF)[,1:2]

#整理最佳模型解释率bar<- data.frame(variable = c( "整体"  , "草地", "耕地" ,"林地"),                 Exp = c(86.9,53.03,81.72,1.73))#两列数据，一列分组，一列最优模型解释率bar$variable<- factor(bar$variable, levels = bar$variable)#指定分组排列顺序，实为可视化时x轴的展示顺序，可针对其 levels展开修改barplot<- ggplot(bar, aes(variable, Exp))+  geom_bar(stat = "identity", width=0.5,fill = "steelblue")+  scale_y_continuous(expand = c(0,0),limits = c(0,90),breaks = seq(0,80,20))+  scale_x_discrete(expand = c(0.17,0))+  theme_classic(base_line_size = 0.75)+  theme(panel.background = element_blank(),        axis.ticks.x = element_blank(),        axis.title.x = element_blank(),        axis.text.x = element_blank(),        axis.title.y = element_blank(),        axis.text.y = element_text(color = "black", size = 10))+  labs(title = "Explained variation (%)", size = 12,face="bold")barplot

#相关性和重要性数据整理data_cor <- rbind(data.raw,data.raw)data_cor$group <- c(rep("整体",nrow(data.raw)),data.raw$土地利用类型)data_cor$group <- factor(data_cor$group,levels = c("整体","草地","耕地","林地"))r <- data.frame(env=colnames(data_cor)[-c(1,3,11)],                整体=cor(data_cor[data_cor$group=="整体",2:10],method ="spearman")[2,-2],                草地=cor(data_cor[data_cor$group=="草地",2:10],method ="spearman")[2,-2],                耕地=cor(data_cor[data_cor$group=="耕地",2:10],method ="spearman")[2,-2],                林地=cor(data_cor[data_cor$group=="林地",2:10],method ="spearman")[2,-2])r_plot <- melt(r,id = "env", value.name = "correlation")#整理结果以便作图r_plot$variable<- factor(r_plot$variable, levels = r_plot$variable[c(1,9,17,25)])r_plot$env <- factor(r_plot$env,levels = c( "土壤剖面深度", "Soc","N(μg/g)", "P(μg/g)",                                         "TFe2O3", "SM","Fe","容重", "田间持水率" ))circle<- data.frame(env=colnames(data_cor)[-c(1,3,11)]) %>%  left_join(data.frame(env = row.names(importance(all_RF)),                       整体 = ifelse(importance(all_RF)[,2]<0.05,                                    importance(all_RF)[,1], NA))) %>%  left_join(data.frame(env = row.names(importance(grass_RF)),                       草地 = ifelse(importance(grass_RF)[,2]<0.05,                                         importance(grass_RF)[,1], NA))) %>%  left_join(data.frame(env = row.names(importance(land_RF)),                       耕地 = ifelse(importance(land_RF)[,2]<0.05,                                             importance(land_RF)[,1], NA))) %>%  left_join(data.frame(env = row.names(importance(forest_RF)),                       林地 = ifelse(importance(forest_RF)[,2]<0.05,                                     importance(forest_RF)[,1], NA)))%>%melt(id = "env", value.name = "Importance") ;circlecircle$env<- factor(circle$env, levels=rev(c( "土壤剖面深度", "Soc","N(μg/g)", "P(μg/g)",                                  "TFe2O3", "SM","Fe","容重", "田间持水率" )))#是否rev不重要circle#导出结果write.table(circle, file ="D:/result/5/5.1随机森林特征重要性.txt" , sep = "\t", quote = F)write.table(r_plot, file ="D:/result/5/5.1随机森林环境因子相关性.txt" , sep = "\t", quote = F)write.table(bar, file ="D:/result/5/5.1随机森林模型解释率.txt" , sep = "\t", quote = F)#r_plot <- read.table("5.1随机森林环境因子相关性.txt",sep = "\t",header = T,row.names=1)#circle <- read.table("5.1随机森林特征重要性.txt",sep = "\t",header = T,row.names=1)##绘图时遇中文无法显示问题，且中文过多时无法出图，解决方法如下install.packages("sysfonts")install.packages("showtextdb")install.packages("showtext")library(sysfonts)library(showtextdb)library(showtext)showtext_auto()###热图加重要性图绘制heatmap<-  ggplot()+geom_tile(data = r_plot, aes(variable,env,fill = correlation))+ #热图  scale_fill_gradientn(colours = rev(c("#e0745a", "#fbfbf7", "#63a1cb")),  #自定义颜色                       limits = c(-1, 1))+  geom_point(data = circle, aes(x = variable, y = env,                                 size = Importance), shape = 21)+  scale_size_area(breaks = seq(2, 17, 5),                  limits = c(1, 18), max_size = 8)+#圆圈大小  theme(panel.grid = element_blank(), panel.background = element_rect(color = "black",fill=NA), legend.key = element_blank(),         axis.text.x = element_text(size = 10,colour = "black", angle =45, hjust = 1, vjust = 1),         axis.text.y = element_text(size = 10,colour = "black"), axis.ticks = element_line(color = 'black'),        legend.title = element_text(size = 12,face="bold"),        legend.text =element_text(size = 10,colour = "black"))+  scale_x_discrete(expand = c(0.17, 0)) +  scale_y_discrete(expand = c(0.078, 0)) +  guides(size=guide_legend(order = 1),fill=guide_colorbar(order=2))+  labs(y = '', x = '', fill = 'Correlation',size="Importance(%)")heatmap

plot <- barplot + heatmap +plot_layout(ncol = 1, heights = c(1, 4))#拼图，一列，前者在上面，后者在下面，高度比为1:4

#保存结果文件，分别为pdf，pptx格式graph2ppt(file = "D:/result/5/5.1随机森林出图.ppt",           height = 13/2.54, width = 8.3/2.54, append = T)ggsave("D:/result/5/5.1随机森林出图.pdf",       plot,width = 12,height = 14.6,units = "cm")