本结内容属于思农数据工程产物,欢迎大家批评指正。
本次想要复现的图片是焦硕老师于2020年发表在Global Change Biology上的Abundant fungi adapt tobroaderenvironmental gradients than rare fungi in agricultural fields中的figure1b,如下图所示,其主要用来研究不同环境因子与丰富和稀有物种的相对丰度的关系。此图有众多复现,但是大多数不涉及到随机森林模型的调参,这与标题中的best random forest model有一定出入,因此本次想分享下自己学习过程中所认知的bestrandomforest model,并结合手头上的土壤示例数据进行还原。
#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,下图为示例数据简要展示,建议正式作图不要在数据中出现中文,可能导出时出现乱码问题,关于此问题后面也给出了解决方案
这里主要探寻影响土壤SOC的关键指标,分组变量为土地利用类型,因变量为Soc,剩余环境因子为自变量,这里与例图基本一致,只是把因变量微生物相对丰度换成了土壤SOC,如果想做例图一样的研究,仅需把soc列改为某一类微生物在不同样本中的相对丰度,剩余列为对应样本的环境变量即可 。我理解的是随机森林做回归,因变量只能是一列数据,关于beta多样性等两两之间的指标,想要借助于随机森林模型探寻其影响因素,我的想法是因变量为关注指标,自变量就是对应两两样本之间环境变量差值的绝对值,仅供参考。
紧接着我们就通过for循环确定示例数据的最佳随机森林模型,我目前的理解中影响模型解释率比较重要的就是两个指标,一个是mtry,即每棵树在进行节点分裂时所考虑的变量个数;另外一个就是ntree,即决策树的数量,过少会欠拟合,过多会过拟合。在随机森林模型中,决策树数量默认为500,回归mtry默认为自变量个数的三分之一,分类默认为自变量个数的平方根。最后就是关于随机森林是否应该区分训练集和测试集,付费得来的知识告诉我随机森林分析使用2/3的样本用于构建决策树,1/3样本作为袋外数据(out-of-bag,OOB)进行模型验证。在树集足够大的情况下,OOB方法等同于cross-validation,随机森林一般不用分训练集和测试集进行分析,那这里我们就不再区分训练集和测试集了,直接使用全部样本。
#逐步建立最优随机森林模型
## 将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]#最优参数下的变量重要性及显著性
接下来的代码就是依次按照上面的步骤分别探寻不同草地类型,其Soc的重要影响因子,不再一一解释
#草地样本最佳随机森林模型
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_grass
set.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_land
set.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_forest
set.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]
接下来针对上述结果进行可视化,重现例图,主要参考Datasphere\:R语言绘制随机森林重要性+相关性热图https://mp.weixin.qq.com/s/nnsH3OqOjLY5KKYh93_FGA
#整理最佳模型解释率
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") ;circle
circle$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")
主要借用于焦老师的分析思路,对土壤soc的影响因素展开了探究,最后的结果图也与目标还原图基本类似,本次分享结束。
作者:思农生信团队