Cox构建预测模型（6）：如何用R语言绘制ROC曲线？（附全套代码）

临床预测模型作为临床研究的“高阶玩法”，不仅仅是改变临床实践的重要途径，更是发表高分SCI文章的热门选择。

但不论零基础的小白，还是已经了解过临床预测模型的同学，刚开始都会一头雾水。简单概括，Cox回归预测模型的基础统计策略大致可以概括为“一表四图”，即均衡性表、列线图、校准图、ROC图、DCA图。

之前的文章中我们已经为大家介绍过均衡性表、列线图、校准图的绘制方法，今天继续复现seer公共数据库文章，重点介绍预测模型验证三件套之ROC曲线的绘制，除了R语言复现，同样用风暴统计进行操作，对新手小白或希望更加便捷完成统计分析工作的读者十分友好！

主要内容包括： 一、文献解读 二、利用R语言复现 三、利用在线网站复现   四、小结

一、文献解读

案例文献是沈阳医学院公共卫生学院学者基于SEER数据库的一项回顾性研究，旨在建立一个列线图来预测老年恶性骨肿瘤（MBT）患者的总生存期（OS）。

1. 摘要

背景：恶性骨肿瘤（MBT）是老年患者死亡的原因之一。我们研究的目的是建立一个列线图来预测老年MBT患者的总生存期（OS）。

方法：从SEER数据库下载了2004年至2018年所有老年MBT患者的临床病理数据。他们被随机分配到训练集（70%）和验证集（30%）。采用单因素和多因素Cox回归分析确定老年MBT患者的独立危险因素。基于这些危险因素构建列线图，以预测老年MBT患者的1年，3年和5年OS。然后，利用一致性指数（C指数）、校准曲线和受试者工作曲线下面积（AUC）来评价预测模型的准确性和判别力。决策曲线分析（DCA）用于评估列线图的临床潜在应用价值。根据列线图上的分数，将患者分为高风险组和低风险组。Kaplan-Meier（K-M）曲线用于测试两名患者之间的生存差异。

结果：从SEER数据库下载了2004年至2018年所有老年MBT患者的临床病理数据。他们被随机分配到训练集（70%）和验证集（30%）。采用单因素和多因素Cox回归分析确定老年MBT患者的独立危险因素。基于这些危险因素构建列线图，以预测老年MBT患者的1年，3年和5年OS。然后，利用一致性指数（C指数）、校准曲线和受试者工作曲线下面积（AUC）来评价预测模型的准确性和判别力。决策曲线分析（DCA）用于评估列线图的临床潜在应用价值。根据列线图上的分数，将患者分为高风险组和低风险组。Kaplan-Meier（K-M）曲线用于测试两名患者之间的生存差异。

结论：我们建立了一个新的列线图来预测老年MBT患者的1年，3年，5年的OS。该预测模型可以帮助医生和患者制定治疗计划和后续策略。、

2. 数据介绍

文献共纳入1641名2004-2018诊断的老年MBT患者的临床病理数据。暴露因素包括年龄、种族、性别、诊断年份、组织学类型、分级、分期、原发位置、TNM分期、肿瘤大小、是否化疗、是否放疗及手术方式。

3.研究结果

这篇文献构建Cox回归预测模型的统计思路十分清晰。首先按照7：3将数据集进行拆分获得训练集与验证集，然后做均衡性检验，比较训练集和验证集的差异性，再做单因素和多因素Cox回归，筛选变量构建列线图预测模型，最后对模型进行评价。

文献中通过ROC曲线对预测模型进行区分度验证，A图是训练集1年、3年、5年患者的总生存期的ROC曲线；B图是验证集1年、3年、5年患者的总生存期的ROC曲线。ROC曲线下面积又叫AUC面积，这个值越大说明预测模型的而判别区分能力越好。

二、利用R语言复现

下面我们利用从SEERStat中提取的相关数据进行复现。根据文献中的纳入排除标准，提取涉及的相关数据，最终共纳入1，574名患者（原文献1，641）。介于SEERStat数据库会有更新，因此提取的样本量与原文会有所出入，这里请大家多关注统计方法的运用！本次用到的是R版本是4.3.1。

1. 安装加载R包

if(!require("timeROC")) install.packages("timeROC")
library("timeROC")

2.完成前期准备工作

在绘制ROC曲线之前，必须要完成数据的导入、拆分数据集的工作、预测因子的筛选，校准图的绘制是分别在训练集与验证集独立进行的。

以上这些内容可参见：

Cox构建预测模型（2）：如何用R语言拆分数据集并作均衡性检验

Cox构建预测模型（3）：如何用R语言做cox先单后多回归分析！（附全套代码）

Cox构建预测模型（4）:如何用R语言绘制cox预测模型列线图？（附全套代码）

