R語言 | randomForest包的隨機森林回歸模型以及對重要變量的選擇

R包randomForest的隨機森林回歸模型以及對重要變量的選擇關於隨機森林（random forest），前文「隨機森林分類以及對重要變量的選擇」中已經對其基本原理作了簡單概括。在前文中，響應變量是一組類別變量（代表了樣本的分組信息），此時隨機森林用於執行監督分類的功能，模型的精度在於通過選定的預測變量組合將樣本正確劃分歸類的概率。
而當響應變量是一組連續變量時，可將回歸森林用於回歸。此時，回歸模型的精度體現了選定的預測變量對響應變量數值改變有關方差的整體解釋率（可以理解為回歸的R2）。

  

首先來看一篇與植物根系微生物組有關研究，將文中涉及到隨機森林回歸的部分內容節選出來以幫助大家理解。

Edwards等（2018）分析了來自兩個地區（加利福尼亞州、阿肯色州）以及三個連續種植季（2014、2015、2016）的水稻根系微生物組，意在闡述水稻根系微生物群落的時間動態特徵，並解釋驅動微生物群落組建的成因。文中部分內容使用到隨機森林的回歸模型建立微生物豐度和水稻生長時期的關係，用以識別重要的時間響應微生物類群。

 

作者首先分析了水稻根系微生物群落的組成結構，顯示其受到地理區域（土壤環境）、植物生長時期以及植物根系不同部位的多重影響。隨後分析各個因素的相對效應，發現水稻生長過程中能夠逐漸對其根系菌群進行選擇，植物生長時期的效應越來越明顯，根系菌群朝向一組保守的物種組成趨勢發展，同時地理環境效應的比重逐漸減弱。主要體現在：

（1）不同地理區域種植的水稻，隨著水稻生長時期的推移，根際和根內菌群的相似度逐漸增加，並在後期趨於穩定，與種植季節無關：

（2）對於水稻根際和根內菌群，地區間獨特OTU的相對豐度均隨時間推移而不斷下降，剩餘的OTU中，兩地區間共有OTU的比重上升。

圖A-C，水稻根系微生物群落的PCoA，分別突出了群落組成相似度與根系區域（A）、地理區域（B）或植物生長時間（C）的關係。

圖D，根區、根際和根內樣本中，微生物群落相似度隨植物生長的關係。

 

總之存在一類具有顯著時間特徵的微生物類群，它們的豐度變化與水稻的生長時期密切相關，無論是增加（根系富集）或減少（根系排斥）。隨後，作者期望尋找可用於區分水稻生長時期的微生物OTU組合，包括根際和根內菌群，就使用了隨機森林。

作者選擇了一部分的樣本（加利福尼亞州2014年和阿肯色州2016年的部分數據）用作訓練集，通過隨機森林分別將根際和根內所有微生物OTU的相對豐度與水稻生長時期進行了回歸，並根據各個OTU對模型預測準確性的相對貢獻按重要性進行排名。十折交叉驗證用於評估模型中OTU數量與模型誤差的關係，並按重要性排名對OTU進行選擇，刪除了大部分不重要OTU後，最終保留的OTU用於構建最終的隨機森林回歸模型。

結果中，共確定了85個最重要的與水稻生長時期密切關聯的OTU。訓練集的數據顯示，通過這85個OTU構建的隨機森林回歸模型分別解釋了根際和根內菌群對植物年齡有關方差的91.5％和88.4％，表明它們的豐度組合能夠準確預測水稻的生長時期。此外，測試集以及驗證集數據（涉及了訓練集樣本外的其它採樣時間、地理區域的樣本）也一致地表明，通過這85個OTU構建的隨機森林回歸模型同樣準確預測了水稻的生長時期，凸顯了水稻根際和根內存在特定的微生物類群，在跨地理區域和季節的水稻植物的生命周期中表現出一致的模式。

圖A，橫坐標是水稻的實際生長時期，縱坐標是隨機森林模型通過根際或根內微生物豐度預測的水稻生長時期，模型具有很高的精度。

圖B，重要的85個OTU在水稻根際或根內的豐度與時間的關係。

