小伙伴做生存分析的時候避免不了構(gòu)建Cox回歸風(fēng)險模型，小果告訴大家一個新的方式，小伙伴可以不再用Cox回歸了，我們?nèi)ピ囋嚈C械學(xué)習(xí)的方法之一隨機森林的生存分析也就是randomForestSRC。我們先了解一下方法， 隨機生存森林呢，其實是通過訓(xùn)練大量生存樹，以表決的形式，從個體樹之中加權(quán)選舉出最終的預(yù)測結(jié)果。

構(gòu)建隨機生存森林呢一般有下面幾個流程

Ⅰ. 模型通過“自助法”（Bootstrap）將原始數(shù)據(jù)以有放回的形式隨機抽取樣本，建立樣本子集，并將每個樣本中37%的數(shù)據(jù)作為袋外數(shù)據(jù)（Out-of-Bag Data）排除在外； Ⅱ. 對每一個樣本隨機選擇特征構(gòu)建其對應(yīng)的生存樹； Ⅲ. 利用Nelson-Aalen法估計隨機生存森林模型的總累積風(fēng)險； Ⅳ. 使用袋外數(shù)據(jù)計算模型準(zhǔn)確度。

下面我們開始吧，我們先按照一下R包，沒有的小伙伴先安裝一下哦

if (!require(randomForestSRC)) install.packages("randomForestSRC") if (!require(survival)) install.packages("survival") ? library(randomForestSRC) library(survival) 我們先讀取下示例數(shù)據(jù)看看，這是R包自帶的

下面我們用這幾個數(shù)據(jù)進(jìn)行學(xué)習(xí)和展示

data(veteran, package = "randomForestSRC") head(veteran)

data(pbc, package = "randomForestSRC") head(pbc)

data(wihs, package = "randomForestSRC") head(wihs)

我們使用肺癌的數(shù)據(jù)集進(jìn)行分析 肺癌的兩種治療方案的是隨機試驗 例子一 v.obj <- rfsrc(Surv(time, status) ~ ., data = veteran, ntree = 100, nsplit = 10, ???????????????na.action = "na.impute", tree.err = TRUE, importance = TRUE, block.size = 1) ###設(shè)置樹的個數(shù)為3： ## plot tree number 3 plot(get.tree(v.obj, 3))

###繪制訓(xùn)練森林的結(jié)果 ## plot results of trained forest plot(v.obj)

我們直接繪制前10個個體的生存曲線 #繪制前十個個體生存曲線 matplot(v.obj$time.interest, 100 * t(v.obj$survival.oob[1:10, ]), xlab = "Time", ????????ylab = "Survival", type = "l", lty = 1)

##使用函數(shù)plot.survival繪制前10個個體的生存曲線 plot.survival(v.obj, subset = 1:10)

下面我們###快速優(yōu)化老數(shù)據(jù)的節(jié)點大小##最優(yōu)的生存節(jié)點大小比其他家族# tune.nodesize(Surv(time, status) ~ ., veteran)下面我們###快速優(yōu)化老數(shù)據(jù)的節(jié)點大小##最優(yōu)的生存節(jié)點大小比其他家族# ? tune.nodesize(Surv(time, status) ~ ., veteran)

例子2 vd <- veteran vd$celltype = factor(vd$celltype) vd$diagtime = factor(vd$diagtime) vd.obj <- rfsrc(Surv(time, status) ~ ., vd, ntree = 100, nodesize = 5) plot(get.tree(vd.obj, 3))

下面我們使用疾病的例子演示 原發(fā)性膽汁性肝硬化PBC 來看一下 pbc.obj <- rfsrc(Surv(days, status) ~ ., pbc) print(pbc.obj)

# pbc.obj2 <- rfsrc(Surv(days, status) ~ ., pbc, nsplit = 10, na.action = "na.impute") ? ##與上面相同，但迭代丟失的數(shù)據(jù)算法 pbc.obj3 <- rfsrc(Surv(days, status) ~ ., pbc, na.action = "na.impute", nimpute = 3) ? ## 估算數(shù)據(jù)的快速方法（不進(jìn)行推斷） pbc.imp <- impute(Surv(days, status) ~ ., pbc, splitrule = "random") 我們對兩種方法比較一下 比較RF-SRC和Cox回歸，說明性能的c指數(shù)和Brier評分措施，假設(shè)加載了“pec”和“survival” ##比較RF-SRC和COX require("survival") require("pec") require("prodlim") ## pec所需的預(yù)測函數(shù) predictSurvProb.rfsrc <- function(object, newdata, times, ...) { ??ptemp <- predict(object, newdata = newdata, ...)$survival ??pos <- sindex(jump.times = object$time.interest, eval.times = times) ??p <- cbind(1, ptemp)[, pos + 1] ??if (NROW(p) != NROW(newdata) || NCOL(p) != length(times)) ????stop("Prediction failed") ??p } ## 數(shù)據(jù)，公式規(guī)范 data(pbc, package = "randomForestSRC") pbc.na <- na.omit(pbc)?##remove NA's surv.f <- as.formula(Surv(days, status) ~ .) pec.f <- as.formula(Hist(days, status) ~ 1) ? ## 運行cox/rfsrc模型進(jìn)行說明我們使用了少量的樹 cox.obj <- coxph(surv.f, data = pbc.na, x = TRUE) rfsrc.obj <- rfsrc(surv.f, pbc.na, ntree = 150) ## 計算期望Brier分?jǐn)?shù)的bootstrap交叉驗證估計 ? set.seed(17743) prederror.pbc <- pec(list(cox.obj, rfsrc.obj), data = pbc.na, formula = pec.f, splitMethod = "bootcv", ?????????????????????B = 50) print(prederror.pbc) ## ? plot(prederror.pbc)

