利用随机，提高timeROC

叉叉滴同学 叉叉滴同学的生信笔记 2023-06-29 07:00 发表于广西

之前给同学们分享了利用随机提高ROC的AUC值

不过那个不是的timeROC的AUC

👉利用随机提高ROC曲线的AUC

今天就给同学们分享利用随机，找到timeROC的最优AUC值

废话不多说，直接整👇👇

1.文件准备

准备生存时间、生存状态及基因表达数据的合并文件

图片

2.加载包

library(pROC)library(limma)library(tidyverse)library(future.apply)library(neuralnet)library(caret)library(NeuralNetTools)library(caTools)library(survival)library("glmnet")library(survival)library(survminer)library(timeROC)library(limma)library(future.apply)library(parallel)library(caret)

3.读取文件，，并将生存时间单位转换为年

#读取输入文件rt=read.table(file ="exptime.txt" , header=T, sep="\t", check.names=F,row.names = 1)data=as.data.frame(rt)data$futime=data$futime/365

图片

4.设置主要参数及阈值

#设置随机次数counts=1000#设置训练集与测试集比例pr=0.6   #【0.6就是训练集占比60%】#设置出图阈值plotedge=0.4samplecounts=as.data.frame(1:counts)#【输入你要研究的基因】gene="CES1"

出图的AUC阈值，为了方便演示在这里我设置为了0.4【多数情况下应设置为0.75或更高】

5.进入多线程循环

#开启多线程循环plan(multisession)future_apply(X = samplecounts,MARGIN = 1,FUN = function(seed){  library(pROC)  library(limma)  library(tidyverse)  library(future.apply)  library(neuralnet)  library(caret)  library(NeuralNetTools)  library(caTools)  library(survival)  library("glmnet")  library(survival)  library(survminer)  library(timeROC)  library(limma)  library(future.apply)  library(parallel)  library(caret)

由于利用future_apply之后很多包需要重新加载，因此直接将所有的包放在了future_apply之后

6.将矩阵随机分为train组和test组

split_T <- sample.split(rownames(data), SplitRatio = pr)trainData <- subset(data, split_T =="TRUE")%>%as.data.frame()testData <- subset(data, split_T =="FALSE")%>%as.data.frame()trainData_gene=trainData[,gene,drop=F]trainData_gene=data.frame(futime=trainData$futime,fustat=trainData$fustat,gene=trainData_gene[,1])rownames(trainData_gene)=rownames(trainData)testData_gene=testData[,gene,drop=F]testData_gene=data.frame(futime=testData$futime,fustat=testData$fustat,gene=testData_gene[,1])rownames(testData_gene)=rownames(testData)

图片

图片

可以看到，train组数据随机分到了200个样本，test组数据随机分到了134个样本

7.进行timeROC分析

##【计算timeROC】ROC_train=timeROC(T=trainData_gene$futime,delta=trainData_gene$fustat,                  marker=trainData_gene$gene,cause=1,                  weighting='aalen',                  times=c(1,3,5),ROC=TRUE)AUC_1_year_train=ROC_train[["AUC"]][["t=1"]]AUC_3_year_train=ROC_train[["AUC"]][["t=3"]]AUC_5_year_train=ROC_train[["AUC"]][["t=5"]]ROC_test=timeROC(T=testData_gene$futime,delta=testData_gene$fustat,                       marker=testData_gene$gene,cause=1,                       weighting='aalen',                       times=c(1,3,5),ROC=TRUE)AUC_1_year_test=ROC_test[["AUC"]][["t=1"]]AUC_3_year_test=ROC_test[["AUC"]][["t=3"]]AUC_5_year_test=ROC_test[["AUC"]][["t=5"]]

这里选择了1，3，5年进行分析，同学们也可以修改为其他年份进行分析

8.判断是否符合阈值并绘制图片

这么多次的循环，每次都输出图片的话岂不是得筛选到天荒地老