圖C，85個OTU根據豐度的時間響應模式，可分為早期定植（豐度早期很高，但隨時間降低）或晚期定植（豐度早期很低，但隨時間增加）以及複合型（早期晚期定植區分不明顯）類群，柱形圖展示了3種類群中OTU平均豐度變化趨勢。

 

文中後續還有內容涉及到隨機森林模型的應用，這裡就不再展示了。節選的部分章節想必也足以能夠幫助大家了解這個方法的實際應用了。

此外，該文獻的補充材料中提供了有關OTU豐度表、分析過程的R代碼等非常全面的數據，大家若有興趣參考原文即可。

接下來，展示隨機森林回歸及對重要變量選擇在R語言中的實現方法。

通過R包randomForest的隨機森林執行回歸

  

對於隨機森林應用在類別型響應變量中的分類功能，前文「隨機森林分類模型以及對重要變量的選擇」中，已經以R包randomForest中的方法為例展示了如何通過隨機森林構建分類預測模型（分類模型的訓練和測試），以及篩選對區分已知分類具有高度鑑別模式的變量組合（評估變量的相對重要性）。

在下文中，將響應變量更換為連續型變量，繼續展示R包randomForest的隨機森林回歸方法以及實現對重要變量的選擇。其實無論執行的分類還是回歸，過程中使用的方法都是一樣的，R函數還是那些R函數，區別只在於對部分函數的參數設置以及結果的解讀方式上。

下文的測試數據，R代碼等的百度盤連結（提取碼，fljz）：

https://pan.baidu.com/s/1Lfk21hGrWDehodWchBiIIg

若百度盤失效，也可在GitHub的備份中獲取：

https://github.com/lyao222lll/sheng-xin-xiao-bai-yu

 

植物根系菌群結構與植物生長密切相關，到目前為止，已經在許多研究中都有報導了，這已經是個共識。因此，下文的示例數據也同樣來自某植物根際區域細菌群落組成的16S擴增子測序數據，類似地，仿照上文文獻中的過程，通過隨機森林建立微生物與植物生長時期的響應關係，根據微生物豐度推斷植物生長時期，並尋找一組重要的時間特徵類群。

示例數據「otu_table.txt」中，共記錄了45個連續生長時間中植物根際土壤樣本中細菌OTU的相對豐度信息。

「plant_age.txt」中，記錄了這45個根際土壤樣本對應的植物生長時間，時間單位是天。

   

將OTU豐度數據讀入R中，可以事先作一些預處理，例如剔除低豐度類群等。

隨後，對應OTU豐度樣本與植物生長時間的關係，加載randomForest包並運行隨機森林。

##數據預處理
#讀取 OTUs 豐度表
otu <- read.delim('otu_table.txt', row.names = 1)

#過濾低豐度 OTUs 類群，它們對分類貢獻度低，且影響計算效率
#例如剔除總豐度低於 0.05% 的值
otu <- otu[which(rowSums(otu) >= 0.0005), ]

#合併有關於植物生長時間的信息
plant <- read.delim('plant_age.txt', row.names = 1)

otu <- data.frame(t(otu))
otu <- otu[rownames(plant), ]
otu <- cbind(otu, plant)

#為了方便後續評估隨機森林模型的性能
#將總數據集分為訓練集（佔 70%）和測試集（佔 30%）
set.seed(123)
train <- sample(nrow(otu), nrow(otu)*0.7)
otu_train <- otu[train, ]
otu_test <- otu[-train, ]

##randomForest 包的隨機森林
library(randomForest)

#隨機森林計算（默認生成 500 棵決策樹），詳情 ?randomForest
set.seed(123)
otu_train.forest <- randomForest(plant_age~., data = otu_train, importance = TRUE)
otu_train.forest

結果中，% Var explained體現了預測變量（用於回歸的所有OTU）對響應變量（植物年齡）有關方差的整體解釋率。在本示例中，剔除了低豐度的OTU後，剩餘的OTU（約2600個左右）解釋了約89.17%的總方差，可以理解為該回歸的R2=0.8917，相當可觀的一個數值，表明該植物根際細菌的群落結構隨植物生長密切相關。

 

查看該模型的預測性能，可以看到具有較高的精度。

#使用訓練集，查看預測精度
plant_predict <- predict(otu_train.forest, otu_train)

plot(otu_train$plant_age, plant_predict, main = '訓練集',
xlab = 'Plant age (days)', ylab = 'Predict')
abline(1, 1)