下面我們計算Cox回歸 #Cox ## 計算cox回歸的袋外C指數(shù)并與rfsrc進(jìn)行比較 rfsrc.obj <- rfsrc(surv.f, pbc.na) cat("out-of-bag Cox Analysis ...", "\n") ## out-of-bag Cox Analysis ... cox.err <- sapply(1:100, function(b) { ??if (b%%10 == 0) ????cat("cox bootstrap:", b, "\n") ??train <- sample(1:nrow(pbc.na), nrow(pbc.na), replace = TRUE) ??cox.obj <- tryCatch({ ????coxph(surv.f, pbc.na[train, ]) ??}, error = function(ex) { ????NULL ??}) ??if (!is.null(cox.obj)) { ????get.cindex(pbc.na$days[-train], pbc.na$status[-train], predict(cox.obj, pbc.na[-train, ????])) ??} else NA }) ? cat("E:\\SXG") ## ## OOB error rates cat("\tRSF???????????: ", rfsrc.obj$err.rate[rfsrc.obj$ntree], "\n") ## RSF???????????:?0.1714494 cat("\tCox regression : ", mean(cox.err, na.rm = TRUE), "\n")

下面我們對生存競爭風(fēng)險比例模型進(jìn)行構(gòu)建 來看一下 #生存競爭風(fēng)險比例 wihs.obj <- rfsrc(Surv(time, status) ~ ., wihs, nsplit = 3, ntree = 100) plot.competing.risk(wihs.obj)

cif <- wihs.obj$cif.oob Time <- wihs.obj$time.interest idu <- wihs$idu cif.haart <- cbind(apply(cif[, , 1][idu == 0, ], 2, mean), apply(cif[, , 1][idu == ??????????????????????????????????????????????????????????????????????????????1, ], 2, mean)) cif.aids <- cbind(apply(cif[, , 2][idu == 0, ], 2, mean), apply(cif[, , 2][idu == ?????????????????????????????????????????????????????????????????????????????1, ], 2, mean)) matplot(Time, cbind(cif.haart, cif.aids), type = "l", lty = c(1, 2, 1, 2), col = c(4, ???????????????????????????????????????????????????????????????????????????????????4, 2, 2), lwd = 3, ylab = "Cumulative Incidence") legend("bottomright", legend = c("HAART (Non-IDU)", "HAART (IDU)", "AIDS (Non-IDU)", ?????????????????????????????????"AIDS (IDU)"), lty = c(1, 2, 1, 2), col = c(4, 4, 2, 2), lwd = 3, cex = 0.6)#圖注

小伙伴可能不太懂結(jié)果的含義，小果這里簡單解讀一些

隨機生存森林可以對變量重要性進(jìn)行排名，VIMP法和最小深度法是最常用的方法：變量VIMP值小于0說明該變量降低了預(yù)測的準(zhǔn)確性，而當(dāng)VIMP值大于0則說明該變量提高了預(yù)測的準(zhǔn)確性；最小深度法通過計算運行到最終節(jié)點時的最小深度來給出各變量對于結(jié)局事件的重要性。 上述為綜合兩種方法的散點圖，其中，藍(lán)色點代表VIMP值大于0，紅色則代表VIMP值小于0；在紅色對角虛線上的點代表兩種方法對該變量的排名相同，高于對角虛線的點代表其VIMP排名更高，低于對角虛線的點則代表其最小深度排名更高。相較于Cox比例風(fēng)險回歸模型等傳統(tǒng)生存分析方法，隨機生存森林模型的預(yù)測準(zhǔn)確度至少等同或優(yōu)于傳統(tǒng)生存分析方法。 隨機森林呢其實有很多的優(yōu)勢，小果這里就不一一介紹了，小伙伴學(xué)會了嗎，將自己的數(shù)據(jù)去做個實驗，動手試試這個新方法吧！

小伙伴多多理解代碼的意義哦~