预测模型ROC、timeROC曲线，timeAUC顶刊全套路与美图代码

二、ROC(AUC)

受试者工作特征曲线（receiver operating characteristic curve, 简称ROC曲线）是一种评价分类方法/模型的可视化工具，是预测模型评价体系中最核心的指标。ROC曲线下面积AUC值，越高越高。通常发文需要AUC>0.75，高分如PNAS杂志渠道0.99，低分如WJG杂志只有0.56也见刊了。

ROC曲线其实代表了无数个分类器。ROC曲线的横坐标和纵坐标其实是不相关性的，所以不能把ROC曲线当做一个函数曲线来理解，应该把ROC曲线看成无数个点，每个点都代表一个分类器，每个点对应的xy轴代表了这个分类器的性能。ROC曲线就是分类器性能随着阈值(cutoff)的变化而变化的过程。对于ROC曲线，一个重要的特征是它的曲线下面积（AUC），AUC=0.5为随机分类（蓝色random assume线条），识别能力为0，面积越接近于1识别能力越强，面积等于1为完全识别。如图所示，黑色model1曲线AUC>红色model2曲线 → 区分能力model1>model2。左上角越近,证明分类器效果越好。实际上，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）也规范了预测模型的报告过程及质量评价。而依赖time的ROC更复杂，主要依赖是survivalROC 和timeROC包进行时间依赖ROC绘制。

三、常规ROC款式

我阅读了pubmed几百篇SCI，截取收集了全部AUC、ROC套路，部分总结如下，完全版在文末。

四、timeROC timeAUC款式

与time相关ROC曲线类似多组ROC的绘制。更多时候文章中用的是 timeAUC。我阅读了pubmed几百篇SCI，截取收集了全部timeAUC、timeROC套路，部分总结如下，完全版在文末。

五、cutoff诊断试验指标体系

用于确定最佳的指标阈值（也称为cutoff值），以获得最佳的分类效果。通过ROC曲线，可以找到最佳的阈值，使得在敏感性和特异性之间达到最佳的平衡点，从而实现最佳的分类效果。基于ROC寻找cutoff是计算指标体系是个大工程，目前的R包存在各种局限。比如只有点估计值，没有置信区间，只有常规指标没有不平衡指标，基于经验法没有更优秀的bootstrap估计等等。A款智荟代码，一键bootstrap最优算法计算各种指标点估计与95%CI，自动存为table方便写作。

六、智荟ROC-AUC款式

如何正确快速美观的绘制ROC是个技术活。尤其是预测模型，分为训练验证，同数据不同数据，展示效果不一样的时候，还有放大缩小图等等。以及根据cutoff计算的后续IDI NRI等等，都是坑。我开发了许多一键输出R代码，部分总结如下，针对上面pubmed中SCI套路，一一对应，自动绘制套路曲线与p值。大千世界，一念之间，掌握随意。

七、发散思维

1、ROC分常规与timeROC，同理AUC和timeAUC，前者的AUC=Cindex，后者需要注意！timeAUC≠ cindex 而经常= harrel Cindex，timeAUC面积值计算还有如下几种。

2.ROC常用做平衡数据预测模型，存在不平衡问题时候，需要PR曲线互为补充。原创图形如下。

八、智荟R代码

福利来咯，我总结的ROC代码如下，当然非2025最新代码，最棒在A款。

#准备数据
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"))

#----混淆矩阵→ 正态分布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)
#引入 epiR3升级宏
rval <- epi.tests.new(data, conf.level = 0.95)
# 3位数诊断试验数据
print.epi.tests(rval)
#更全指标
summary.epi.tests(rval)

#------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
 
 
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"))