TCGA+biomarker——ROC曲線

通常情況下，通過以下幾種指標來對模型進行評價。

1）區分度：採用指標C-index和ROC曲線來評價區分度，一般文章都是二選一。

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曲線圖解

解釋：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=樣本總數。

(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曲線是如何繪製的。

如上圖所示，共計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, 說明該模型分類器的準確性處於中等水平。

目的2——比較不同模型的優劣

下圖中有3個模型的ROC曲線，根據AUC面積大小，可直觀展示出3個模型之間的優劣，A模型>B模型>C模型。

但如下圖，上圖中兩條ROC曲線相交於一點，AUC值幾乎一樣：當需要高Sensitivity時，模型A比B好；當需要高Speciticity時，模型B比A好。

目的3——確定最佳臨界點

所謂找到最優臨界點，就是保證真陽性率（敏感高）高的同時假陽性率（特異性）要儘量的小，建立max（TPR+（1-FPR））的模型。例如：交叉點坐標為（95.8%，14.6%），相應的臨界值即可出來，說明敏感度是95.8%，特異度是84.6%。

如何繪製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 rateplot(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的線#計算最佳cutoffcutoff<-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.inputtime_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年的ROCplot(time_ROC,time=3,col="red",title=FALSE,lwd=2) #將生成一條兩倍於 默認寬度的線條#在此基礎上添加5年的ROCplot(time_ROC,time=5,col="blue",add=TRUE,title=FALSE,lwd=2)#add 10年的ROCplot(time_ROC,time=10,col="black",add=TRUE,title=FALSE,lwd=2)#添加圖例?legendlegend("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$TPsummary(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的FPtime_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基線表
長按識別二維碼，關注"YJY技能修煉"