論文

Drivers and trends of global soil microbial carbon over two decades

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

這個(gè)里面有很多地圖的圖

還有自定義圖例形狀的代碼

數(shù)據(jù)和代碼

https://github.com/gpatoine/drivers_trends_microbial_carbon

這里有隨機(jī)森林模型 然后對(duì)變量重要性進(jìn)行排序的代碼，今天的推文我們重復(fù)一下論文中的這部分內(nèi)容，目前能夠利用代碼和數(shù)據(jù)運(yùn)行得到結(jié)果，但是還不明白原理和代碼中參數(shù)的具體作用。今天的內(nèi)容只是對(duì)運(yùn)行過程的記錄。

部分示例數(shù)據(jù)集截圖

image.png

前10個(gè)變量是用來構(gòu)建模型的變量，其中有一個(gè)是分類變量，其他都是數(shù)值型數(shù)據(jù)，最后一列Cmic是因變量

讀取數(shù)據(jù)

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

構(gòu)建隨機(jī)森林模型

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

這一步需要的時(shí)間還是相對(duì)比較長的

代碼中各個(gè)參數(shù)都是什么意思還需要仔細(xì)看看

輸出模型的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畫兩個(gè)柱形圖來展示

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

后面還有代碼是將這個(gè)隨機(jī)森林模型重復(fù)運(yùn)行100次，使用到了map()和map_dfr()函數(shù)，這兩個(gè)函數(shù)還得仔細(xì)學(xué)習(xí)一下用法

關(guān)于這個(gè)代碼感興趣的可以去看看原文提供的代碼

示例數(shù)據(jù)和代碼可以給公眾號(hào)推文點(diǎn)贊，點(diǎn)擊在看，最后留言獲取

歡迎大家關(guān)注我的公眾號(hào)

小明的數(shù)據(jù)分析筆記本

小明的數(shù)據(jù)分析筆記本 公眾號(hào) 主要分享：1、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子；2、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)、基因組學(xué)、群體遺傳學(xué)文獻(xiàn)閱讀筆記；3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！