機(jī)器學(xué)習(xí)算法-隨機(jī)森林初探（1）
機(jī)器學(xué)習(xí)算法-隨機(jī)森林之理論概述
表達(dá)數(shù)據(jù)集來源于 https://file.biolab.si/biolab/supp/bi-cancer/projections/。
為了展示隨機(jī)森林的能力，我們用一套早期的淋巴瘤基因表達(dá)芯片數(shù)據(jù)集，包含77個(gè)樣品，2個(gè)分組和7070個(gè)變量。
讀入數(shù)據(jù)
expr_file <- "DLBCL.expr.txt"
metadata_file <- "DLBCL.metadata.txt"
# 每個(gè)基因表達(dá)值是內(nèi)部比較，只要是樣品之間標(biāo)準(zhǔn)化的數(shù)據(jù)即可，其它什么轉(zhuǎn)換都關(guān)系不大
expr_mat <- read.table(expr_file, row.names = 1, header = T, sep="\t")
metadata <- read.table(metadata_file, row.names=1, header=T, sep="\t")
dim(expr_mat)
## [1] 7070 77基因表達(dá)表
expr_mat[1:4,1:5]
## DLBCL_1 DLBCL_2 DLBCL_3 DLBCL_4 DLBCL_5
## A28102 -1 25 73 267 16
## AB000114_at -45 -17 91 41 24
## AB000115_at 176 531 257 202 187
## AB000220_at 97 353 80 138 39Metadata表
head(metadata)
## class
## DLBCL_1 DLBCL
## DLBCL_2 DLBCL
## DLBCL_3 DLBCL
## DLBCL_4 DLBCL
## DLBCL_5 DLBCL
## DLBCL_6 DLBCL樣品篩選和排序
對(duì)讀入的表達(dá)數(shù)據(jù)進(jìn)行轉(zhuǎn)置。通常我們是一行一個(gè)基因，一列一個(gè)樣品。在構(gòu)建模型時(shí)，數(shù)據(jù)通常是反過來的，一列一個(gè)基因，一行一個(gè)樣品。每一列代表一個(gè)變量
(variable)，每一行代表一個(gè)案例
(case)。這樣更方便提取每個(gè)變量，且易于把模型中的x,y放到一個(gè)矩陣中。
樣本表和表達(dá)表中的樣本順序?qū)R一致也是需要確保的一個(gè)操作。
# 表達(dá)數(shù)據(jù)轉(zhuǎn)置
# 習(xí)慣上我們是一行一個(gè)基因，一列一個(gè)樣品
# 做機(jī)器學(xué)習(xí)時(shí)，大部分?jǐn)?shù)據(jù)都是反過來的，一列一個(gè)基因，一行一個(gè)樣品
# 每一列代表一個(gè)變量
expr_mat <- t(expr_mat)
expr_mat_sampleL <- rownames(expr_mat)
metadata_sampleL <- rownames(metadata)
common_sampleL <- intersect(expr_mat_sampleL, metadata_sampleL)
# 保證表達(dá)表樣品與METAdata樣品順序和數(shù)目完全一致
expr_mat <- expr_mat[common_sampleL,,drop=F]
metadata <- metadata[common_sampleL,,drop=F]判斷是分類還是回歸
如果group對(duì)應(yīng)的列為數(shù)字，轉(zhuǎn)換為數(shù)值型 - 做回歸
如果group對(duì)應(yīng)的列為分組，轉(zhuǎn)換為因子型 - 做分類
# R4.0之后默認(rèn)讀入的不是factor，需要做一個(gè)轉(zhuǎn)換
# devtools::install_github("Tong-Chen/YSX")
library(YSX)
# 此處的class根據(jù)需要修改
group = "class"
# 如果group對(duì)應(yīng)的列為數(shù)字，轉(zhuǎn)換為數(shù)值型 - 做回歸
# 如果group對(duì)應(yīng)的列為分組，轉(zhuǎn)換為因子型 - 做分類
if(numCheck(metadata[[group]])){
if (!is.numeric(metadata[[group]])) {
metadata[[group]] <- mixedToFloat(metadata[[group]])
}
} else{
metadata[[group]] <- as.factor(metadata[[group]])
}隨機(jī)森林初步分析
library(randomForest)
# 查看參數(shù)是個(gè)好習(xí)慣
# 有了前面的基礎(chǔ)概述，再看每個(gè)參數(shù)的含義就明確了很多
# 也知道該怎么調(diào)了
# 每個(gè)人要解決的問題不同，通常不是別人用什么參數(shù)，自己就跟著用什么參數(shù)
# 尤其是到下游分析時(shí)
# ?randomForest加載包之后，直接分析一下，看到結(jié)果再調(diào)參。（竟然被awk生成的隨機(jī)數(shù)給整蒙了，也談隨機(jī)數(shù)生成種子）
# 設(shè)置隨機(jī)數(shù)種子，具體含義見 https://mp.weixin.qq.com/s/6plxo-E8qCdlzCgN8E90zg
set.seed(304)
# 直接使用默認(rèn)參數(shù)
rf <- randomForest(expr_mat, metadata$class)查看下初步結(jié)果,
隨機(jī)森林類型判斷為分類，構(gòu)建了500棵樹，每次決策時(shí)使用了84個(gè)基因，OOB估計(jì)的錯(cuò)誤率是12.99%，挺高的。
分類效果評(píng)估矩陣Confusion matrix,顯示DLBCL組的分類錯(cuò)誤率為0.034，FL組的分類錯(cuò)誤率為0.42。這是隨機(jī)森林的默認(rèn)操作，樣本量多的分類錯(cuò)誤率會(huì)低
(后面我們?cè)僬{(diào)整)。
rf
##
## Call:
## randomForest(x = expr_mat, y = metadata$class)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 84
##
## OOB estimate of error rate: 12.99%
## Confusion matrix:
## DLBCL FL class.error
## DLBCL 56 2 0.03448276
## FL 8 11 0.42105263隨機(jī)森林參數(shù)調(diào)試 - 決策樹的數(shù)目
增加決策樹的數(shù)目到1000測(cè)試下分類率是否會(huì)降低。OOB估計(jì)的錯(cuò)誤率是11.69%，還是不低。
分類效果評(píng)估矩陣Confusion matrix,顯示DLBCL組的分類錯(cuò)誤率為0.017，FL組的分類錯(cuò)誤率為0.42。也都略有降低。
set.seed(304)
rf1000 <- randomForest(expr_mat, metadata$class, ntree=1000)
rf1000
##
## Call:
## randomForest(x = expr_mat, y = metadata$class, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 84
##
## OOB estimate of error rate: 11.69%
## Confusion matrix:
## DLBCL FL class.error
## DLBCL 57 1 0.01724138
## FL 8 11 0.42105263增加決策樹的數(shù)目到2000測(cè)試下分類率是否會(huì)降低。OOB估計(jì)的錯(cuò)誤率是11.69%沒有變化。
分類效果評(píng)估矩陣Confusion matrix,顯示DLBCL組的分類錯(cuò)誤率為0.017，FL組的分類錯(cuò)誤率為0.42。FL組樣品量少，分類效果還是不好。
set.seed(304)
rf2000 <- randomForest(expr_mat, metadata$class, ntree=2000)
rf2000
##
## Call:
## randomForest(x = expr_mat, y = metadata$class, ntree = 2000)
## Type of random forest: classification
## Number of trees: 2000
## No. of variables tried at each split: 84
##
## OOB estimate of error rate: 11.69%
## Confusion matrix:
## DLBCL FL class.error
## DLBCL 57 1 0.01724138
## FL 8 11 0.42105263繪制下從第1到2000棵樹時(shí)，OOB的變化趨勢(shì)是怎樣的？從這張圖可以看到，600棵樹之后，基本沒有影響了。
library(ggplot2)
YSX::sp_lines(as.data.frame(rf2000$err.rate), manual_color_vector="Set2",
x_label="Number of trees", y_label="Error rate",line_size=0.6,
width=6, height=4
)
隨機(jī)森林參數(shù)調(diào)試 - 用于決策的變量的數(shù)目 (m)
影響隨機(jī)森林的第二個(gè)參數(shù)(m):
構(gòu)建每棵決策樹時(shí)隨機(jī)抽取的變量的數(shù)目。（假設(shè)你已閱讀過了機(jī)器學(xué)習(xí)算法-隨機(jī)森林之理論概述）。在randomForest函數(shù)中有一個(gè)參數(shù)mtry即是做這個(gè)的。在處理分類問題時(shí)，其默認(rèn)值為sqrt(p)；處理回歸問題時(shí)，其默認(rèn)值為p/3；p是總的變量數(shù)目。
我們先人工測(cè)試幾個(gè)不同的mtry看下效果。mtry=20時(shí)錯(cuò)誤率為15.58%。
# 增加樹的數(shù)目沒有給出好的結(jié)果，這里還是用的默認(rèn)的500棵樹以便獲得較快的運(yùn)行速度
set.seed(304)
rf_mtry20 <- randomForest(expr_mat, metadata$class, mtry=20)
rf_mtry20
##
## Call:
## randomForest(x = expr_mat, y = metadata$class, mtry = 20)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 20
##
## OOB estimate of error rate: 15.58%
## Confusion matrix:
## DLBCL FL class.error
## DLBCL 57 1 0.01724138
## FL 11 8 0.57894737我們先人工測(cè)試幾個(gè)不同的mtry看下效果。mtry=50時(shí)錯(cuò)誤率為11.69%(默認(rèn)使用了84個(gè)變量，準(zhǔn)確率為12.99)。
# 增加樹的數(shù)目沒有給出好的結(jié)果，這里還是用的默認(rèn)的500棵樹以便獲得較快的運(yùn)行速度
set.seed(304)
rf_mtry50 <- randomForest(expr_mat, metadata$class, mtry=50)
rf_mtry50
##
## Call:
## randomForest(x = expr_mat, y = metadata$class, mtry = 50)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 50
##
## OOB estimate of error rate: 11.69%
## Confusion matrix:
## DLBCL FL class.error
## DLBCL 57 1 0.01724138
## FL 8 11 0.42105263我們先人工測(cè)試幾個(gè)不同的mtry看下效果。mtry=100時(shí)錯(cuò)誤率為11.69%。
# 增加樹的數(shù)目沒有給出好的結(jié)果，這里還是用的默認(rèn)的500棵樹以便獲得較快的運(yùn)行速度
set.seed(304)
rf_mtry100 <- randomForest(expr_mat, metadata$class, mtry=100)
rf_mtry100
##
## Call:
## randomForest(x = expr_mat, y = metadata$class, mtry = 100)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 100
##
## OOB estimate of error rate: 11.69%
## Confusion matrix:
## DLBCL FL class.error
## DLBCL 57 1 0.01724138
## FL 8 11 0.42105263一個(gè)個(gè)測(cè)試也不是辦法，tuneRF給我們提供了一個(gè)根據(jù)OOB值迭代鑒定最合適的mtry值的函數(shù)。(測(cè)試幾次，不太好用，改一個(gè)stepfactor，結(jié)果變化很大)
# mtryStart: 從多少個(gè)變量開始嘗試，可用默認(rèn)值。程序會(huì)自動(dòng)向更多變量或更少變量迭代。
# ntreeTry: 迭代時(shí)構(gòu)建多少棵樹，這里設(shè)置為500，與我們上面獲得的效果比較好的樹一致。
# stepFactor: 迭代步長，mtryStart向更多變量迭代時(shí)乘以此值；mtryStart向更少變量迭代時(shí)除以此值。
# improve：如果迭代不能給OOB帶來給定值的改善，則停止迭代。
set.seed(304)
tuneRF(expr_mat, metadata$class, ntreeTry=500, stepFactor=1.1, improve=1e-5)
## mtry = 84 OOB error = 12.99%
## Searching left ...
## mtry = 77 OOB error = 10.39%
## 0.2 1e-05
## mtry = 70 OOB error = 11.69%
## -0.125 1e-05
## Searching right ...
## mtry = 92 OOB error = 11.69%
## -0.125 1e-05
## mtry OOBError
## 70.OOB 70 0.1168831
## 77.OOB 77 0.1038961
## 84.OOB 84 0.1298701
## 92.OOB 92 0.1168831randomForest自帶了另一個(gè)函數(shù)rfcv，通過嵌套交叉驗(yàn)證方式評(píng)估了根據(jù)變量重要性降低預(yù)測(cè)變量后的模型的性能。
result = rfcv(expr_mat, metadata$class, cv.fold=10)
result$error.cv
## 7070 3535 1768 884 442 221 110 55 28 14 7
## 0.10389610 0.11688312 0.10389610 0.09090909 0.09090909 0.09090909 0.09090909 0.10389610 0.10389610 0.12987013 0.14285714
## 3 1
## 0.15584416 0.15584416后面我們?cè)敿?xì)講下交叉驗(yàn)證這個(gè)最常用的評(píng)估模型的方法。
往期精品(點(diǎn)擊圖片直達(dá)文字對(duì)應(yīng)教程)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
后臺(tái)回復(fù)“生信寶典福利第一波”或點(diǎn)擊閱讀原文獲取教程合集
?
（請(qǐng)備注姓名-學(xué)校/企業(yè)-職務(wù)等）