顶刊JCO的ROC曲线绘制秘诀，让你的临床研究再上一个level！

碎碎念专栏帷幕：多指标联合ROC曲线

大家好，我是风。这里是风风的碎碎念专栏。我们知道，在临床研究中，单一一个指标往往较难准确对患者的结局进行判断，为了更加精确判断患者的预后，我们往往会利用多个指标来构建模型，得到一个新的指标，从而期望得到更加准确的价值来判断患者预后。我的导师曾经让我写一个Meta分析，给的方向是根据8个临床指标来联合预测特定肿瘤患者的预后情况，一开始我的做法就是构建临床预后模型然后建立一个新的指标，这种做法大家比较熟悉的就是风险模型的RiskScore值了对吧？这时候我导又问了一个新的问题，如果这个新指标预测能力不佳，有没可能再跟目前公认的金标准指标联合，从而得到更准确的预测效果呢？

听上去是非常典型的临床思维啊对吧？但是该怎么做呢？如果是你，会怎么做呢？后来我记得是参考了一篇Journal of Clinical Oncology(IF 32.596)的文章，做法是联合新的marker和金标准绘制双因素的ROC曲线。

具体要如何操作呢？我们一起来看一下：

多指标ROC曲线

ROC曲线大家都非常熟悉了，这里我不再多做介绍，我们直接绘制一个多指标的ROC曲线，这次我们用一个新的R包——timeROC包，首先是当然是查看数据啦，我们使用一份包含了新marker和免疫检查位点表达量的数据来预测患者的预后：

# install.packages("timeROC")Sys.setenv(LANGUAGE = "en") #显示英文报错信息options(stringsAsFactors = FALSE) #禁止chr转成factorlibrary(timeROC)library(survival)library(tidyverse)# 首先读取数据（起名太麻烦我们就用标志性的feng吧）feng_36 <- read_tsv("MulROC.txt")head(feng_36)## # A tibble: 6 x 10## ID status time_year Marker1 LAG3 IDO1 PDCD1 CD44 CTLA4 IDO2## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>## 1 SA595602 0 1.58 -0.0796 16.3 12.4 20.0 11.8 11.8 11.8## 2 SA595606 0 0.583 -0.172 15.3 11.8 20.8 12.1 11.8 11.8## 3 SA595608 0 0.5 -1.02 15.8 11.8 20.8 12.2 11.8 11.8## 4 SA595600 0 2 -0.239 14.3 12.1 20.6 11.8 11.8 12.1## 5 SA595590 0 0.833 0.594 16.3 18.5 15.0 12.0 11.8 12.0## 6 SA595598 0 0.583 -0.273 16.7 15.0 20.0 11.8 11.8 11.8rocCol=c("#E41A1C", "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F")# ID 患者编号# status 生存状态，0为生存，1为死亡# time_year 生存时间，以年为单位# Marker1 构建的新marker，可以是RiskScore，也可以是其他指标# LAG3 IDO1 PDCD1等等都为免疫检查位点分子表达量

# 接着我们绘制ROC曲线，相较于survivalROC，timeROC不需要提前拟合生存分析，而且可以直接指定生存分析时间，从而一次绘制多条ROC曲线，比如我们要在一张图上展示1、2、3年的多指标ROC曲线，这里使用Marker1、LAG3、IDO1进行绘制：

