生存分析-随机森林实验与代码.docx

1、.随机森林模型在生存分析中的应用【摘要】 目的：本文探讨随机森林方法用于高维度、强相关、小样本的生存资料分析时， 可以起到变量筛选的作用。 方法：以乳腺癌数据集构建乳腺癌转移风险评估模型为实例进行实证分析， 使用随机森林模型进行变量选择， 然后拟合 cox 回归模型。 结果： 随机森林模型通过对变量的选择，有效的解决数据维度高且强相关的情况，得到了较高的 AUC 值。一、数据说明该乳腺癌数据集来自于NCBI，有 77 个观测值以及 22286 个基因变量。通过筛选选取 454 个基因变量。将数据随机分为训练集合测试集， 其中 2/3 为训练集，1/3 为测试集。绘制 K-M 曲线图：二、随机森

2、林模型随机森林由许多的决策树组成，因为这些决策树的形成采用了随机的方法，因此也叫做随机决策树。 随机森林中的树之间是没有关联的。当测试数据进入随机森林时，其实就是让每一颗决策树进行分类，最后取所有决策树中分类结果最多的那类为最终的结果。 因此随机森林是一个包含多个决策树的分类器，并且其输出的类别是由个别树输出的类别的众数而定。.使用 randomForestSRC 包得到的随机森林模型具有以下性质：Numberof deaths:27Numberof trees: 800Minimumterminalnodesize: 3Averageno. of terminalnodes:14.4275N

3、o.of variablestriedat eachsplit: 3Totalno. of variables:452Analysis: RSFFamily: survSplittingrule: logrankError rate: 19.87%发现直接使用随机森林得到的模型，预测误差很大，达到了19.8%,进一步考虑使用随机森林模型进行变量选择，结果如下：our.rf$rfsrc.refit.objSamplesize: 52Numberof deaths:19Numberof trees: 500Minimumterminalnodesize: 2Averageno. of termin

4、alnodes:11.554No.of variablestriedat eachsplit: 3Totalno. of variables:9Analysis: RSFFamily: survSplittingrule: logrank*random*Numberof randomsplitpoints:10Error rate: 11.4%our.rf$topvars1 213821_s_at219778_at204690_at220788_s_at202202_s_at6 211603_s_at213055_at219336_s_at37892_at一共选取了 9 个变量，同时误差只有1

5、1.4%接下来，使用这些变量做cox 回归，剔除模型中不显著（0.01 ）的变量，最终参与模型建立的变量共有4 个。模型结果如下：exp(coef)exp(-coef)lower .95upper .95218150_at1.65410.60460.1108624.6800200914_x_at0.99151.00860.340942.8833220788_s_at0.26493.77500.059441.1805201398_s_at1.74570.57290.331099.2038201719_s_at2.47080.40470.938086.5081.202945_at0.41182.4

6、2840.039904.2499203261_at3.15020.31740.3364129.4983203757_s_at0.78611.27200.616561.0024205068_s_at0.10739.31800.022230.5181最后选取六个变量拟合生存模型，绘制生存曲线如下：下面绘制 ROC 曲线，分别在训练集和测试集上绘制ROC 曲线，结果如下：训练集：测试集：.由于测试集上的样本过少， 所以得到的 AUC 值波动大，考虑使用 bootstrap 多次计算训练集上的AUC 值并求平均来测试模型的效果：AUC at 1 year：0.8039456AUC at 3 year：

7、0.6956907AUC at 5 year：0.7024846由此可以看到， 随机森林通过删除贡献较低的变量，完成变量选择的工作， 在测试集上具有较高的AUC 值，但是比 lasso-cox 模型得到的 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 = train,ntree = 800,mtry = 3,nodesiz

8、e = 3,splitrule = logrank)disease.rfour.rf- var.select(object=disease.rf, vdv,method = vh.vimp, nrep = 50)our.rf$rfsrc.refit.objour.rf$topvarsindex-numeric(var.rf$modelsize)for(i in 1:var.rf$modelsize)indexi-which(names(dat)=var.rf$topvarsi)data-dat,c(1,2,index)i-sample(1:77,52)train-datai,test-data

9、-i,mod.brea-coxph(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-coxph(Surv(time,status).,data=train_data)summary(mod.brea1)names(coef(mod.brea1)plot(survfit(mod.brea1),x

10、lab=Time,ylab=Proportion,main=CoxModel,=TRUE,col=c(black,red,red),ylim=c(0.6,1)index0-numeric(length(coef(mod.brea1)coefficients-coef(mod.brea1)name-gsub(,names(coefficients)for(j in 1:length(index0)index0j-which(names(dat)=namej)library(survivalROC)riskscore-as.matrix(dati,index0)%*% as.mat

11、rix(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,status=train$status,m

12、arker=riskscore,predict.time=5,span = 0.25*(nrow(train)(-0.20)a-matrix(data=c(y1,y3,y5,y1$AUC,y3$AUC,y5$AUC),nrow=3,ncol=2);aplot(y1$FP,y1$TP,type=l,xlab=FalsePositiveRate,ylab=TruePositiveRate,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)l

13、egend(bottomright,bty=n,legend=c(AUC at 1 year:0.9271,AUC at 3years:0.8621,AUCat5years:0.8263),col=c(green,red,blue),lty=c(1,2,3),cex=0.9)abline(0,1)riskscore-as.matrix(dat-i,index0)%*% as.matrix(coefficients)y1-survivalROC(Stime=test$time,status=test$status,marker=riskscore,predict.time=1,span = 0.

14、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-matrix(data=c(y1,y3,y5,y1$AUC,y3$AUC,y5$AUC),nrow=

15、3,ncol=2);a.plot(y1$FP,y1$TP,type=l,xlab=FalsePositiveRate,ylab=TruePositiveRate,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 3years:0.7611,AUCat5years:0.7611),col=c(green,red,bl

16、ue),lty=c(1,2,3),cex=0.9)abline(0,1)a-matrix(0,30,3)for (c in 1:30)i-sample(1:77,52)train-datai,test-data-i,mod.brea-coxph(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-c

17、oxph(Surv(time,status).,data=train_data)names(coef(mod.brea1)index0-numeric(length(coef(mod.brea1).coefficients-coef(mod.brea1)name-gsub(,names(coefficients)for(j in 1:length(index0)index0j-which(names(dat)=namej)riskscore-as.matrix(dat-i,index0)%*% as.matrix(coefficients)y1-survivalROC(Stime=test$time,status=test$status,marker=riskscore,predict.time=1,span = 0.25*(