我们在看临床模型类文献的时候,虽然常看到用 X-tile 寻找变量的最佳cutoff值,但是有时候也会看到有的文章是用ROC曲线来寻找最佳cutoff值的,下面我们一探究竟吧,注本期所用的连续型变量为riskscore,而构建riskscore的基因表达量也是连续型变量。
载入必要的R包
rm(list=ls())
library(tibble)
library(dplyr)
library(tidyr)
library(survivalROC)
载入数据
clinical<-data.table::fread(file ="risk.txt",data.table = F)
head(clinical)[1:5,1:5]
## 将第一列变成行名
clinical <- clinical %>%
column_to_rownames("sample")
处理时间
### 将生存时间转化成年(之前已经为年了,所以除以1。如果为天则除以365。月则除以12)
clinical$futime=clinical$futime/1
predict_time1=3
predict_time2=5
绘制3年的ROC曲线
Auc_Text=c()
you_roc=survivalROC(Stime=clinical$futime, status=clinical$fustat, marker = clinical$riskScore, predict.time =predict_time1, method="KM")
plot(you_roc$FP, you_roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col="#DAA520",
xlab="False positive rate",
ylab="True positive rate",
lwd = 3, cex.main=1.5, cex.lab=1.3, cex.axis=1.3, font=1.3)
Auc_Text=c(Auc_Text,paste0("riskScore"," (AUC=",sprintf("%.3f",you_roc$AUC),")"))
abline(0,1)
legend("bottomright",Auc_Text,lwd=2,bty="n",col="#DAA520")
#dev.off()
### 利用ROC曲线计算最佳cutoff
cutoff_3year<-you_roc$cut.values[which.max(you_roc$TP-you_roc$FP)]
cutoff_3year
绘制5年的ROC曲线
Auc_Text=c()
you_roc=survivalROC(Stime=clinical$futime, status=clinical$fustat, marker = clinical$riskScore, predict.time =predict_time2, method="KM")
plot(you_roc$FP, you_roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col="#EE3B3B",
xlab="False positive rate",
ylab="True positive rate",
lwd = 3, cex.main=1.5, cex.lab=1.3, cex.axis=1.3, font=1.3)
Auc_Text=c(Auc_Text,paste0("riskScore"," (AUC=",sprintf("%.3f",you_roc$AUC),")"))
abline(0,1)
legend("bottomright",Auc_Text,lwd=2,bty="n",col="#EE3B3B")
#dev.off()
#利用ROC曲线计算最佳cutoff
cutoff_5year<-you_roc$cut.values[which.max(you_roc$TP-you_roc$FP)]
cutoff_5year
那么当构建riskscore的元素为分期、性别、亚型等分类变量时,又应该怎么用ROC曲线寻找riskscore的最佳cutoff呢?且听下回分解。
更多详细可见:https://mp.weixin.qq.com/s/0wZTazrIl2EQ1__9gNdywQ