论文
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
这里有随机森林模型 然后对变量重要性进行排序的代码，今天的推文我们重复一下论文中的这部分内容，目前能够利用代码和数据运行得到结果，但是还不明白原理和代码中参数的具体作用。今天的内容只是对运行过程的记录。
部分示例数据集截图
前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
还可以用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))
后面还有代码是将这个随机森林模型重复运行100次，使用到了map()和map_dfr()函数，这两个函数还得仔细学习一下用法
关于这个代码感兴趣的可以去看看原文提供的代码
示例数据和代码可以给公众号推文点赞，点击在看，最后留言获取
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