#使用測試集，評估預測性能
plant_predict <- predict(otu_train.forest, otu_test)

plot(otu_test$plant_age, plant_predict, main = '測試集',
xlab = 'Plant age (days)', ylab = 'Predict')
abline(1, 1)

   

但是，並非這所有的2600餘個OTU都對回歸的精度具有可觀的貢獻。有些OTU豐度的時間特徵並不明顯，可能在回歸中產生較大噪聲，對模型精度帶來較高的誤差。因此，最好將低貢獻的OTU去除。

基於已經構建好的隨機森林回歸模型，可以從中評估OTU的重要性。將OTU按重要程度高低進行排名後，選擇排名靠前的一部分OTU，這些重要的OTU是明顯的與植物生長時間密切關聯的一些細菌類群。

##OTU 的重要性評估
#查看表示每個預測變量（細菌 OTU）重要性的得分
#summary(otu_train.forest)
importance_otu <- otu_train.forest$importance
head(importance_otu)

#或者使用函數 importance()
importance_otu <- data.frame(importance(otu_train.forest), check.names = FALSE)
head(importance_otu)

#作圖展示 top30 重要的 OTUs
varImpPlot(otu_train.forest, n.var = min(30, nrow(otu_train.forest$importance)),
main = 'Top 30 - variable importance')

「%IncMSE」即increase in mean squared error，通過對每一個預測變量隨機賦值，如果該預測變量更為重要，那麼其值被隨機替換後模型預測的誤差會增大。因此，該值越大表示該變量的重要性越大；

「IncNodePurity」即increase in node purity，通過殘差平方和來度量，代表了每個變量對分類樹每個節點上觀測值的異質性的影響，從而比較變量的重要性。該值越大表示該變量的重要性越大。

對於「%IncMSE」或「IncNodePurity」，二選一作為判斷預測變量重要性的指標。需注意的是，二者的排名存在一定的差異。

#可以根據某種重要性的高低排個序，例如根據「IncNodePurity」指標
importance_otu <- importance_otu[order(importance_otu$IncNodePurity, decreasing = TRUE), ]
head(importance_otu)

#輸出表格
#write.table(importance_otu, 'importance_otu.txt', sep = '\t', col.names = NA, quote = FALSE)

 

交叉驗證

那麼，最終選擇多少重要的預測變量（本示例為OTUs）是更合適的呢？

可通過執行十折交叉驗證，根據交叉驗證曲線對OTU進行取捨。交叉驗證法的作用就是嘗試利用不同的訓練集/驗證集劃分來對模型做多組不同的訓練/驗證，來應對單獨測試結果過於片面以及訓練數據不足的問題。此處使用訓練集本身進行交叉驗證。

##交叉驗證輔助評估選擇特定數量的 OTU
#5 次重複十折交叉驗證
set.seed(123)
otu_train.cv <- replicate(5, rfcv(otu_train[-ncol(otu_train)], otu_train$plant_age, cv.fold = 10, step = 1.5), simplify = FALSE)
otu_train.cv

#提取驗證結果繪圖
otu_train.cv <- data.frame(sapply(otu_train.cv, '[[', 'error.cv'))
otu_train.cv$otus <- rownames(otu_train.cv)
otu_train.cv <- reshape2::melt(otu_train.cv, id = 'otus')
otu_train.cv$otus <- as.numeric(as.character(otu_train.cv$otus))

otu_train.cv.mean <- aggregate(otu_train.cv$value, by = list(otu_train.cv$otus), FUN = mean)
head(otu_train.cv.mean, 10)

#擬合線圖
library(ggplot2)

ggplot(otu_train.cv.mean, aes(Group.1, x)) +
geom_line() +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) +
labs(title = '',x = 'Number of OTUs', y = 'Cross-validation error')

交叉驗證曲線展示了模型誤差與用於擬合的OTU數量之間的關係。誤差首先會隨OTUs數量的增加而減少，開始時下降非常明顯，但到了特定範圍處，下降幅度將不再有顯著變化，甚至有所增加。

根據交叉驗證曲線，提示保留9-13個重要的OTU即可以獲得理想的回歸結果，因為此時的誤差達到最小。

 

