2024年4月3日发(作者：)

随机森林模型在生存分析中的应用

【摘要】目的：本文探讨随机森林方法用于高维度、强相关、小样本的生

存资料分析时，可以起到变量筛选的作用。方法：以乳腺癌数据集构建乳腺癌转

移风险评估模型为实例进行实证分析，使用随机森林模型进行变量选择，然后拟

合cox回归模型。结果：随机森林模型通过对变量的选择，有效的解决数据维

度高且强相关的情况，得到了较高的AUC值。

一、数据说明

该乳腺癌数据集来自于NCBI，有77个观测值以及22286个基因变量。通过

筛选选取454个基因变量。将数据随机分为训练集合测试集，其中2/3为训练集，

1/3为测试集。绘制K-M曲线图：

二、随机森林模型

随机森林由许多的决策树组成，因为这些决策树的形成采用了随机的方法，

因此也叫做随机决策树。随机森林中的树之间是没有关联的。当测试数据进入随

机森林时，其实就是让每一颗决策树进行分类，最后取所有决策树中分类结果最

多的那类为最终的结果。因此随机森林是一个包含多个决策树的分类器，并且其

输出的类别是由个别树输出的类别的众数而定。

使用randomForestSRC包得到的随机森林模型具有以下性质：

Numberofdeaths:27

Numberoftrees:800

Minimumterminalnodesize:3

inalnodes:14.4275

ablestriedateachsplit:3

ables:452

Analysis:RSF

Family:surv

Splittingrule:logrank

Errorrate:19.87%

发现直接使用随机森林得到的模型，预测误差很大，达到了19.8%,

进一步考虑使用随机森林模型进行变量选择，结果如下：

>$

Samplesize:52

Numberofdeaths:19

Numberoftrees:500

Minimumterminalnodesize:2

inalnodes:11.554

ablestriedateachsplit:3

ables:9

Analysis:RSF

Family:surv

Splittingrule:logrank*random*

Numberofrandomsplitpoints:10

Errorrate:11.4%

>$topvars

[1]"213821_s_at""219778_at"

at"

[6]"211603_s_at""213055_at""219336_s_at""37892_at"

"204690_at""220788_s_at""202202_s_

一共选取了

9

个变量，同时误差只有

11.4%

接下来，使用这些变量做

cox

回归，剔除模型中不显著（>0.01）的变量，最终

参与模型建立的变量共有4个。模型结果如下：

exp(coef)exp(-coef)lower.95upper.95

`218150_at`

`200914_x_at`

`220788_s_at`

`201398_s_at`

`201719_s_at`

`202945_at`

`203261_at`

`203757_s_at`

`205068_s_at`

1.6541

0.9915

0.2649

1.7457

2.4708

0.4118

3.1502

0.7861

0.1073

0.6046

1.0086

3.7750

0.5729

0.4047

2.4284

0.3174

1.2720

9.3180

0.11086

0.34094

0.05944

0.33109

0.93808

0.03990

0.33641

0.61656

0.02223

24.6800

2.8833

1.1805

9.2038

6.5081

4.2499

29.4983

1.0024

0.5181

最后选取六个变量拟合生存模型，绘制生存曲线如下：

下面绘制ROC曲线，分别在训练集和测试集上绘制ROC曲线，结果如下：

训练集：

测试集：

由于测试集上的样本过少，所以得到的

AUC

值波动大，考虑使用

bootstrap

多次计算训练集上的

AUC

值并求平均来测试模型的效果：

AUCat1year

：

0.8039456

AUCat3year

：

0.6956907

AUCat5year

：

0.7024846

由此可以看到，随机森林通过删除贡献较低的变量，完成变量选择的工作，在测

试集上具有较高的

AUC

值，但是比

lasso-cox

模型得到的

AUC

略低。

附录：

load("~/R/")

library(survival)

(10)

i<-sample(1:77,52)

train<-dat[i,]

test<-dat[-i,]

library(randomForestSRC)

<-rfsrc(Surv(time,status)~.,data=train,

ntree=800,mtry=3,

nodesize=3,splitrule="logrank")

<-(object=,vdv,

method="",nrep=50)

$

$topvars

index<-numeric($modelsize)

for(iin1:$modelsize){

index[i]<-which(names(dat)==$topvars[i])

}

data<-dat[,c(1,2,index)]

i<-sample(1:77,52)

train<-data[i,]

test<-data[-i,]

<-coxph(Surv(time,status)~.,data=train)

train_data<-train[,c(1,2,which(summary()$coefficients[,5]<=0.1)+2)]

tset_data<-test[,c(1,2,which(summary()$coefficients[,5]<=0.1)+2)]

1<-coxph(Surv(time,status)~.,data=train_data)

summary(1)

names(coef(1))

plot(survfit(1),xlab="Time",ylab="Proportion",main="Cox

Model",=TRUE,col=c("black","red","red"),ylim=c(0.6,1))

index0<-numeric(length(coef(1)))

coefficients<-coef(1)

name<-gsub("`","",names(coefficients))

for(jin1:length(index0)){

index0[j]<-which(names(dat)==name[j])

}

library(survivalROC)

riskscore<-(dat[i,index0])%*%(coefficients)

y1<-survivalROC(Stime=train$time,status=train$status,marker=riskscore,

=1,span=0.25*(nrow(train))^(-0.20))

y3<-survivalROC(Stime=train$time,status=train$status,marker=riskscore,

=3,span=0.25*(nrow(train))^(-0.20))

y5<-survivalROC(Stime=train$time,status=train$status,marker=riskscore,

=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);a

plot(y1$FP,y1$TP,type="l",xlab="FalsePositiveRate",ylab="TruePositive

Rate",main="Time-dependentROCcurve",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("AUCat1year:0.9271","AUCat3

years:0.8621","AUCat5

years:0.8263"),col=c("green","red","blue"),lty=c(1,2,3),cex=0.9)

abline(0,1)

riskscore<-(dat[-i,index0])%*%(coefficients)

y1<-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=

1,span=0.25*(nrow(train))^(-0.20))

y3<-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=

3,span=0.25*(nrow(train))^(-0.20))

y5<-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=

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);a

plot(y1$FP,y1$TP,type="l",xlab="FalsePositiveRate",ylab="TruePositive

Rate",main="Time-dependentROCcurve",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("AUCat1year:0.8761","AUCat3

years:0.7611","AUCat5

years:0.7611"),col=c("green","red","blue"),lty=c(1,2,3),cex=0.9)

abline(0,1)

a<-matrix(0,30,3)

for(cin1:30){

i<-sample(1:77,52)

train<-data[i,]

test<-data[-i,]

<-coxph(Surv(time,status)~.,data=train)

train_data<-train[,c(1,2,which(summary()$coefficients[,5]<=0.1)+2)]

tset_data<-test[,c(1,2,which(summary()$coefficients[,5]<=0.1)+2)]

1<-coxph(Surv(time,status)~.,data=train_data)

names(coef(1))

index0<-numeric(length(coef(1)))

coefficients<-coef(1)

name<-gsub("`","",names(coefficients))

for(jin1:length(index0)){

index0[j]<-which(names(dat)==name[j])

}

riskscore<-(dat[-i,index0])%*%(coefficients)

y1<-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=

1,span=0.25*(nrow(train))^(-0.20))

y3<-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=

3,span=0.25*(nrow(train))^(-0.20))

y5<-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=

5,span=0.25*(nrow(train))^(-0.20))

a[c,]<-c(y1$AUC,y3$AUC,y5$AUC)

}