隨機森林工作過程可概括如下：
（1）假設訓練集中共有N個物件、M個變數，從訓練集中隨機有放回地抽取N個物件構建決策樹；
（2）在每一個節點隨機抽取m<M個變數，將其作為分割該節點的候選變數，每一個節點處的變數數應一致；
（3）完整生成所有決策樹，無需剪枝（最小節點為1）；
（4）重複（1）–（3）過程，獲得大量決策樹；終端節點的所屬類別由節點對應的眾數類別決定；
（5）對於新的觀測點，用所有的樹對其進行分類，其類別由多數決定原則生成。
相較於其它分類方法，隨機森林通常具有如下優勢：
分類準確率通常更高；
能夠有效處理具有高維特徵（多元）的資料集，而且不需要降維；
在處理大資料集時也具有優勢；
可應用於具有大量缺失值的資料中；
能夠在分類的同時度量變數對分類的相對重要性。
本篇使用微生物群落研究中的16S擴增子測序資料，展示R包randomForest中的隨機森林方法。
示例資料，R程式碼的百度盤連結：
https://pan.baidu.com/s/10MWBfjBnYIzf6Cx2Zd9CjA
資料集
示例檔案“otu_table.txt”為來自16S測序所獲得的細菌OTUs丰度表格，共計120個樣本，其中60個來自環境c（c組），60個來自環境h（h組）。
接下來使用該資料：
（1）任一OTUs的丰度都很難作為判別兩種不同環境的標準，因此接下來綜合考慮所有OTUs的丰度並進行建模，目的是找到能夠穩定區分兩種環境的代表性OTUs組合（作為生物標誌物）；
（2）透過代表性OTUs的丰度構建預測模型，即僅透過這些OTUs的丰度就能夠判斷樣本分類。
首先讀入資料預處理。
#讀取 OTUs 丰度表
otu <- read.table('otu_table.txt', sep = '\t', row.names = 1, header = TRUE, fill = TRUE)
#過濾低丰度 OTUs 類群，它們對分類貢獻度低，且影響計算效率
#120 個樣本，就按 OTUs 丰度的行和不小於 120 為準吧
otu <- otu[which(rowSums(otu) >= 120), ]
#合併分組，得到能夠被 randomForest 識別計算的格式
group <- read.table('group.txt', sep = '\t', row.names = 1, header = TRUE, fill = TRUE)
otu <- data.frame(t(otu))
otu_group <- cbind(otu, group)
#將總資料集分為訓練集（佔 70%）和測試集（佔 30%）
set.seed(123)
select_train <- sample(120, 120*0.7)
otu_train <- otu_group[select_train, ]
otu_test <- otu_group[-select_train, ]
隨機森林構建
模型擬合
randomForest包方法的細節介紹可參考：
https://www.stat.berkeley.edu/~breiman/RandomForests/
#randomForest 包的隨機森林
library(randomForest)
#隨機森林計算（預設生成 500 棵決策樹），詳情 ?randomForest
set.seed(123)
otu_train.forest <- randomForest(groups ~ ., data = otu_train, importance = TRUE)
otu_train.forest
plot(margin(otu_train.forest, otu_train$groups), main = '觀測值被判斷正確的機率圖')
randomForest()函式從訓練集中有放回地隨機抽取84個觀測點，在每棵樹的每個節點隨機抽取36個變數，從而生成了500棵經典決策樹。
生成樹時沒有用到的樣本點所對應的類別可由生成的樹估計，與其真實類別比較即可得到袋外預測（out-of-bag，OOB）誤差，即OOB estimate of error rate，可用於反映分類器的錯誤率。此處為為1.19%，顯示分類器模型的精準度是很高的，可以有效識別兩類分組。
Confusion matrix比較了預測分類與真實分類的情況，class.error代表了錯誤分類的樣本比例，這裡是很低的：c 組的41個樣本中40個正確分類，h組的43個樣本全部正確分類。
機率圖顯示絕大部分樣本的分類具有非常高的正確率。
若識別模糊，則會出現偏離。
分類器效能測試
不妨使用構建好的分類器分類訓練集樣本，檢視判別的樣本分類情況。
#訓練集自身測試
train_predict <- predict(otu_train.forest, otu_train)
compare_train <- table(train_predict, otu_train$groups)
compare_train
sum(diag(compare_train)/sum(compare_train))
擬合的分類模型返回來重新識別訓練集資料時，甚至糾正了在擬合時的錯誤劃分。
接下來使用測試集資料，進一步評估分類器效能。
#使用測試集評估
test_predict <- predict(otu_train.forest, otu_test)
compare_test <- table(otu_test$groups, test_predict, dnn = c('Actual', 'Predicted'))
compare_test
顯示初步構建的隨機森林分類器能夠準確分類訓練集外的樣本。
尋找代表性的OTUs組合
變數重要性
隨機森林除了分類器外的另一常用功能是識別重要的變數，即計算變數的相對重要程度。
在這裡，就是期望尋找能夠穩定區分兩種環境的代表性OTUs組合（作為生物標誌物）。
##關鍵 OTUs 識別
#查看錶示每個變數（OTUs）重要性的得分
#summary(otu_train.forest)
importance_otu <- otu_train.forest$importance
head(importance_otu)
#或者使用函式 importance()
importance_otu <- data.frame(importance(otu_train.forest))
head(importance_otu)
#可以根據某種重要性的高低排個序，例如根據“Mean Decrease Accuracy”指標
importance_otu <- importance_otu[order(importance_otu$MeanDecreaseAccuracy, decreasing = TRUE), ]
head(importance_otu)
#輸出表格
#write.table(importance_otu, 'importance_otu.txt', sep = '\t', col.names = NA, quote = FALSE)
此處“Mean Decrease Accuracy”和“Mean Decrease Gini”為隨機森林模型中的兩個重要指標。其中，“mean decrease accuracy”表示隨機森林預測準確性的降低程度，該值越大表示該變數的重要性越大；“mean decrease gini”計算每個變數對分類樹每個節點上觀測值的異質性的影響，從而比較變數的重要性。該值越大表示該變數的重要性越大。
到這一步，可從中篩選一些關鍵OTUs作為代表物種，作為有效區分兩種環境的生物標誌物。
該圖展示了其中top30關鍵的OTUs，將它們劃分為“關鍵OTUs”的依據為模型中的兩個重要指標（兩個指標下各自包含30個OTUs，預設由高往低排）。
#作圖展示 top30 重要的 OTUs
varImpPlot(otu_train.forest, n.var = min(30, nrow(otu_train.forest$importance)), main = 'Top 30 – variable importance')
交叉驗證
那麼，選擇多少重要的變數（本示例為OTUs）是更合適的呢？
可根據計算得到的各OUTs重要性的值（如“Mean Decrease Accuracy”），將OTUs由高往低排序後，透過執行重複5次的十折交叉驗證，根據交叉驗證曲線對OTU進行取捨。交叉驗證法的作用就是嘗試利用不同的訓練集/驗證集劃分來對模型做多組不同的訓練/驗證，來應對單獨測試結果過於片面以及訓練資料不足的問題。此處使用訓練集本身進行交叉驗證。
##交叉驗證幫助選擇特定數量的 OTUs
#5 次重複十折交叉驗證
set.seed(123)
otu_train.cv <- replicate(5, rfcv(otu_train[-ncol(otu_train)], otu_train$group, 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))
#擬合線圖
library(ggplot2)
library(splines) #用於在 geom_smooth() 中新增擬合線，或者使用 geom_line() 替代 geom_smooth() 繪製普通折線
p <- ggplot(otu_train.cv, aes(otus, value)) +
geom_smooth(se = FALSE, method = 'glm', formula = y~ns(x, 6)) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) +
labs(title = '',x = 'Number of OTUs', y = 'Cross-validation error')
p
交叉驗證曲線展示了模型誤差與用於擬合的OTUs數量之間的關係。誤差首先會隨OTUs數量的增加而減少，開始時下降非常明顯，但到了特定範圍處，下降幅度將不再有顯著變化，甚至有所增加。再根據簡約性原則，大致選擇重要性排名前30的OTUs就可以了。
#大約提取前 30 個重要的 OTUs
p + geom_vline(xintercept = 30)
#根據 OTUs 重要性排序後選擇，例如根據“Mean Decrease Accuracy”指標
importance_otu <- importance_otu[order(importance_otu$MeanDecreaseAccuracy, decreasing = TRUE), ]
head(importance_otu)
#輸出表格
#write.table(importance_otu[1:30, ], 'importance_otu_top30.txt', sep = '\t', col.names = NA, quote = FALSE)
簡約模型
如上的交叉驗證曲線可反映出，並非所有變數都有助於預測模型，重要性排名靠後的變數的效應不明顯甚至存在著負效應。就本文的示例而言，有些OTUs對於分類的貢獻度並不高，有些可能在組間區別不大甚至會增加錯誤率。
因此，對於一開始構建的隨機森林分類器，很多變數其實是可以剔除的。不妨就以上述選擇的前30個最重要的OTUs代替原資料集中所有的OTUs進行建模，一方面助於簡化分類器模型，另一方面還可提升分類精度。
##簡約分類器
#選擇 top30 重要的 OTUs，例如上述已經根據“Mean Decrease Accuracy”排名獲得
otu_select <- rownames(importance_otu)[1:30]
#資料子集的訓練集和測試集
otu_train_top30 <- otu_train[ ,c(otu_select, 'groups')]
otu_test_top30 <- otu_test[ ,c(otu_select, 'groups')]
#隨機森林計算（預設生成 500 棵決策樹），詳情 ?randomForest
set.seed(123)
otu_train.forest_30 <- randomForest(groups ~ ., data = otu_train_top30, importance = TRUE)
otu_train.forest_30
plot(margin(otu_train.forest_30, otu_test_top30$groups), main = '觀測值被判斷正確的機率圖')
#訓練集自身測試
train_predict <- predict(otu_train.forest_30, otu_train_top30)
compare_train <- table(train_predict, otu_train_top30$groups)
compare_train
#使用測試集評估
test_predict <- predict(otu_train.forest_30, otu_test_top30)
compare_test <- table(otu_test_top30$groups, test_predict, dnn = c('Actual', 'Predicted'))
compare_test
與上文使用所有OTUs構建的分類器相比，OOB estimate of error rate降低，且Confusion matrix中也無錯誤分類（先前是有一個錯誤的），表現為精度提高。
再使用訓練集和測試集評估分類器效能。
#訓練集自身測試
train_predict <- predict(otu_train.forest_30, otu_train_top30)
compare_train <- table(train_predict, otu_train_top30$groups)
compare_train
#使用測試集評估
test_predict <- predict(otu_train.forest_30, otu_test_top30)
compare_test <- table(otu_test_top30$groups, test_predict, dnn = c('Actual', 'Predicted'))
compare_test
即便只根據30個（重要的）OTUs的丰度判斷樣本分類，也是能夠準確劃分的。
##NMDS 排序圖中展示分類
#NMDS 降維
nmds <- vegan::metaMDS(otu, distance = 'bray')
result <- nmds$points
result <- as.data.frame(cbind(result, rownames(result)))
#獲得上述的分類預測結果
predict_group <- c(train_predict, test_predict)
predict_group <- as.character(predict_group[rownames(result)])
#作圖
colnames(result)[1:3] <- c('NMDS1', 'NMDS2', 'samples')
result$NMDS1 <- as.numeric(as.character(result$NMDS1))
result$NMDS2 <- as.numeric(as.character(result$NMDS2))
result$samples <- as.character(result$samples)
result <- cbind(result, predict_group)
head(result)
ggplot(result, aes(NMDS1, NMDS2, color = predict_group)) +
geom_polygon(data = plyr::ddply(result, 'predict_group', function(df) df[chull(df[[1]], df[[2]]), ]), fill = NA) +
geom_point()