因此，根據計算得到的各OUT重要性的值（如「IncNodePurity」），將OTU由高往低排序後，最後大約選擇前9-13個OTU就可以了。

#提示保留 9-13 個重要的 OTU，可以使隨機森林回歸的精度最大化
#首先根據某種重要性的高低排個序，例如根據「IncNodePurity」指標
importance_otu <- importance_otu[order(importance_otu$IncNodePurity, decreasing = TRUE), ]

#然後取出排名靠前的 OTU，例如 top10 最重要的 OTU
importance_otu.select <- importance_otu[1:10, ]
importance_otu.select

#輸出表格
#write.table(importance_otu.select, 'importance_otu.select.txt', sep = '\t', col.names = NA, quote = FALSE)

#有關 OTU 的物種分類信息等其它細節，可後續自行補充

#不妨簡單查看下這些重要的 OTU 豐度與植物生長時間的關係
#可以看到趨勢非常明顯，包括根際富集或排斥等都有涉及
otu_id.select <- rownames(importance_otu.select)
otu.select <- otu[ ,c(otu_id.select, 'plant_age')]
otu.select <- reshape2::melt(otu.select, id = 'plant_age')

ggplot(otu.select, aes(x = plant_age, y = value)) +
geom_point() +
geom_smooth() +
facet_wrap(~variable, ncol = 3, scale = 'free_y') +
labs(title = '',x = 'Plant age (days)', y = 'Relative abundance')

   

現在來看，只使用選擇的10個重要的OTU構建回歸模型，效果如何。

##只包含 10 個重要預測變量的簡約回歸
otu.select <- otu[ ,c(otu_id.select, 'plant_age')]

#為了方便後續評估隨機森林模型的性能，將總數據集分為訓練集（佔 70%）和測試集（佔 30%）
set.seed(123)
train <- sample(nrow(otu.select), nrow(otu.select)*0.7)
otu_train.select <- otu.select[train, ]
otu_test.select <- otu.select[-train, ]

#隨機森林計算（默認生成 500 棵決策樹），詳情 ?randomForest
set.seed(123)
otu_train.select.forest <- randomForest(plant_age~., data = otu_train.select, importance = TRUE)
otu_train.select.forest

#使用訓練集，查看預測精度
plant_predict <- predict(otu_train.select.forest, otu_train.select)

plot(otu_train.select$plant_age, plant_predict, main = '訓練集',
xlab = 'Plant age (days)', ylab = 'Predict')
abline(1, 1)

#使用測試集，評估預測性能
plant_predict <- predict(otu_train.select.forest, otu_test.select)

plot(otu_test.select$plant_age, plant_predict, main = '測試集',
xlab = 'Plant age (days)', ylab = 'Predict')
abline(1, 1)

結果顯示，與先前使用全部的OTU（去除低豐度後約2600多個）相比，只使用這10個更重要OTU的隨機森林回歸模型，獲得了更出色的效果。

一是體現在% Var explained，模型中預測變量（這10個重要的OTU）對響應變量（植物年齡）有關方差的整體解釋率達到了96.49%，大於先前的89.17%，這得益於排除了不重要或者高噪聲的OTU。二是預測性能更佳，特別是對於測試集，與先前的分布圖相比較（見上文），植物實際生長時間與通過根際菌群預測的植物年齡更趨一致。

     

  

上文過程中，為了評估隨機森林回歸模型的預測性能，將數據拆分為訓練集和測試集分開討論。但如果目的不涉及預測，只是為了尋找重要的預測變量，完全可以直接使用所有數據代入回歸，無需再區分訓練集與測試集，這也是很多文獻中經常提到的方法。

例如另一篇類似的研究，Zhang等（2018）尋找水稻根系中與植物生長期密切相關的biomarkers時，同樣使用的隨機森林回歸。但區別是文中沒有涉及預測，關注的重點在於識別重要的微生物類群。

    

  

Edwards J, Santosmedellin C, Liechty Z, et al. Compositional shifts in root-associated bacterial and archaeal microbiota track the plant life cycle in field-grown rice. PLOS Biology, 2018, 16(2).Zhang J, Zhang N, Liu Y, et al. Root microbiota shift in rice correlates with resident time in the field and developmental stage. Science China-life Sciences, 2018, 61(6): 613-621.