TCGA+biomarker——ROC曲线
通常情况下，通过以下几种指标来对模型进行评价。
1）区分度：采用指标C-index和ROC曲线来评价区分度，一般文章都是二选一。
- C-index即一致性指数（index of concordance），通过评估模型预测结果与实际观察结果的符合程度，以评价模型的预测准确性。
- ROC曲线，展示特异性和敏感性，ROC曲线下的面积被称为AUC，它介于0.5和1之间，作为数值可以直观的评价分类器的好坏，值越大越好。
2）校准度：通常采用校准曲线（calibration curve）来进行评价一致性/校准度，即预测值和真实值之间的差异
3）DCA：决策曲线（DCA）用来帮助确定高风险的患者进行干预、低风险的患者避免过度医疗。
前期已经介绍过C-index，Calibration curve，DCA，今天再来说说ROC曲线......
ROC曲线简介
Receiver Operating Characteristic Curve，简称ROC曲线，又称为感受性曲线。它是根据一系列不同的二分类方式（分界值或决定阈），以真阳性率（敏感性）为纵坐标，假阳性率（1-特异性）为横坐标绘制的曲线。AUC（Area Under Curve）被定义为ROC曲线下的面积。我们往往使用AUC值作为模型的评价标准是因为很多时候ROC曲线并不能清晰的说明哪个分类器的效果更好，而作为一个数值，对应AUC更大的分类器效果更好。
ROC曲线案例
- 展示不同时间的ROC曲线
image
- 展示不同模型的ROC曲线
image
-
不同数据集的ROC曲线
image
ROC曲线图解
image
解释：ROC曲线图是反映敏感性与特异性之间关系的曲线。横坐标X轴为特异性，也称为假阳性率（误报率：1-Specificity ），X轴越接近零准确率越高；纵坐标Y轴称为敏感度，也称为真阳性率（敏感度：Sensitivity），Y轴越大代表准确率越好。通过评估不同阈值的真阳性和假阳性，可以构建一条曲线，该曲线从左下方延伸到右上方，并向左上方弯曲。在正、负类之间没有判别力的分类器将形成一条对角线，在假阳性率0和真阳性率0（坐标：0,0）与假阳性率 1和真阳性率（坐标：1,1）。
根据曲线位置，把整个图划分成了两部分，曲线下方部分的面积被称为AUC（Area Under Curve），用来表示预测准确性，AUC值越高，也就是曲线下方面积越大，说明预测准确率越高。曲线越接近左上角（X越小，Y越大），预测准确率越高。ROC曲线下的面积值AUC正常在0.5和1之间。在AUC>0.5的情况下，AUC越接近于1，说明诊断效果越好。AUC在 0.5～0.7时有较低准确性，AUC在0.7～0.9时有一定准确性，AUC在0.9以上时有较高准确性。AUC=0.5时，说明诊断方法完全不起作用，无诊断价值。AUC<0.5不符合真实情况，在实际中极少出现。
ROC曲线绘制原理
为了更好理解ROC曲线绘图思路，我们需要了解X轴和Y轴数据是如何计算的？
“混淆矩阵”：对于二分类问题，可将样本根据其真实类别与学习器预测类别的组合划分为TP(true positive)、FP(false positive)、TN(true negative)、FN(false negative)四种情况，TP+FP+TN+FN=样本总数。
[图片上传中...(image-3930ac-1594171865147-3)]
(1) 真阳性(True Positive，TP)：检测阳性，且实际阳性；正确肯定的匹配数目；
(2) 假阳性(False Positive，FP)：检测阳性，但实际阴性；误报，给出的匹配是不正确的；
(3) 真阴性(True Negative，TN)：检测阴性，且实际阴性；正确拒绝的非匹配数目；
(4) 假阴性(False Negative，FN)：检测阴性，但实际阳性；漏报，没有正确找到的匹配的数目。
X轴—假阳率（False Positive Rate）的计算方法是：假阳（False Positive）预测的总数除以假阳和真阴（True Negative）的总和。False Positive Rate = False Positives /（False Positives + True Negatives）
Y轴—真阳率（True Positive Rate）是用真阳性（True Positive）预测的总数除以真阳和假阴（False Negative）之和得出的分数。真阳性率称为敏感性或召回率。True Positive Rate = True Positives /（True Positives +False Negatives ）
下面以一个具体的例子来详细了解ROC曲线是如何绘制的。
image
如上图所示，共计10例样本，根据event分成两类，event=1的表示阳性样本，event=0的表示阴性样本，Riskscore表示测量的某个指标值。
1）首先，依据Riskscore值从大到小对这10例样本进行排序（上图已经是按此规则排过序的）；
2）接下来，依次把Riskscore值作为阈值（即阈值依次为21.34，21.31，19.42.........11.33）,共计10个阈值，当Riskscore值大于等于此阈值时被认为是阳性，否则被认为是阴性。比如，当以Riskscore=17.17作为阈值，那么前4例样本被分类成阳性样本，剩余6例样本被分类成阴性样本。
据此，我们可以得到：
真阳性(True Positive，TP)：检测阳性，且实际阳性；TP=2
假阳性(False Positive，FP)：检测阳性，但实际阴性；FP=2
真阴性(True Negative，TN)：检测阴性，且实际阴性；TN=2
假阴性(False Negative，FN)：检测阴性，但实际阳性；FN=4
那么，根据以上指标可计算：
FPR = FP/（FP+ TN）=2/2+2=0.50
TPR = TP /（TP +FN）=2/2+4=0.33
这样，依次把这10个Riskscore值作为阈值，我们就可以得到10组（TPR,FPR）值，把这10组（TPR,FPR）绘制出来得到的曲线就是ROC曲线。
ROC曲线的目的
目的1——评估某个模型的好坏！
在AUC>0.5的情况下，AUC越接近于1，说明诊断效果越好。AUC在 0.5～0.7时有较低准确性，AUC在0.7～0.9时有一定准确性，AUC在0.9以上时有较高准确性。下图表明AUC=0.79, 说明该模型分类器的准确性处于中等水平。
image
目的2——比较不同模型的优劣
下图中有3个模型的ROC曲线，根据AUC面积大小，可直观展示出3个模型之间的优劣，A模型>B模型>C模型。
image
但如下图，上图中两条ROC曲线相交于一点，AUC值几乎一样：当需要高Sensitivity时，模型A比B好；当需要高Speciticity时，模型B比A好。
image
目的3——确定最佳临界点
所谓找到最优临界点，就是保证真阳性率（敏感高）高的同时假阳性率（特异性）要尽量的小，建立max（TPR+（1-FPR））的模型。例如：交叉点坐标为（95.8%，14.6%），相应的临界值即可出来，说明敏感度是95.8%，特异度是84.6%。
image
如何绘制ROC曲线？
以下代码源自生信技能树公共号中“学徒带你7步3251行代码+300行注释完成TCGA数据库挖掘实战全文复现”这篇教程，写的非常清楚，采用了三种方式绘制ROC曲线，每一种都很美观！在此基础代码上，也可以灵活运用，绘制不同模型、不同数据集的ROC曲线。
rm(list = ls())
options(stringsAsFactors = F)
#时间依赖的ROC曲线
#载入数据
KM.input<-read.csv(file = "KM.input.csv",header = T)
head(KM.input)
# X event time_year RiskScore
# 1 TCGA-A1-A0SE-01A-11R-A085-13 0 3.6191781 0.3686996
# 2 TCGA-A1-A0SH-01A-11R-A085-13 0 3.9369863 1.8872376
# 3 TCGA-A1-A0SJ-01A-11R-A085-13 0 1.1397260 0.9330429
# 4 TCGA-A1-A0SK-01A-12R-A085-13 1 2.6493151 0.5337272
# 5 TCGA-A1-A0SM-01A-11R-A085-13 0 0.6630137 3.0229704
# 6 TCGA-A1-A0SO-01A-22R-A085-13 0 2.3342466 0.8083318
#数据包括样本名称、事件(生：0,死：1),生存时间，风险值。该数据源自cox风险比例模型
#--------------------------------------survivalROC包---------------------------------------------------
#1.载入R包
#install.packages("survivalROC")
#install.packages("timeROC") #2个包都可以绘制生存时间依赖的ROC曲线
#1.使用SurvivalROC （不能显示置信区间和SD）
library(survival)
library(survivalROC)
?survivalROC #查看这个函数的格式
#输入数据
Survival_ROC_input<-KM.input
#这里我预测五年的生存率
survival_ROC<-survivalROC(Stime=Survival_ROC_input$time_year, #生存时间，Event time or censoring time for subjects
status=Survival_ROC_input$event, #生存状态,dead or alive
marker=Survival_ROC_input$RiskScore, #风险得分，Predictor or marker value
predict.time=5, #预测5年的生存时间
method="KM" #使用KM法进行拟合，默认的方法是method="NNE"
)
survival_ROC_AUC<-round(survival_ROC$AUC,3)#ROC曲线的AUC保留3位小数（文章保留了3位）
#画图
#x轴为False positive rate，y轴为True positive rate
plot(survival_ROC$FP,survival_ROC$TP,type="l",xlim=c(0,1),ylim=c(0,1),
xlab="False positive rate",
ylab="True positive rate",
main=paste0("ROC Curve", " (", "AUC = ",survival_ROC_AUC," )"), #标题
cex.main=1.5,#放大标题
cex.lab=1.3,#坐标轴标签（名称）的缩放倍数
cex.axis=1.3, font=1.2, #放大轴上的数字
lwd=1.5, #指定线条宽度
col="red" #红色
)
abline(a=0,b=1) #添加一条y=x的线
#计算最佳cutoff
cutoff<-survival_ROC$cut.values[which.max(survival_ROC$TP-survival_ROC$FP)]
cutoff
#---------------------------------------timeROC包------------------------------------------------------------------------
#2.使用timeROC(可以计算置信区间和SD）
#Time ROC可以时计算多个时间的AUC
#文章没有说明计算的几年的，我这里计算3，5，10年
library(timeROC)
library(survival)
?timeROC #看一下说明书
#输入数据
time_ROC_input<-KM.input
time_ROC<-timeROC(T=time_ROC_input$time_year, #生存时间(dead和alive的生存时间).
delta=time_ROC_input$event, #生存结局，Censored的样本必须用0表示
marker=time_ROC_input$RiskScore, #预测的变量，这里是风险评分，在没有特殊情况下，值越大，越容易发生危险事件
cause=1, #阳性结局的赋值（必须是1或2），也就是dead的赋值，这里dead是1表示的
weighting = "marginal", #权重计算方法，这是默认方法，选择KM计算删失分布，weighting="aalen" [选用COX]，weighting="aalen" [选用Aalen]
times = c(3,5,10), #计算3、5、10年的ROC曲线
ROC=TRUE,
iid=TRUE #计算AUC
)
time_ROC #查看结果，可以看到，还包括了SE
#绘制ROC曲线啦
summary(time_ROC) #返回12个参数
time_ROC$AUC
#绘制3年的ROC
plot(time_ROC,time=3,col="red",title=FALSE,lwd=2) #将生成一条两倍于 默认宽度的线条
#在此基础上添加5年的ROC
plot(time_ROC,time=5,col="blue",add=TRUE,title=FALSE,lwd=2)
#add 10年的ROC
plot(time_ROC,time=10,col="black",add=TRUE,title=FALSE,lwd=2)
#添加图例
?legend
legend("bottomright", #图例设置在右下角
c(paste0("AUC at 3 years = ", round(time_ROC$AUC[1],3)),
paste0("AUC at 5 years = ", round(time_ROC$AUC[2],3)),
paste0("AUC at 10 years = ", round(time_ROC$AUC[3],3))),
col=c("red","blue","black"),lwd=2,bty="n")
#3. 在timeroc包分析基础上使用ggplot2来画图
library(ggplot2)
time_ROC$TP
summary(time_ROC)
#整理成数据框，gpplot需要数据框
time_ROC.res<-data.frame(TP_3year=time_ROC$TP[,1], #获取3年的ROC的TP
FP_3year=time_ROC$FP[,1], #获取3年的ROC的FP
TP_5year=time_ROC$TP[,2], #获取5年的ROC的TP
FP_5year=time_ROC$FP[,2], #获取5年的ROC的FP
TP_10year=time_ROC$TP[,3], #获取10年的ROC的TP
FP_10year=time_ROC$FP[,3]) #获取10年的ROC的FP
time_ROC$AUC #这里放了3，5，10年的AUC，后续分别取子集time_ROC$AUC[[1]]，time_ROC$AUC[[2]]，time_ROC$AUC[[3]]
TimeROC_plot<-ggplot()+
geom_line(data=time_ROC.res,aes(x=FP_3year,y=TP_3year),size=1,color="red")+
geom_line(data=time_ROC.res,aes(x=FP_5year,y=TP_5year),size=1,color="blue")+
geom_line(data=time_ROC.res,aes(x=FP_10year,y=TP_10year),size=1,color="black")+
geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey",size=1, linetype = 2 #添加虚线
)+
theme_bw()+
annotate("text",x=0.75,y=0.25,size=4.5,
label=paste0("AUC at 3 years = ",round(time_ROC$AUC[[1]],3)),color="red")+
annotate("text",x=0.75,y=0.15,size=4.5,
label=paste0("AUC at 5 years = ",round(time_ROC$AUC[[2]],3)),color="blue")+
annotate("text",x=0.75,y=0.05,size=4.5,
label=paste0("AUC at 10 years = ",round(time_ROC$AUC[[3]],3)),color="black")+
labs(x="False positive rate",y="True positive rate")+
theme(axis.text=element_text(face="bold", size=11, color="black"),#加粗刻度标签
axis.title=element_text(face="bold", size=14, color="black")) #加粗xy轴标签名字
TimeROC_plot
往期回顾
TCGA+biomarker——常见结果展示
TCGA+biomarker——Sample基线表
TCGA+biomarker——单因素Cox回归
TCGA+biomarker——多因素Cox回归
TCGA+biomarker——Cox回归森林图
TCGA+biomarker——Calibration curve
TCGA+biomarker——C-index
TCGA+biomarker——决策曲线分析法（DCA
更多内容可关注公共号“YJY技能修炼”~~~