feng_multi_ROC <- timeROC(feng_36$time_year, # 生存时间 feng_36$status, # 生存状态 feng_36$Marker1, # 绘制的指标 other_markers = as.matrix(feng_36[,c("LAG3","IDO1")]), # 矫正的协变量 cause = 1, # 结局事件对应的值，这里为死亡对应的“1” weighting = "cox", # 计算权重，可以是cox（采用Cox法），也可以设置为marginal（采用KM法） times = c(1,2,3), # 绘制ROC曲线的时间节点，这里会1，2，3年ROC曲线 ROC = TRUE)feng_multi_ROC # 查看ROC结果## Time-dependent-Roc curve estimated using IPCW (n=196, without competing risks). ## Cases Survivors Censored AUC (%)## t=1 7 175 14 70.37## t=2 21 106 69 69.11## t=3 26 50 120 76.97## ## Method used for estimating IPCW:cox ## ## Total computation time : 0.83 secs.# 得到了结果之后，接下来就是绘制ROC曲线进行可视化了，使用plot画图plot(feng_multi_ROC, time = 1, title = FALSE, col = rocCol[1], lwd = 2)plot(feng_multi_ROC, time = 2, title = FALSE, col = rocCol[2], lwd = 2, add = TRUE)plot(feng_multi_ROC, time = 3, title = FALSE, col = rocCol[3], lwd = 2, add = TRUE)# 加上图例legend("bottomright", c(paste0("1 Year AUC = ", signif(feng_multi_ROC$AUC[1],2)), paste0("2 Year AUC = ", signif(feng_multi_ROC$AUC[2],2)), paste0("3 Year AUC = ", signif(feng_multi_ROC$AUC[3],2))), col = rocCol, bty = 'n', lwd = 2, cex = 0.8)

多ROC曲线对比

在绘制完ROC曲线之后，我们往往想将新的结果与其他单一指标进行对比，同样可以使用timeROC进行绘制，通过构建新指标ROC曲线，最后再用plotAUCcurve可视化多个时间节点的AUC值，具体如下：

feng_multi_ROC <- timeROC(feng_36$time_year, feng_36$status, feng_36$Marker1, other_markers = as.matrix(feng_36[,c("LAG3","IDO1","PDCD1")]), cause = 1, weighting = "marginal", times = c(0,0.5,1,1.5,2,2.5,3,3.5,4), ROC = TRUE)feng_multi_ROC # 查看ROC结果## Time-dependent-Roc curve estimated using IPCW (n=196, without competing risks). ## Cases Survivors Censored AUC (%)## t=0 0 196 0 NA## t=1 7 175 14 70.31## t=2 21 106 69 68.61## t=3 26 50 120 77.99## t=4 28 16 152 81.53## ## Method used for estimating IPCW:marginal ## ## Total computation time : 0.01 secs.# 构建其他指标的ROC分析feng_1_ROC <- timeROC(feng_36$time_year, feng_36$status, feng_36$LAG3, cause=1, weighting="marginal", times=c(0,0.5,1,1.5,2,2.5,3,3.5,4), iid=TRUE)feng_1_ROC## Time-dependent-Roc curve estimated using IPCW (n=196, without competing risks). ## Cases Survivors Censored AUC (%) se## t=0 0 196 0 NA NA## t=1 7 175 14 79.72 6.15## t=2 21 106 69 63.23 6.83## t=3 26 50 120 69.40 6.78## t=4 28 16 152 64.62 11.42## ## Method used for estimating IPCW:marginal ## ## Total computation time : 0.18 secs.feng_2_ROC <- timeROC(feng_36$time_year, feng_36$status, feng_36$IDO1, cause=1, weighting="marginal", times=c(0,0.5,1,1.5,2,2.5,3,3.5,4), iid=TRUE)feng_2_ROC## Time-dependent-Roc curve estimated using IPCW (n=196, without competing risks). ## Cases Survivors Censored AUC (%) se## t=0 0 196 0 NA NA## t=1 7 175 14 51.91 10.35## t=2 21 106 69 56.06 6.68## t=3 26 50 120 61.69 6.24## t=4 28 16 152 59.69 7.96## ## Method used for estimating IPCW:marginal ## ## Total computation time : 0.23 secs.feng_3_ROC <- timeROC(feng_36$time_year, feng_36$status, feng_36$PDCD1, cause=1, weighting="marginal", times=c(0,0.5,1,1.5,2,2.5,3,3.5,4), iid=TRUE)feng_3_ROC## Time-dependent-Roc curve estimated using IPCW (n=196, without competing risks). ## Cases Survivors Censored AUC (%) se## t=0 0 196 0 NA NA## t=1 7 175 14 28.63 10.69## t=2 21 106 69 44.35 7.75## t=3 26 50 120 30.57 6.94## t=4 28 16 152 40.59 9.28## ## Method used for estimating IPCW:marginal ## ## Total computation time : 0.08 secs.feng_4_ROC <- timeROC(feng_36$time_year, feng_36$status, feng_36$CTLA4, cause=1, weighting="marginal", times=c(0,0.5,1,1.5,2,2.5,3,3.5,4), iid=TRUE)feng_4_ROC## Time-dependent-Roc curve estimated using IPCW (n=196, without competing risks). ## Cases Survivors Censored AUC (%) se## t=0 0 196 0 NA NA## t=1 7 175 14 46.29 0.99## t=2 21 106 69 46.70 1.21## t=3 26 50 120 48.22 2.93## t=4 28 16 152 45.52 4.50## ## Method used for estimating IPCW:marginal ## ## Total computation time : 0.07 secs.feng_5_ROC <- timeROC(feng_36$time_year, feng_36$status, feng_36$CD44, cause=1, weighting="marginal", times=c(0,0.5,1,1.5,2,2.5,3,3.5,4), iid=TRUE)feng_5_ROC## Time-dependent-Roc curve estimated using IPCW (n=196, without competing risks). ## Cases Survivors Censored AUC (%) se## t=0 0 196 0 NA NA## t=1 7 175 14 38.00 4.95## t=2 21 106 69 45.38 5.55## t=3 26 50 120 45.52 5.08## t=4 28 16 152 50.61 9.21## ## Method used for estimating IPCW:marginal ## ## Total computation time : 0.07 secs.feng_6_ROC <- timeROC(feng_36$time_year, feng_36$status, feng_36$IDO2, cause=1, weighting="marginal", times=c(0,0.5,1,1.5,2,2.5,3,3.5,4), iid=TRUE)feng_6_ROC## Time-dependent-Roc curve estimated using IPCW (n=196, without competing risks). ## Cases Survivors Censored AUC (%) se## t=0 0 196 0 NA NA## t=1 7 175 14 49.40 7.26## t=2 21 106 69 50.57 4.92## t=3 26 50 120 50.99 4.03## t=4 28 16 152 49.10 4.94## ## Method used for estimating IPCW:marginal ## ## Total computation time : 0.06 secs.最后将多个ROC曲线的AUC结果具现在一张图上，这次我们用plotAUCcurve可视化：plotAUCcurve(feng_multi_ROC, conf.int=FALSE, col=rocCol[1])plotAUCcurve(feng_1_ROC, conf.int=FALSE, col=rocCol[2], add=TRUE)plotAUCcurve(feng_2_ROC, conf.int=FALSE, col=rocCol[3], add=TRUE)plotAUCcurve(feng_3_ROC, conf.int=FALSE, col=rocCol[4], add=TRUE)plotAUCcurve(feng_4_ROC, conf.int=FALSE, col=rocCol[5], add=TRUE)plotAUCcurve(feng_5_ROC, conf.int=FALSE, col=rocCol[6], add=TRUE)plotAUCcurve(feng_6_ROC, conf.int=FALSE, col=rocCol[7], add=TRUE)# 图例设置legend("bottomright", c("Marker1","LAG3","IDO1","PDCD1","CD44","CTLA4","IDO2"), col = rocCol, bty = 'n', lty = 1, lwd = 2, cex = 0.8)

观察结果，可以发现联合多个指标在早期比单独使用LAG3预测效果差，但是在晚期则预测效果明显由于其他单个指标的预测效果，说明联合Markers等多个指标对患者远期生存的预测效果更好，这样我们多指标ROC曲线以及ROC曲线对比也就到此完成啦，将上面的图形放进文章，相信会让你的模型更有说服力，让你的临床研究也更容易打动审稿人。

好了，这样我们今天内容到此结束，这也是春节之前风风专栏最后一次推文啦。

后台回复“feng36”获得本次代码和数据，下一期我们将开启全新的篇章，我们下次再见吧！*^_^*