生存分析-隨機(jī)森林實(shí)驗(yàn)與代碼

1、-. z.隨機(jī)森林模型在生存分析中的應(yīng)用【摘要】目的：本文探討隨機(jī)森林方法用于高維度、強(qiáng)相關(guān)、小樣本的生存資料分析時(shí)，可以起到變量篩選的作用。方法：以乳腺癌數(shù)據(jù)集構(gòu)建乳腺癌轉(zhuǎn)移風(fēng)險(xiǎn)評(píng)估模型為實(shí)例進(jìn)展實(shí)證分析，使用隨機(jī)森林模型進(jìn)展變量選擇，然后擬合co*回歸模型。結(jié)果：隨機(jī)森林模型通過(guò)對(duì)變量的選擇，有效的解決數(shù)據(jù)維度高且強(qiáng)相關(guān)的情況，得到了較高的AUC值。一、數(shù)據(jù)說(shuō)明該乳腺癌數(shù)據(jù)集來(lái)自于NCBI，有77個(gè)觀測(cè)值以及22286個(gè)基因變量。通過(guò)篩選選取454個(gè)基因變量。將數(shù)據(jù)隨機(jī)分為訓(xùn)練集合測(cè)試集，其中2/3為訓(xùn)練集，1/3為測(cè)試集。繪制K-M曲線圖：二、隨機(jī)森林模型隨機(jī)森林由許多的決策樹(shù)組成，因?yàn)?/p>

2、這些決策樹(shù)的形成采用了隨機(jī)的方法，因此也叫做隨機(jī)決策樹(shù)。隨機(jī)森林中的樹(shù)之間是沒(méi)有關(guān)聯(lián)的。當(dāng)測(cè)試數(shù)據(jù)進(jìn)入隨機(jī)森林時(shí)，其實(shí)就是讓每一顆決策樹(shù)進(jìn)展分類，最后取所有決策樹(shù)中分類結(jié)果最多的那類為最終的結(jié)果。因此隨機(jī)森林是一個(gè)包含多個(gè)決策樹(shù)的分類器，并且其輸出的類別是由個(gè)別樹(shù)輸出的類別的眾數(shù)而定。使用 randomForestSRC包得到的隨機(jī)森林模型具有以下性質(zhì)： Number of deaths: 27 Number of trees: 800 Minimum terminal node size: 3 Average no. of terminal nodes: 14.4275No. of vari

3、ables tried at each split: 3 Total no. of variables: 452 Analysis: RSF Family: surv Splitting rule: logrank Error rate: 19.87%發(fā)現(xiàn)直接使用隨機(jī)森林得到的模型，預(yù)測(cè)誤差很大，到達(dá)了19.8%,進(jìn)一步考慮使用隨機(jī)森林模型進(jìn)展變量選擇，結(jié)果如下： Sample size: 52 Number of deaths: 19 Number of trees: 500 Minimum terminal node size: 2 Average no. of terminal node

4、s: 11.554No. of variables tried at each split: 3 Total no. of variables: 9 Analysis: RSF Family: surv Splitting rule: logrank *random* Number of random split points: 10 Error rate: 11.4% our.rf$topvars1 213821_s_at 219778_at 204690_at 220788_s_at 202202_s_at6 211603_s_at 213055_at 219336_s_at 37892_

5、at 一共選取了9個(gè)變量，同時(shí)誤差只有11.4%接下來(lái)，使用這些變量做co*回歸，剔除模型中不顯著0.01的變量，最終參與模型建立的變量共有4個(gè)。模型結(jié)果如下： e*p(coef) e*p(-coef) lower .95 upper .95218150_at 1.6541 0.6046 0.11086 24.6800200914_*_at 0.9915 1.0086 0.34094 2.8833220788_s_at 0.2649 3.7750 0.05944 1.1805201398_s_at 1.7457 0.5729 0.33109 9.2038201719_s_at 2.4708 0

6、.4047 0.93808 6.5081202945_at 0.4118 2.4284 0.03990 4.2499203261_at 3.1502 0.3174 0.33641 29.4983203757_s_at 0.7861 1.2720 0.61656 1.0024205068_s_at 0.1073 9.3180 0.02223 0.5181最后選取六個(gè)變量擬合生存模型，繪制生存曲線如下：下面繪制ROC曲線，分別在訓(xùn)練集和測(cè)試集上繪制ROC曲線，結(jié)果如下：訓(xùn)練集：測(cè)試集：由于測(cè)試集上的樣本過(guò)少，所以得到的AUC值波動(dòng)大，考慮使用bootstrap屢次計(jì)算訓(xùn)練集上的AUC值并求平均來(lái)測(cè)

