ROC曲线，套路全解析（一）
本章思维导图：
前言
ROC(receiver operating characteristic curve) 受试者工作曲线，又称为感受性曲线（sensitivity curve）是医学诊断试验、预测模型性能区分度评价的最核心指标(见表1)1。ROC曲线其实代表了无数个分类器。ROC曲线的横坐标和纵坐标其实是不相关性的，所以不能把ROC曲线当做一个函数曲线来理解，应该把ROC曲线看成无数个点，每个点都代表一个分类器，每个点对应的xy轴代表了这个分类器的性能。ROC曲线就是分类器性能随着阈值(cutoff)的变化而变化的过程。对于ROC曲线，一个重要的特征是它的曲线下面积（AUC），AUC=0.5为随机分类（蓝色random assume线条），识别能力为0，面积越接近于1识别能力越强，面积等于1为完全识别。如图所示，黑色model1曲线AUC>红色model2曲线 → 区分能力model1>model2（见图1）。可见曲线距离左上角越近,证明分类器效果越好。实际上，AUC>0.9是及其优秀，AUC 0.8-0.9 优秀，0.7-0.8 可接受，0.5-0.7 较差模型性能。要计算ROC，对于设计类型为诊断试验的研究，金标准分组Y一般来说是二分类，试验变量可以为连续型、分类、有序分类。对于设计类型为预测模型的研究，结局Y可以依赖随访时间变量time成为预后（prognosis）预测，也可以不依赖随访time成为诊断（diagnosis）预测，后者与诊断试验统计分析方法一致2。TRIPOD申明（Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis）也规范了预测模型的报告过程及质量评价（见图2）3。而依赖time的ROC更复杂4，主要依赖是survivalROC 和timeROC包进行时间依赖ROC绘制。本章节阐述非时间依赖型ROC（Logistic回归），下一章节阐述时间依赖型（Kaplan-Meier Cox回归）。
表1 预测模型模型性能评价区分度、校准度、重分类度
图1 ROC曲线是无数个cutoff点集合而成
图2 TRIPOD申明与ROC
ROC切点
ROC自动最佳切点，有多种计算方法，最主要是依据Youden指数。Youden指数依据ROC指标体系中的敏感度、特异度来获得，最佳切点算法还可以“closest.topleft”，具体产生过程和指标概念见表2、表3。其中灵敏度、特异度、准确率、阳性预测值、阴性预测值是除了AUC外最常用的5个ROC评价指标。ROC自动切点除了pROC外有OptimalCutpoints、Epi包，但是笔者测试后面2个有bug，结果有时候会不够精准与pROC及其他软件结果不一致，最好不用。当然，也可以MedCalc、SPSS模块根据ROC结果坐标表格手动计算youden指数。SAS软件可自行编程。
表2 诊断四格表两种形态
表3 ROC诊断指标体系
曲线比较
ROC曲线比较是采用AUC95%CI是否有重叠来简单判断异。诊断试验中同一对象，2种方法配对计算P 值，即Delong-method。预测模型中，训练集同一结局指标，不同预测因素也是配对Delong-method，而如果要比较训练集和验证集ROC性能差异，此时二者为独立数据集则应用Z检验、bootstrap、venkatrama法。单独ROC曲线与45°的chance line (或者 random assume line)比较有差异，只能证明此ROC有一定区分度。ROC曲线有原始曲线和平滑曲线两种，计算的AUC和图形状态均不一样。ROC 95%CI估计有参数、半参数、非参数法。常规的ROC曲线用的是非参数法，如empirical method，有些用参数法，如binormal method5。绘制ROC曲线有ROCR、Rplot、pROC等，首推pROC包6。pROC包也同时支持多分类ROC，另外python的micro-average和macro-average 也支持多分类ROC。也有文章提到校正混杂后ROC曲线，用ROCt包7，目前已经下架。笔者查看原始代码，感觉就是建模多了个混杂比较校正不校正两条曲线的差异。
论文ROC
初步检索pubmed，ROC曲线展示大致有如下几种。经典款、带error bar、带95%CI、平滑ROC。pROC包示例代码见，https://web.expasy.org/pROC/screenshots.html。
图3 主流ROC曲线形态
发散思维
与ROC曲线近似的，还有PR曲线。极度偏态的数据集上，Precision-Recall(PR) Curve可表现更全面8。PR其实就是x轴TPR或recall或Sen，也叫查全率，y轴PPV或precision，也就是查准率。ROC的X轴是FPR,Y轴是TPR。FPR体现的是有多少负类被错抓成了正类，TPR体现的是有多少正类被正确的分类为正类, Precison衡量的是模型判别为正类的样本中，有多少确实就是正类。一条曲线在ROC曲线中压过另一条曲线，那么他在PR曲线中也会相同的全面优于另一条曲线（如图4）。
图4 ROC（左）与PR曲线（右）
ROC代码
pacman::p_load(caret,rms,forswang1gn,data.table,SparseM,TH.data,readxl,Hmisc,quickReg,rms,lattice,survival,Formula,ggplot2,pscl)
setwd("D:/写书/一键ROC")
load("D:/写书/一键ROC/rhc.rdata")
#5735人，按照Y划分7:3，效果比简单随机好
set.seed(1234)
library(caret)
train.id=createDataPartition(rhc$death,p=0.7)$Resample1
#划分train,test数据集
train=rhc[train.id,]
test=rhc[-train.id,]
#y不能设为factor，其他随意
#快速获取变量名向量c("")
dput(names(rhc))
#surv2md1 存活概率反向
conVars <- c("age","aps1","crea1","surv2md1")
catVars <- c("sex","swang1")
Vars<-c(conVars,catVars)
train$y<-train$death_01
#分类变量转为→factor因子，后面dummy哑变量用
rhc[catVars] <- lapply(rhc[catVars], factor)
#--ROC计算包有ROCR,mlr3,plotROC、pROC,本章选最优解pROC
#------MODEL1 p1 多因素logistic获取p1-----------------
#logistic生成LR似然值用来后续nomogram;rm或glm函数
#此处建模可用ML genmod Cox等等，预测模型只是为得到预测概率p
F1<-glm(y~aps1+sex+swang1,train, family = binomial(link = "logit"))
train$p1 <- predict(F1,newdata =train, type = "response")
test$p1 <- predict(F1,newdata =test, type = "response")
#------MODEL2 F2= model1+新指标-----------------
F2<-glm(y~aps1+sex+swang1+aps1+crea1+age,train, family = binomial(link = "logit"))
train$p2 <- predict(F2,newdata =train, type = "response")
test$p2 <- predict(F2,newdata =test, type = "response")
#------ROC 曲线比较
library(pROC)
#smooth = T采用平滑技术,F 就和其他软件一样,但是平滑后计算CI要用smooth.roc ，慎用.
roc1_train <- roc(death~p1, smooth = F, percent = TRUE, data = train, auc= TRUE, plot = TRUE)
roc2_train <- roc(death~p2, smooth = F, percent = TRUE, data = train, auc= TRUE, plot = TRUE)
#roc，lines还可多次添加曲线
plot(roc1_train)
lines.roc(roc2_train, col= "red")
legend(x = 40, y = 40,fill = c("black", "red"),legend = c("model1", "model2"), cex = 1)
#-------------------------计算PPV NPV SEN SPE等指标，不同R包测试youden
#AUC95%ci，empirical法；method = "bin" 或者"non"结果近似；
library(ROCit)
rocit_emp <- rocit(score = train$p1,class = train$y, method = "emp")
ciAUC(rocit_emp)
#------ROC 原始曲线比较
library(pROC)
#当p1换成单指标时，即为单指标ROC，p1可能logistic或cox或ML生成
roc1_train <- roc(death~p1, smooth = F, percent = TRUE, data = train, auc= TRUE, plot = TRUE)
roc2_train <- roc(death~p2, smooth = F, percent = TRUE, data = train, auc= TRUE, plot = TRUE)
#roc，lines还可多次添加曲线;R中默认会画出增加4%的坐标轴
plot(roc1_train,xlim=c(100,0),ylim=c(0,100),xaxs="i",yaxs="i")
lines.roc(roc2_train, col= "red")
legend(x = 40, y = 40,fill = c("black", "red"),legend = c("model1", "model2"), cex = 1)
#方法1，切点0.651，AUC 0.621；thres="best"默认youden，还可以closest.topleft
auc(roc1_train)
coords(roc1_train, "best", ret=c("threshold", "specificity", "sensitivity", "accuracy","ppv", "npv"), transpose = FALSE)
#ci.auc用R版本3.6.1 来建造的,R4.1出错
#默认delong-method 计算95%CI
#ci(roc1_train)
ci.auc(roc1_train,conf.level=0.95, method=c("bootstrap"),boot.n=1000)
#不能设置seed()
#ppv npv bootstrap 95%CI 由于多次计算bootstrap可能会超过1，bugs;注意bootstrap仅得到2.5% 50%,97.5%界值
ci.coords(roc1_train, x="best", input = "threshold", ret=c("sensitivity", "specificity", "ppv", "npv","accuracy"), conf.level=0.95, boot.n=100,boot.stratified=TRUE)
#方法2，切点0.651，有时候不一致与方法1
library(OptimalCutpoints)
optimal.cutpoint.Youden <- optimal.cutpoints(X="p1", status ="death_01", tag.healthy = 0, methods = "Youden", data =train)
summary(optimal.cutpoint.Youden)
plot(optimal.cutpoint.Youden)
#方法3，切点0.651,有时候不一致与方法1
library(Epi)
ROC(form= y~aps1+sex+swang1, data=train, plot="ROC", PV=TRUE, MX=TRUE)
#----ROC平滑处理得不到“threshold”；平滑方法有binormal", "density", "fitdistr"等
roc1_train_smooth <- roc(death~p1, smooth = T,smooth.method="binormal", percent = TRUE, data = train, auc= TRUE, plot = TRUE)
roc2_train_smooth <- roc(death~p2, smooth = T,smooth.method="binormal",percent = TRUE, data = train, auc= TRUE, plot = TRUE)
plot(roc1_train_smooth, col= "brown")
lines.roc(roc2_train_smooth, col= "orange")
legend(x = 40, y = 40,fill = c("brown", "orange"),legend = c("model1", "model2"), cex = 1)
#切点不详，AUC 0.621
auc(roc1_train_smooth)
ci.auc(roc1_train_smooth,conf.level=0.95, method=c("bootstrap"),boot.n=1000)
#ci.coords(roc1_train_smooth, x="best", input = "threshold", ret=c("sensitivity", "specificity", "ppv", "npv","accuracy"))
图5 原始ROC曲线(左)与平滑ROC曲线(右)
#----------------混淆矩阵→ 正态分布95%CI 诊断指标体系 library(ROCR)后手动计算，
#备份 P>0.5 empirical roc 分2类，但是一般不用。
#方法1，ROCR+epiR包混搭计算
library(ROCR)
integer<-as.integer(train$p1>=0.651)
logit.perf <- table(integer,train$y)
logit.perf
#integer 0 1
#0 876 1130
#1 534 1476
#epiR包，需要手动输入刚刚的confusion matrix
# 需要确定顺序！+ - + - , 四格表诊断试验输出95%CI，
library(epiR)
data <- as.table(matrix(c(1476,534,1130,876), nrow = 2, byrow = TRUE))
rval <- epi.tests(data, conf.level = 0.95)
print(rval)
#==================引入 epiR升级宏
source("./epi.tests.new-epiR3位数诊断试验.R")
rval <- epi.tests.new(data, conf.level = 0.95)
# 3位数诊断试验数据
print.epi.tests(rval)
#更全指标
summary.epi.tests(rval)
#方法2，手动计算 F1值等，BER=1/2 * (FPR+FN/(FN+TP))，仅点估计无95%CI
#logit.perf 是01,01排序
performance <- function(table,n=4){
if(!all(dim(table)==c(2,2)))
stop('Must be a 2×2 table')
tn=table[1,1]
fp=table[2,1]
fn=table[1,2]
tp=table[2,2]
sensitivity = tp/(tp+fn)
specificity = tn/(tn+fp)
ppv = tp/(tp+fp)
npv = tn/(tn+fn)
hitrate = (tp+tn)/(tp+tn+fp+fn)
F1_score=2*(ppv*sensitivity)/(ppv+sensitivity)
result <- paste('Sensitivity = ',round(sensitivity,n),
'\nSpecificity = ',round(specificity,n),
'\nPositive Predictive Value =',round(ppv,n),
'\nNeswang1tive Predictive Value =',round(npv,n),
'\nAccuracy = ',round(hitrate,n),'\n',sep='',
'F1_score = ',round(F1_score,n))
cat(result)
}
#缺陷是没有CI
performance(logit.perf)
表4 Empirical ROC性能指标体系
#------ROC 配对曲线比较，验证集内
library(pROC)
roc1_train <- roc(death~p1, smooth = F, percent = TRUE, data = train, auc= TRUE, plot = TRUE)
roc2_train <- roc(death~p2, smooth = F, percent = TRUE, data = train, auc= TRUE, plot = TRUE)
#验证集, 获取ROC曲线差异
roc1_test <- roc(death~p1, smooth = F, percent = TRUE, data = test)
roc2_test <- roc(death~p2, smooth =F, percent = TRUE, data = test)
#========roc1_train + roc1_test 同一数据配对ROC
#配对比较用 delong ，paired = TRUE, method = "delong"
roc.test(roc1_train, roc2_train,paired = TRUE, method = "delong")
#0.621,0.679，p< 0.001
# 输出 pdf roc12_train.pdf
pdf("roc12_train.pdf",family="Times", width=6,height=6)
plot(roc1_train,xaxs = "i", yaxs = "i",col= "black")
lines.roc(roc2_train, col= "red")
legend("bottomright", legend=c("model1train","model2train"), col=c("black","red"),lty=1,lwd=2)
dev.off()
#关闭图形
#========roc1_train + roc1_test 配对对比ROC
# 非配对ROC 用bootstrap 或者 method="venkatrama";#smooth = T或F 与train数据需要保持统一
roc.test(roc1_train, roc1_test, method = "bootstrap",boot.n = 100)
#0.621,0.596，p=0.143
# 输出 pdf roc1_traintest.pdf
pdf("roc1_traintest2.pdf",family="Times", width=6,height=6)
plot(roc1_train,xaxs = "i", yaxs = "i") # 制作第一条ROC曲线,原来起
lines.roc(roc1_test, col= "red")
legend("bottomright", legend=c("train","validation"), col=c("black","red"),lty=1,lwd=2)
dev.off()
#关闭图形
#多组曲线ggplot2格式
g2 <- ggroc(list(m1=roc1_train, m2=roc2_train, m3=roc1_test, m4=roc2_test))
g2
图6 基于同一数据集配对ROC与不同数据集独立ROC曲线的比较
图7 ggplot2样式ROC
library(pROC)
rocobj <- plot.roc(train$death, train$p1,main="Confidence intervals", percent=TRUE, ci=TRUE, print.auc=TRUE)
#bootstrap绘制灵敏度置信区间，每5一次共20组CI连成蓝色sen区间
ciobj <- ci.se(rocobj, specificities=seq(0, 100, 5))
plot(ciobj, type="shape", col="#1c61b6AA")
#youden max切点处加标记
plot(ci(rocobj, of="thresholds", thresholds="best"))
图8 带灵敏度置信区间的ROC
小结
掌握pROC包参数，绘制AUC值95%CI（c-index）与模型性能评价5件套指标（sen、spe、acc、ppv、npv）。掌握多种ROC曲线形态，及曲线差异性检验。
原创不易，转载请说明来自本公众号。
参考文献
Verdonschot JAJ, Hazebroek MR, Wang P, Sanders-van Wijk S, Merken JJ, Adriaansen YA, van den Wijngaard A, Krapels IPC, Brunner-La Rocca HP, Brunner HG, Heymans SRB. Clinical Phenotype and Genotype Associations With Improvement in Left Ventricular Function in Dilated Cardiomyopathy. Circ Heart Fail 2018;11(11):e005220.doi:10.1161/CIRCHEARTFAILURE.118.005220 van Smeden M, Reitsma JB, Riley RD, Collins GS, Moons KG. Clinical prediction models: diagnosis versus prognosis. J Clin Epidemiol 2021;132:142-145.doi:10.1016/j.jclinepi.2021.01.009 Collins GS, Reitsma JB, Altman DG, Moons KG. Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD): the TRIPOD statement. BMJ 2015;350:g7594.doi:10.1136/bmj.g7594 Kamarudin AN, Cox T, Kolamunnage-Dona R. Time-dependent ROC curve analysis in medical research: current methods and applications. BMC Med Res Methodol 2017;17(1):53.doi:10.1186/s12874-017-0332-6 Colak E, Mutlu F, Bal C, Oner S, Ozdamar K, Gok B, Cavusoglu Y. Comparison of semiparametric, parametric, and nonparametric ROC analysis for continuous diagnostic tests using a simulation study and acute coronary syndrome data. Comput Math Methods Med 2012;2012:698320.doi:10.1155/2012/698320 Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, Muller M. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics 2011;12:77.doi:10.1186/1471-2105-12-77 https://rdrr.io/cran/ROCt/man/adjusted.ROCt.html. Yan B, Guan D, Wang C, Wang J, He B, Qin J, Boheler KR, Lu A, Zhang G, Zhu H. An integrative method to decode regulatory logics in gene transcription. Nat Commun 2017;8(1):1044.doi:10.1038/s41467-017-01193-0