生存分析随机森林实验与代码文档格式.docx
《生存分析随机森林实验与代码文档格式.docx》由会员分享,可在线阅读,更多相关《生存分析随机森林实验与代码文档格式.docx(8页珍藏版)》请在冰点文库上搜索。
27
Numberoftrees:
800
Minimumterminalnodesize:
3
Averageno.ofterminalnodes:
No.ofvariablestriedateachsplit:
Totalno.ofvariables:
452
Analysis:
RSF
Family:
surv
Splittingrule:
logrank
Errorrate:
%
发现直接使用随机森林得到的模型,预测误差很大,达到了%,进一步考虑使用随机森林模型进行变量选择,结果如下:
>
$
Samplesize:
52
Numberofdeaths:
19
500
2
9
logrank*random*
Numberofrandomsplitpoints:
10
$topvars
[1]"
213821_s_at"
"
219778_at"
204690_at"
220788_s_at"
202202_s_at"
[6]"
211603_s_at"
213055_at"
219336_s_at"
37892_at"
一共选取了9个变量,同时误差只有%
接下来,使用这些变量做cox回归,剔除模型中不显著(>
)的变量,最终参与模型建立的变量共有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`
最后选取六个变量拟合生存模型,绘制生存曲线如下:
下面绘制ROC曲线,分别在训练集和测试集上绘制ROC曲线,结果如下:
训练集:
测试集:
由于测试集上的样本过少,所以得到的AUC值波动大,考虑使用bootstrap多次计算训练集上的AUC值并求平均来测试模型的效果:
AUCat1year:
AUCat3year:
AUCat5year:
由此可以看到,随机森林通过删除贡献较低的变量,完成变量选择的工作,在测试集上具有较高的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)
$in1:
$modelsize){
index[i]<
-which(names(dat)==$topvars[i])
}
data<
-dat[,c(1,2,index)]
-data[i,]
-data[-i,]
-coxph(Surv(time,status)~.,data=train)
train_data<
-train[,c(1,2,which(summary$coefficients[,5]<
=+2)]
tset_data<
-test[,c(1,2,which(summary$coefficients[,5]<
-coxph(Surv(time,status)~.,data=train_data)
summary
names(coef)
plot(survfit,xlab="
Time"
ylab="
Proportion"
main="
CoxModel"
=TRUE,col=c("
black"
"
red"
),ylim=c,1))
index0<
-numeric(length(coef))
coefficients<
-coef
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=*(nrow(train))^)
y3<
-survivalROC(Stime=train$time,status=train$status,marker=riskscore,=3,span=*(nrow(train))^)
y5<
-survivalROC(Stime=train$time,status=train$status,marker=riskscore,=5,span=*(nrow(train))^)
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"
TruePositiveRate"
Time-dependentROCcurve"
col="
green"
)
lines(y3$FP,y3$TP,col="
lty=2)
lines(y5$FP,y5$TP,col="
blue"
lty=3)
legend("
bottomright"
bty="
n"
legend=c("
AUCat1year:
AUCat3years:
AUCat5years:
),col=c("
),lty=c(1,2,3),cex=
abline(0,1)
(dat[-i,index0])%*%(coefficients)
-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=1,span=*(nrow(train))^)
-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=3,span=*(nrow(train))^)
-survivalROC(Stime=test$time,status=test$status,marker=riskscore,=5,span=*(nrow(train))^)
-matrix(0,30,3)
for(cin1:
30){
i<
train<
test<
<
train_data<
tset_data<
names(coef)
index0<
coefficients<
name<
for(jin1:
}
riskscore<
y1<
y3<
y5<
a[c,]<
-c(y1$AUC,y3$AUC,y5$AUC)