7、試模型的效果：AUC at 1 year：0.8039456AUC at 3 year：0.6956907AUC at 5 year：0.7024846由此可以看到，隨機(jī)森林通過(guò)刪除奉獻(xiàn)較低的變量，完成變量選擇的工作，在測(cè)試集上具有較高的AUC值，但是比lasso-co*模型得到的AUC略低。附錄：load(/R/brea.rda)library(survival)set.seed(10)i-sample(1:77,52)train-dati,test-dat-i,library(randomForestSRC)disease.rf-rfsrc(Surv(time,status).,data

8、= train, ntree = 800,mtry = 3, nodesize = 3,splitrule = logrank)disease.rfour.rf- var.select(object=disease.rf, vdv, method = vh.vimp, nrep = 50)our.rf$topvarsinde*-numeric(var.rf$modelsize)for(i in 1:var.rf$modelsize) inde*i-which(names(dat)=var.rf$topvarsi)data-dat,c(1,2,inde*)i-sample(1:77,52)tra

9、in-datai,test-data-i,mod.brea-co*ph(Surv(time,status).,data=train)train_data-train,c(1,2,which(summary(mod.brea)$coefficients,5=0.1)+2)tset_data-test,c(1,2,which(summary(mod.brea)$coefficients,5=0.1)+2)mod.brea1-co*ph(Surv(time,status).,data=train_data)summary(mod.brea1)names(coef(mod.brea1)plot(sur

10、vfit(mod.brea1),*lab=Time,ylab = Proportion,main=Co* Model,=TRUE,col=c(black,red,red),ylim=c(0.6,1)inde*0-numeric(length(coef(mod.brea1)coefficients-coef(mod.brea1)name-gsub(,names(coefficients)for(j in 1:length(inde*0) inde*0j-which(names(dat)=namej)library(survivalROC)riskscore-as.matri*(d

11、ati,inde*0)%*% as.matri*(coefficients)y1-survivalROC(Stime=train$time,status=train$status,marker=riskscore,predict.time=1,span = 0.25*(nrow(train)(-0.20)y3-survivalROC(Stime=train$time,status=train$status,marker=riskscore,predict.time=3,span = 0.25*(nrow(train)(-0.20)y5-survivalROC(Stime=train$time,

12、status=train$status,marker=riskscore,predict.time=5,span = 0.25*(nrow(train)(-0.20)a-matri*(data=c(y1,y3,y5,y1$AUC,y3$AUC,y5$AUC),nrow=3,ncol=2);aplot(y1$FP,y1$TP,type=l,*lab=False Positive Rate,ylab = True Positive Rate,main=Time-dependent ROC curve,col=green) lines(y3$FP,y3$TP,col=red,lty=2)lines(

13、y5$FP,y5$TP,col=blue,lty=3)legend(bottomright,bty=n,legend = c(AUC at 1 year:0.9271,AUC at 3 years:0.8621,AUC at 5 years:0.8263),col=c(green,red,blue),lty=c(1,2,3),ce*=0.9)abline(0,1)riskscore-as.matri*(dat-i,inde*0)%*% as.matri*(coefficients)y1-survivalROC(Stime=test$time,status=test$status,marker=

14、riskscore,predict.time=1,span = 0.25*(nrow(train)(-0.20)y3-survivalROC(Stime=test$time,status=test$status,marker=riskscore,predict.time=3,span = 0.25*(nrow(train)(-0.20)y5-survivalROC(Stime=test$time,status=test$status,marker=riskscore,predict.time=5,span = 0.25*(nrow(train)(-0.20)a-matri*(data=c(y1

15、,y3,y5,y1$AUC,y3$AUC,y5$AUC),nrow=3,ncol=2);aplot(y1$FP,y1$TP,type=l,*lab=False Positive Rate,ylab = True Positive Rate,main=Time-dependent ROC curve,col=green) lines(y3$FP,y3$TP,col=red,lty=2)lines(y5$FP,y5$TP,col=blue,lty=3)legend(bottomright,bty=n,legend = c(AUC at 1 year:0.8761,AUC at 3 years:0.

16、7611,AUC at 5 years:0.7611),col=c(green,red,blue),lty=c(1,2,3),ce*=0.9)abline(0,1)a-matri*(0,30,3)for (c in 1:30) i-sample(1:77,52) train-datai, test-data-i, mod.brea-co*ph(Surv(time,status).,data=train) train_data-train,c(1,2,which(summary(mod.brea)$coefficients,5=0.1)+2) tset_data-test,c(1,2,which

17、(summary(mod.brea)$coefficients,5=0.1)+2) mod.brea1-co*ph(Surv(time,status).,data=train_data) names(coef(mod.brea1) inde*0-numeric(length(coef(mod.brea1) coefficients-coef(mod.brea1) name-gsub(,names(coefficients) for(j in 1:length(inde*0) inde*0j-which(names(dat)=namej) riskscore-as.matri*(dat-i,inde*0)%*% as.matri*(coefficients) y1-survivalROC(Stime=test$time,status=test$status,marker=riskscore,predict.time=1,span = 0.25*