跟着Nature Communications学数据分析：R语言

论文

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))

image.png

后面还有代码是将这个随机森林模型重复运行100次，使用到了map()和map_dfr()函数，这两个函数还得仔细学习一下用法

关于这个代码感兴趣的可以去看看原文提供的代码

示例数据和代码可以给公众号推文点赞，点击在看，最后留言获取

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！