Cox构建预测模型（5）：如何用R语言绘制校准曲线（附全套代码）

3. 绘制ROC曲线

这里根据文献，我们需要绘制1年、3年、5年的ROC曲线，为减少篇幅，只为大家展示训练集ROC曲线的绘制，验证集同理。从下面的代码可以看出通过R语言绘制ROC曲线还是较为复杂的，对于代码的含义，小编也进行了标注，大家可以辅助进行理解。

train$lp <- predict(train_model, newdata = train,type ="lp")
time_roc_train <- timeROC(
  T = train$months,  #结局时间
  delta = train$status,  #生存结局 
  marker = train$lp,  #预测变量
  cause = 1,  #阳性结局指标值
  weighting="marginal",  #权重计算方法，默认marginal，采用KM计算删失分布
  times = c(12, 12*3, 12*5),
  ROC = TRUE,  #保存sensitivities 和 specificties值
  iid = TRUE  #iid = TRUE 才会保存置信区间，但是样本量大了后，耗时耗资源
)


#计算train-1年的AUC
train1y <- paste0("train 1 year AUC (95%CI) = ", 
                  sprintf("%.3f",time_roc_train$AUC[1]) ," (",
                  sprintf("%.3f",confint(time_roc_train, level = 0.95)$CI_AUC[1,1]/100)," - ",
                  sprintf("%.3f",confint(time_roc_train, level = 0.95)$CI_AUC[1,2]/100),")")


#计算train-3年的AUC
train3y <- paste0("train 3 year AUC (95%CI) = ", 
                  sprintf("%.3f",time_roc_train$AUC[2]) ," (",
                  sprintf("%.3f",confint(time_roc_train, level = 0.95)$CI_AUC[2,1]/100)," - ",
                  sprintf("%.3f",confint(time_roc_train, level = 0.95)$CI_AUC[2,2]/100),")")


#计算train-5年的AUC
train5y <- paste0("train 5 year AUC (95%CI) = ", 
                  sprintf("%.3f",time_roc_train$AUC[3]) ," (",
                  sprintf("%.3f",confint(time_roc_train, level = 0.95)$CI_AUC[3,1]/100)," - ",
                  sprintf("%.3f",confint(time_roc_train, level = 0.95)$CI_AUC[3,2]/100),")")


#train AUC 多时间点合一
plot(title="",time_roc_train,col="DodgerBlue",time=12,lty=1,lwd = 2,family = "serif") #绘制ROC曲线       
plot(time_roc_train,time=12*3,lty=1,lwd = 2,add=TRUE,col="LightSeaGreen",family = "serif")#add=TRUE指在前一条曲线基础上新增
plot(time_roc_train,time=12*5,lty=1,lwd = 2,add=TRUE,col="DarkOrange",family = "serif" )
legend("bottomright",c(train1y,train3y,train5y),col=c("DodgerBlue","LightSeaGreen","DarkOrange"),lty=1,lwd=2)

代码解读：计算过AUC值后，将3条ROC曲线绘制在同一张表格中，关键指标为“add=TRUE”，否则仅绘制单条曲线，legend表示注释表格的位置，"bottomright"指右侧底部。

R语言复现得到的ROC曲线如下：

三、利用在线网站复现

如果没有代码基础，或者希望通过更便捷的方式完成统计分析，推荐使用这个智能在线统计分析平台——风暴统计。一键成三线表，便捷又快速。

网址：www.medsta.cn（电脑端浏览器打开） 或medsta.cn（medical statistics缩写）

1.进入网站分析模块

电脑端打开风暴统计平台——“风暴智能统计”模块，点击“临床预测模型（最新）”，进入“cox预测模型”页面。

2.完成前期的准备工作

利用风暴统计平台进行ROC曲线的绘制之前，也需要完成数据导入、拆分数据集的工作哦！

3.绘制ROC曲线

完成前期准备工作之后，进入“Cox预测模型”模块，根据提示选入回归因变量与回归自变量。最最重要的是多因素回归设置，也就是你所期望的自变量筛选方式，可以是先单后多，也可以是逐步回归法，阈值也可以进行调节，十分灵活。

选中之后，下方就直接给出ROC曲线结果，包括训练集与验证集3个不同时间点的6张ROC曲线图，在左侧可以直接调整图形的线条粗细，AUC95%可信区间注释的位置。小白式操作也可以绘制出精美的图形！

四、小结

风暴统计网站给出的图线与AUC面积及95%置信区间和R语言结果完全一致，只是网站目前不能够把同一模型的多条曲线绘制在一张图中。但是与R语言几十行代码相比，风暴统计足够简洁、快速，结果可靠，在实际生活中，还十分值得一试的！