因此我们需要对其添加限制条件👇

train组和test组1，3，5年的AUC均需大于阈值才可输出

##【判断条件】if(AUC_1_year_train>plotedge &    AUC_3_year_train>plotedge &    AUC_5_year_train>plotedge &    AUC_1_year_test>plotedge &    AUC_3_year_test>plotedge &    AUC_5_year_test>plotedge){  pdf(file=paste0(seed,"_Train_ROC.pdf"), width=5, height=5)  plot(ROC_train,time=1,col="green",title=FALSE,lwd=2)  plot(ROC_train,time=3,col="blue",add=TRUE,title=FALSE,lwd=2)  plot(ROC_train,time=5,col="red",add=TRUE,title=FALSE,lwd=2)  legend('bottomright',         c(paste0('AUC at 1 years: ',sprintf("%.03f",ROC_train$AUC[1])),           paste0('AUC at 3 years: ',sprintf("%.03f",ROC_train$AUC[2])),           paste0('AUC at 5 years: ',sprintf("%.03f",ROC_train$AUC[3]))),         col=c("green","blue","red"), lwd=2, bty = 'n')  dev.off()  pdf(file=paste0(seed,"_Test_ROC.pdf"), width=5, height=5)  plot(ROC_test,time=1,col="green",title=FALSE,lwd=2)  plot(ROC_test,time=3,col="blue",add=TRUE,title=FALSE,lwd=2)  plot(ROC_test,time=5,col="red",add=TRUE,title=FALSE,lwd=2)  legend('bottomright',         c(paste0('AUC at 1 years: ',sprintf("%.03f",ROC_test$AUC[1])),           paste0('AUC at 3 years: ',sprintf("%.03f",ROC_test$AUC[2])),           paste0('AUC at 5 years: ',sprintf("%.03f",ROC_test$AUC[3]))),         col=c("green","blue","red"), lwd=2, bty = 'n')  dev.off()  }})

这里为了方便演示，阈值设置为了0.4，因此输出了很多图片

同学们在实际操作中可以适当的提高阈值

图片

图片

代码如上所示

但，如果是比较大的矩阵，想要明确train和test的最佳比例

到底需要循环多少次？？？

这也是一个非常常见且棘手的问题

解决的方式也很简单，就是让代码不停地跑——无人值守

👇设置阈值范围，让代码自己跑

#设置出图阈值#【为了演示直观，所以阈值调到了0.3、0.4】plotedges=c(0.3,0.4,0.5,0.6,0.65,0.9)#【设置阈值，建议多设几个阶梯，如果没跑过循环建议从0.5开始设，梯度尽量小】#【输入你要研究的基因】genes="CES1"workspace=getwd(){source("主代码_train.R",encoding = "UTF-8")job::job(train={XXD_timeROC(datainput,plotedges,genes,workspace)})}

此处，同学们可以依照之前跑单次timeROC时候的经验进行调整。

如果一次都没跑过的话，可以从0.5开始设置一个阈值，随后以0.03左右为一个梯度进行增加

若得到结果，则会构建文件夹并储存图片，文件夹上则显示了关键的信息——阈值及train组的比例

图片

若卡在了某个阈值，则代码会进行不断循环，直到AUC符合阈值便进行输出【循环到符合为止，无需手动重新设置循环】

此外，无人值守为后台运行代码，即使将Rstudio出现bug意外关闭，也不影响代码的运行，它会在后台完成它的使命图片。

同样的，由于代码运行于后台，因此不会像其它代码一般“霸占”Rstudio，同学们仍然可以进行其它分析。尽可能提高效率

利用随机提高timeROC的代码就给同学们分享到这

完整的随机代码如上所示

无人值守的代码可在下方付费获取哦

👇感谢同学们的关注与支持👇

1元 = 10微信豆

可试读95%

1元 = 10微信豆
喜欢此内容的人还喜欢

微信扫一扫
关注该公众号