一套完整的基于随机森林的机器学习流程（特征选择、交叉验证、模型评估））...
機器學習實操(以隨機森林為例)
為了展示隨機森林的操作，我們用一套早期的前列腺癌和癌旁基因表達芯片數據集，包含102個樣品(50個正常，52個腫瘤)，2個分組和9021個變量 (基因)。（https://file.biolab.si/biolab/supp/bi-cancer/projections/info/prostata.html）
數據格式和讀入數據
輸入數據為標準化之后的表達矩陣，基因在行，樣本在列。隨機森林對數值分布沒有假設。每個基因表達值用于分類時是基因內部在不同樣品直接比較，只要是樣品之間標準化的數據即可，其他任何線性轉換如log2，scale等都沒有影響 (數據在：https://gitee.com/ct5869/shengxin-baodian/tree/master/machinelearning)。
樣品表達數據：
prostat.expr.txt
樣品分組信息：
prostat.metadata.txt
基因表達表示例如下：
expr_mat[1:4,1:5]## normal_1 normal_2 normal_3 normal_4 normal_5 ## AADAC 1.3 -1 -7 -4 5 ## AAK1 0.4 0 10 11 8 ## AAMP -0.4 20 -7 -14 12 ## AANAT 143.3 19 397 245 328Metadata表示例如下
head(metadata)## class ## normal_1 normal ## normal_2 normal ## normal_3 normal ## normal_4 normal ## normal_5 normal ## normal_6 normaltable(metadata)## metadata ## normal tumor ## 50 52樣品篩選和排序
對讀入的表達數據進行轉置。通常我們是一行一個基因，一列一個樣品。在構建模型時，數據通常是反過來的，一列一個基因，一行一個樣品。每一列代表一個變量 (variable)，每一行代表一個案例 (case)。這樣更方便提取每個變量，且易于把模型中的x,y放到一個矩陣中。
樣本表和表達表中的樣本順序對齊一致也是需要確保的一個操作。
# 表達數據轉置 # 習慣上我們是一行一個基因，一列一個樣品 # 做機器學習時，大部分數據都是反過來的，一列一個基因，一行一個樣品 # 每一列代表一個變量 expr_mat <- t(expr_mat) expr_mat_sampleL <- rownames(expr_mat) metadata_sampleL <- rownames(metadata)common_sampleL <- intersect(expr_mat_sampleL, metadata_sampleL)# 保證表達表樣品與METAdata樣品順序和數目完全一致 expr_mat <- expr_mat[common_sampleL,,drop=F] metadata <- metadata[common_sampleL,,drop=F]判斷是分類還是回歸
前面讀數據時已經給定了參數stringsAsFactors =T，這一步可以忽略了。
如果group對應的列為數字，轉換為數值型 - 做回歸
如果group對應的列為分組，轉換為因子型 - 做分類
隨機森林一般分析
library(randomForest)# 查看參數是個好習慣 # 有了前面的基礎概述，再看每個參數的含義就明確了很多 # 也知道該怎么調了 # 每個人要解決的問題不同，通常不是別人用什么參數，自己就跟著用什么參數 # 尤其是到下游分析時 # ?randomForest# 查看源碼 # randomForest:::randomForest.default加載包之后，直接分析一下，看到結果再調參。
# 設置隨機數種子，具體含義見 https://mp.weixin.qq.com/s/6plxo-E8qCdlzCgN8E90zg set.seed(304)# 直接使用默認參數 rf <- randomForest(expr_mat, metadata[[group]])查看下初步結果, 隨機森林類型判斷為分類，構建了500棵樹，每次決策時從隨機選擇的94個基因中做最優決策 (mtry)，OOB估計的錯誤率是9.8%，挺高的。
分類效果評估矩陣Confusion matrix,顯示normal組的分類錯誤率為0.06，tumor組的分類錯誤率為0.13。
rf## ## Call: ## randomForest(x = expr_mat, y = metadata[[group]]) ## Type of random forest: classification ## Number of trees: 500 ## No. of variables tried at each split: 94 ## ## OOB estimate of error rate: 9.8% ## Confusion matrix: ## normal tumor class.error ## normal 47 3 0.0600000 ## tumor 7 45 0.1346154隨機森林標準操作流程 (適用于其他機器學習模型)
拆分訓練集和測試集
library(caret) seed <- 1 set.seed(seed) train_index <- createDataPartition(metadata[[group]], p=0.75, list=F) train_data <- expr_mat[train_index,] train_data_group <- metadata[[group]][train_index]test_data <- expr_mat[-train_index,] test_data_group <- metadata[[group]][-train_index]dim(train_data)## [1] 77 9021dim(test_data)## [1] 25 9021Boruta特征選擇鑒定關鍵分類變量
# install.packages("Boruta") library(Boruta) set.seed(1)boruta <- Boruta(x=train_data, y=train_data_group, pValue=0.01, mcAdj=T, maxRuns=300)boruta## Boruta performed 299 iterations in 1.937513 mins. ## 46 attributes confirmed important: ADH5, AGR2, AKR1B1, ANGPT1, ## ANXA2.....ANXA2P1.....ANXA2P3 and 41 more; ## 8943 attributes confirmed unimportant: AADAC, AAK1, AAMP, AANAT, AARS ## and 8938 more; ## 32 tentative attributes left: ALDH2, ATP6V1G1, C16orf45, CDC42BPA, ## COL4A6 and 27 more;查看下變量重要性鑒定結果（實際上面的輸出中也已經有體現了），54個重要的變量，36個可能重要的變量 (tentative variable, 重要性得分與最好的影子變量得分無統計差異)，6,980個不重要的變量。
table(boruta$finalDecision)## ## Tentative Confirmed Rejected ## 32 46 8943繪制鑒定出的變量的重要性。變量少了可以用默認繪圖，變量多時繪制的圖看不清，需要自己整理數據繪圖。
定義一個函數提取每個變量對應的重要性值。
library(dplyr) boruta.imp <- function(x){imp <- reshape2::melt(x$ImpHistory, na.rm=T)[,-1]colnames(imp) <- c("Variable","Importance")imp <- imp[is.finite(imp$Importance),]variableGrp <- data.frame(Variable=names(x$finalDecision), finalDecision=x$finalDecision)showGrp <- data.frame(Variable=c("shadowMax", "shadowMean", "shadowMin"),finalDecision=c("shadowMax", "shadowMean", "shadowMin"))variableGrp <- rbind(variableGrp, showGrp)boruta.variable.imp <- merge(imp, variableGrp, all.x=T)sortedVariable <- boruta.variable.imp %>% group_by(Variable) %>% summarise(median=median(Importance)) %>% arrange(median)sortedVariable <- as.vector(sortedVariable$Variable)boruta.variable.imp$Variable <- factor(boruta.variable.imp$Variable, levels=sortedVariable)invisible(boruta.variable.imp) }boruta.variable.imp <- boruta.imp(boruta)head(boruta.variable.imp)## Variable Importance finalDecision ## 1 AADAC 0 Rejected ## 2 AADAC 0 Rejected ## 3 AADAC 0 Rejected ## 4 AADAC 0 Rejected ## 5 AADAC 0 Rejected ## 6 AADAC 0 Rejected只繪制Confirmed變量。
library(ImageGP)sp_boxplot(boruta.variable.imp, melted=T, xvariable = "Variable", yvariable = "Importance",legend_variable = "finalDecision", legend_variable_order = c("shadowMax", "shadowMean", "shadowMin", "Confirmed"),xtics_angle = 90)提取重要的變量和可能重要的變量
boruta.finalVarsWithTentative <- data.frame(Item=getSelectedAttributes(boruta, withTentative = T), Type="Boruta_with_tentative")看下這些變量的值的分布
caret::featurePlot(train_data[,boruta.finalVarsWithTentative$Item], train_data_group, plot="box")交叉驗證選擇參數并擬合模型
定義一個函數生成一些列用來測試的mtry (一系列不大于總變量數的數值)。
generateTestVariableSet <- function(num_toal_variable){max_power <- ceiling(log10(num_toal_variable))tmp_subset <- c(unlist(sapply(1:max_power, function(x) (1:10)^x, simplify = F)), ceiling(max_power/3))#return(tmp_subset)base::unique(sort(tmp_subset[tmp_subset<num_toal_variable])) } # generateTestVariableSet(78)選擇關鍵特征變量相關的數據
# 提取訓練集的特征變量子集 boruta_train_data <- train_data[, boruta.finalVarsWithTentative$Item] boruta_mtry <- generateTestVariableSet(ncol(boruta_train_data))使用 Caret 進行調參和建模
library(caret) # Create model with default parameters trControl <- trainControl(method="repeatedcv", number=10, repeats=5)seed <- 1 set.seed(seed) # 根據經驗或感覺設置一些待查詢的參數和參數值 tuneGrid <- expand.grid(mtry=boruta_mtry)borutaConfirmed_rf_default <- train(x=boruta_train_data, y=train_data_group, method="rf", tuneGrid = tuneGrid, # metric="Accuracy", #metric='Kappa'trControl=trControl) borutaConfirmed_rf_default## Random Forest ## ## 77 samples ## 78 predictors ## 2 classes: 'normal', 'tumor' ## ## No pre-processing ## Resampling: Cross-Validated (10 fold, repeated 5 times) ## Summary of sample sizes: 71, 69, 69, 69, 69, 69, ... ## Resampling results across tuning parameters: ## ## mtry Accuracy Kappa ## 1 0.9352381 0.8708771 ## 2 0.9352381 0.8708771 ## 3 0.9352381 0.8708771 ## 4 0.9377381 0.8758771 ## 5 0.9377381 0.8758771 ## 6 0.9402381 0.8808771 ## 7 0.9402381 0.8808771 ## 8 0.9452381 0.8908771 ## 9 0.9402381 0.8808771 ## 10 0.9452381 0.8908771 ## 16 0.9452381 0.8908771 ## 25 0.9477381 0.8958771 ## 36 0.9452381 0.8908771 ## 49 0.9402381 0.8808771 ## 64 0.9327381 0.8658771 ## ## Accuracy was used to select the optimal model using the largest value. ## The final value used for the model was mtry = 25.可視化不同參數的準確性分布
plot(borutaConfirmed_rf_default)可視化Top20重要的變量
dotPlot(varImp(borutaConfirmed_rf_default))提取最終選擇的模型，并繪制 ROC 曲線評估模型
borutaConfirmed_rf_default_finalmodel <- borutaConfirmed_rf_default$finalModel先自評，評估模型對訓練集的分類效果
采用訓練數據集評估構建的模型，Accuracy=1; Kappa=1，非常完美。
模型的預測顯著性P-Value [Acc > NIR] : 2.2e-16。其中NIR是No Information Rate，其計算方式為數據集中最大的類包含的數據占總數據集的比例。如某套數據中，分組A有80個樣品，分組B有20個樣品，我們只要猜A，正確率就會有80%，這就是NIR。如果基于這套數據構建的模型準確率也是80%，那么這個看上去準確率較高的模型也沒有意義。confusionMatrix使用binom.test函數檢驗模型的準確性Accuracy是否顯著優于NIR，若P-value<0.05，則表示模型預測準確率顯著高于隨便猜測。
# 獲得模型結果評估矩陣(`confusion matrix`)predictions_train <- predict(borutaConfirmed_rf_default_finalmodel, newdata=train_data) confusionMatrix(predictions_train, train_data_group)## Confusion Matrix and Statistics ## ## Reference ## Prediction normal tumor ## normal 38 0 ## tumor 0 39 ## ## Accuracy : 1 ## 95% CI : (0.9532, 1) ## No Information Rate : 0.5065 ## P-Value [Acc > NIR] : < 2.2e-16 ## ## Kappa : 1 ## ## Mcnemar's Test P-Value : NA ## ## Sensitivity : 1.0000 ## Specificity : 1.0000 ## Pos Pred Value : 1.0000 ## Neg Pred Value : 1.0000 ## Prevalence : 0.4935 ## Detection Rate : 0.4935 ## Detection Prevalence : 0.4935 ## Balanced Accuracy : 1.0000 ## ## 'Positive' Class : normal ##盲評，評估模型應用于測試集時的效果
繪制ROC曲線，計算模型整體的AUC值，并選擇最佳模型。
# 繪制ROC曲線prediction_prob <- predict(borutaConfirmed_rf_default_finalmodel, newdata=test_data, type="prob") library(pROC) roc_curve <- roc(test_data_group, prediction_prob[,1])roc_curve## ## Call: ## roc.default(response = test_data_group, predictor = prediction_prob[, 1]) ## ## Data: prediction_prob[, 1] in 12 controls (test_data_group normal) > 13 cases (test_data_group tumor). ## Area under the curve: 0.9872# roc <- roc(test_data_group, factor(predictions, ordered=T)) # plot(roc)基于默認閾值的盲評
基于默認閾值繪制混淆矩陣并評估模型預測準確度顯著性,結果顯著P-Value [Acc > NIR]<0.05。
# 獲得模型結果評估矩陣(`confusion matrix`)predictions <- predict(borutaConfirmed_rf_default_finalmodel, newdata=test_data) confusionMatrix(predictions, test_data_group)## Confusion Matrix and Statistics ## ## Reference ## Prediction normal tumor ## normal 12 2 ## tumor 0 11 ## ## Accuracy : 0.92 ## 95% CI : (0.7397, 0.9902) ## No Information Rate : 0.52 ## P-Value [Acc > NIR] : 2.222e-05 ## ## Kappa : 0.8408 ## ## Mcnemar's Test P-Value : 0.4795 ## ## Sensitivity : 1.0000 ## Specificity : 0.8462 ## Pos Pred Value : 0.8571 ## Neg Pred Value : 1.0000 ## Prevalence : 0.4800 ## Detection Rate : 0.4800 ## Detection Prevalence : 0.5600 ## Balanced Accuracy : 0.9231 ## ## 'Positive' Class : normal ##選擇模型分類最佳閾值再盲評
youden:max(sensitivities?+?r?×?specificities)
closest.topleft:min((1???sensitivities)2?+?r?×?(1???specificities)2)
r是加權系數，默認是1，其計算方式為r?=?(1???prevalenc**e)/(cos**t?*?prevalenc**e).
best.weights控制加權方式：(cost, prevalence)默認是(1,0.5)，據此算出的r為1。
cost: 假陰性率占假陽性率的比例，容忍更高的假陽性率還是假陰性率
prevalence: 關注的類中的個體所占的比例 (n.cases/(n.controls+n.cases)).
準備數據繪制ROC曲線
library(ggrepel) ROC_data <- data.frame(FPR = 1- roc_curve$specificities, TPR=roc_curve$sensitivities) ROC_data <- ROC_data[with(ROC_data, order(FPR,TPR)),]p <- ggplot(data=ROC_data, mapping=aes(x=FPR, y=TPR)) +geom_step(color="red", size=1, direction = "vh") +geom_segment(aes(x=0, xend=1, y=0, yend=1)) + theme_classic() + xlab("False positive rate") + ylab("True positive rate") + coord_fixed(1) + xlim(0,1) + ylim(0,1) +annotate('text', x=0.5, y=0.25, label=paste('AUC=', round(roc_curve$auc,2))) +geom_point(data=best_thresh, mapping=aes(x=1-specificity, y=sensitivity), color='blue', size=2) + geom_text_repel(data=best_thresh, mapping=aes(x=1.05-specificity, y=sensitivity ,label=best)) p基于選定的最優閾值制作混淆矩陣并評估模型預測準確度顯著性,結果顯著P-Value [Acc > NIR]<0.05。
predict_result <- data.frame(Predict_status=c(T,F), Predict_class=colnames(prediction_prob))head(predict_result)## Predict_status Predict_class ## 1 TRUE normal ## 2 FALSE tumorpredictions2 <- plyr::join(data.frame(Predict_status=prediction_prob[,1] > best_thresh[1,1]), predict_result)predictions2 <- as.factor(predictions2$Predict_class)confusionMatrix(predictions2, test_data_group)## Confusion Matrix and Statistics ## ## Reference ## Prediction normal tumor ## normal 11 0 ## tumor 1 13 ## ## Accuracy : 0.96 ## 95% CI : (0.7965, 0.999) ## No Information Rate : 0.52 ## P-Value [Acc > NIR] : 1.913e-06 ## ## Kappa : 0.9196 ## ## Mcnemar's Test P-Value : 1 ## ## Sensitivity : 0.9167 ## Specificity : 1.0000 ## Pos Pred Value : 1.0000 ## Neg Pred Value : 0.9286 ## Prevalence : 0.4800 ## Detection Rate : 0.4400 ## Detection Prevalence : 0.4400 ## Balanced Accuracy : 0.9583 ## ## 'Positive' Class : normal ##機器學習系列教程
從隨機森林開始，一步步理解決策樹、隨機森林、ROC/AUC、數據集、交叉驗證的概念和實踐。
文字能說清的用文字、圖片能展示的用、描述不清的用公式、公式還不清楚的寫個簡單代碼，一步步理清各個環節和概念。
再到成熟代碼應用、模型調參、模型比較、模型評估，學習整個機器學習需要用到的知識和技能。
一圖感受各種機器學習算法
機器學習算法 - 隨機森林之決策樹初探（1）
機器學習算法-隨機森林之決策樹R 代碼從頭暴力實現（2）
機器學習算法-隨機森林之決策樹R 代碼從頭暴力實現（3）
機器學習算法-隨機森林之理論概述
機器學習算法-隨機森林初探（1）
機器學習 - 隨機森林手動10 折交叉驗證
機器學習 模型評估指標 - ROC曲線和AUC值
機器學習 - 訓練集、驗證集、測試集
一個函數統一238個機器學習R包，這也太贊了吧
基于Caret和RandomForest包進行隨機森林分析的一般步驟 （1）
Caret模型訓練和調參更多參數解讀(2)
基于Caret進行隨機森林隨機調參的4種方式
機器學習第17篇 - 特征變量篩選（1）
機器學習第18篇 - Boruta特征變量篩選（2）
機器學習第19篇 - 機器學習系列補充：數據集準備和更正YSX包
機器學習第20篇 - 基于Boruta選擇的特征變量構建隨機森林
機器學習第21篇 - 特征遞歸消除RFE算法 理論
機器學習第22篇 - RFE篩選出的特征變量竟然是Boruta的4倍之多
機器學習第23篇 - 更多特征變量卻未能帶來隨機森林分類效果的提升
機器學習相關書籍分享
UCI機器學習數據集
送你一個在線機器學習網站，真香！
多套用于機器學習的多種癌癥表達數據集
這個統一了238個機器學習模型R包的參考手冊推薦給你
莫煩Python機器學習
機器學習與人工智能、深度學習有什么關系？終于有人講明白了
往期精品(點擊圖片直達文字對應教程)
機器學習
后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集
總結
以上是生活随笔為你收集整理的一套完整的基于随机森林的机器学习流程（特征选择、交叉验证、模型评估））...的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 监督学习 | CART 分类回归树原理
- 下一篇： m进制数转换为十进制数