跟着Nature Communications学数据分析：R语言做随机森林模型并对变量重要性排序

用户7010445

发布于 2023-01-06 19:45:27

4K00

代码可运行

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

运行总次数：0

代码可运行

论文

Drivers and trends of global soil microbial carbon over two decades

https://www.nature.com/articles/s41467-022-31833-z#data-availability

这个里面有很多地图的图

数据和代码

https://github.com/gpatoine/drivers_trends_microbial_carbon

这里有随机森林模型然后对变量重要性进行排序的代码，今天的推文我们重复一下论文中的这部分内容，目前能够利用代码和数据运行得到结果，但是还不明白原理和代码中参数的具体作用。今天的内容只是对运行过程的记录。

部分示例数据集截图

image.png

前10个变量是用来构建模型的变量，其中有一个是分类变量，其他都是数值型数据，最后一列Cmic是因变量

读取数据

library(readr)
library(tidyverse)
dat<-read_csv("data/20221215/drivers_trends_microbial_carbon-main/rf_example.csv")
dat %>% head()
dat %>% colnames()

构建随机森林模型

library(caret)
set.seed(202)
predictors<-colnames(dat)[1:10]
model <- train(x = dat[,predictors], 
               y = dat$Cmic,
               method = "rf",
               importance = TRUE,
               tuneGrid = expand.grid(mtry = c(2:4)), # length(predictors) or 2:6
               trControl = trainControl(method = "cv", 
                                        number = 20,
                                        p = 0.75,
                                        savePredictions = TRUE))

这一步需要的时间还是相对比较长的

代码中各个参数都是什么意思还需要仔细看看

输出模型的RSEM和R方

model$results %>% as_tibble %>% filter(mtry == model$bestTune %>% unlist) %>% select(RMSE, Rsquared)

棒棒糖图展示模型重要性

varImp(model)

varImp(model) %>% plot
varImp(model, scale = FALSE) %>% plot

image.png

image.png

还可以用ggplot2画两个柱形图来展示

varImp(model)$importance %>% 
  as.data.frame() %>% 
  rownames_to_column("var") %>% 
  arrange(Overall) %>% 
  mutate(var=factor(var,levels = rev(var))) %>% 
  ggplot(aes(x=var,y=Overall))+
  geom_col(aes(fill=var),show.legend = FALSE)+
  theme_bw()+
  labs(x=NULL) -> p1

varImp(model,scale = FALSE)$importance %>% 
  as.data.frame() %>% 
  rownames_to_column("var") %>% 
  arrange(Overall) %>% 
  mutate(var=factor(var,levels = rev(var))) %>% 
  ggplot(aes(x=var,y=Overall))+
  geom_col(aes(fill=var),show.legend = FALSE)+
  theme_bw()+
  labs(x=NULL) -> p2

library(patchwork)
p1+
  theme(axis.text.x = element_text(angle=60,vjust=1,hjust=1))+
  p2+
  theme(axis.text.x = element_text(angle=60,vjust=1,hjust=1))