简介面向医学、生物统计与临床试验中的生存分析场景这份R语言源代码包用于绘制多时间点生存ROC曲线解决传统ROC仅能评估单一时间点的局限适合已有R基础并希望量化生存预测模型性能的研究者使用。压缩包体积仅13KB包含3个文件R脚本负责数据读入、AUC计算和曲线绘制txt文件提供示例数据结构PDF为输出效果参考便于核对运行结果。当前已有351人学习下载。代码以pROC和survival包为基础涵盖数据整理、不同随访时刻ROC曲线叠加、图例标注及AUC对比等关键环节研究者只需准备包含生存时间、事件状态和风险评分的数据即可替换示例文件直接运行快速评估模型在多时间点的区分能力为SCI论文中的预后分析图表提供规范、可复现的实现方案。该工具能显著减少重复编码工作量帮助研究者将精力集中在模型解读与科研结论上。1. 多时间点生存ROC曲线这份R源代码能直接改改就用手头有个临床预测模型的课题模型训练完审稿人一句话让补多时间点生存ROC曲线。常规做法是拿pROC在每个时间点单独画一条ROC画完发现问题很大删失数据根本绕不过去每个时间点能用的样本量都不一样画出来的AUC没法横向比。后来拆到这份R语言源代码资源包里面就三个文件input.txt示例数据、bioR43.multiTimeROC.R主脚本、ROC.pdf效果参照。脚本核心是调用timeROC包一条命令算出多个时间点的AUC再叠画到同一张图上直接输出SCI能用的矢量图。这份资源适合已经有R基础、正在做生存分析模型验证的研究者特别是要应付审稿人“补分析”场景的。拿到手不是让你从头写代码而是把数据套进它的格式改三行参数就能出图。2. 从单时间ROC到多时间生存ROC为什么必须用timeROC而不是pROC2.1 生存ROC和普通ROC的本质区别删失数据怎么处理普通ROC曲线评估的是二分类问题每个样本都有一个确定标签0或1样本之间独立不存在“时间”维度。但生存ROC面对的是随访数据每个样本只有两列关键信息生存时间time和结局状态status。这里有个坑随访结束时还没发生事件的样本状态是删失censored它不代表“永远不发生”只是“在观察窗口内没发生”。如果你硬把删失样本当成0类扔进pROC等于默认这些人在随访期内绝对安全结果会系统性低估模型性能。多时间点生存ROC的解决思路是加权估计。在某个特定时间点t把所有样本划分成两组在t之前发生事件的人是一组在t时刻仍然存活的人是另一组。但对删失样本来说问题来了——有些人在t之前就失访了我们不知道他到底什么时候发生事件。timeROC包的处理方法是IPCW逆删失概率加权给每个样本根据删失分布算一个权重用加权后的伪样本去估计灵敏度sensitivity和特异度specificity从而得到t时刻的ROC曲线和AUC。这套逻辑在timeROC包的文档里写得很清楚核心函数是timeROC()。用pROC做多时间点ROC等于把每个时间点单独拎出来做一次二分类每个时间点的样本集还不一样曲线之间没有一致性。timeROC则是在同一个队列、同一种删失结构下用统一的加权框架估计每个时间点的AUC这才具有横向可比性。这份R源代码选择timeROC方案从方法学上是站得住的。如果你投的期刊有统计审稿人用pROC逐点画图很容易被追问“删失样本怎么处理的”用timeROC自带IPCW权重至少在方法描述上无懈可击。2.2 数据格式是第一道坎先把手里的表改造成生存数据框解压打开input.txt内容是一张典型的纵向随访表格。我按经验整理一下这个格式里最关键的列对应关系如下表。列名类型格式要求说明time数值单位统一为天/月/年生存时间从入组到事件发生或随访截止status数值0或1必须数值型1事件发生死亡/复发/转移0删失或随访结束未事件marker数值连续变量模型预测值比如风险评分、蛋白表达值、signature得分ID字符/数值不与他人重复非必需但建议保留以便核查实际操作里大多数人的原始数据长这样OS_time总生存时间、OS_event是否死亡事件、risk_score模型输出的预测风险。拿到这份资源后第一步不是跑脚本而是把这三列改名为time、status、marker再检查status是不是0/1编码——很多人用Dead/Alive字符串这个后面会报错。# 读取临床数据假设你的文件叫 clin.csv clin - read.csv(clin.csv, stringsAsFactors FALSE) # 检查核心三列的缺失情况和类型 summary(clin$time) table(clin$status) # 确认只有0和1且1代表事件发生 # 如果status是字符串手动转成0/1 clin$status - ifelse(clin$status Dead, 1, 0) # 按资源包的格式约定列名 clin - clin[, c(time, status, marker)] clin - clin[complete.cases(clin), ] # 去掉缺失行 str(clin)逻辑说明这段代码本质是把你的临床原始表标准化成input.txt那种三列结构。summary、table、str三个检查是为了避免后面timeROC函数在运行中途报错。complete.cases这一步很重要——timeROC对缺失值默认是直接报错或者整行剔除与其让它自动丢不如自己在读数据时把处理逻辑写明白知道哪些样本被排除了。参数说明status里事件发生必须记为1删失记为0。如果你的数据反着来AUC会和预期完全相反——所有曲线都跑到对角线下方那是因为事件和删失的编码颠倒了。时间是数值单位混入日期格式如2023-01-01会直接解析失败建议入组前先统一换算成月或天。3. 把脚本跑起来从timeROC核心参数到SCI级别的ROC.pdf3.1 主脚本解读bioR43.multiTimeROC.R做了哪四件事解压看主脚本逻辑很清晰拆开就是四步读取数据、调用timeROC计算多时间点AUC、把多条ROC曲线叠画到一张图、输出PDF。下面这段代码是我按原资源包的功能结构重写的可运行版本保留了核心调用方式你拿到原脚本后对着看结构是一一对应的。# bioR43.multiTimeROC.R 核心逻辑 # 加载依赖包 library(survival) library(timeROC) # 1. 读取示例数据 dat - read.table(input.txt, header TRUE, sep \t) # 确认数据结构 head(dat) # 2. 定义需要评估的时间点单位跟time列一致 # 常见做法取随访时间的四分位数或者根据临床意义选择 eval_times - c(12, 24, 36) # 12/24/36个月 # 3. 计算多时间点AUC res - timeROC( T dat$time, # 生存时间 delta dat$status, # 事件状态1事件0删失 marker dat$marker, # 模型预测值 cause 1, # 竞争风险下指定关注的事件类型 times eval_times, # 评估时间点 iid TRUE, # 是否计算iid分解用于置信区间 weighted TRUE # 是否做IPCW加权生存ROC必须TRUE ) # 4. 查看各时间点AUC和置信区间 print(res$AUC) print(confint(res, level 0.95)$CI_vect) # 5. 生成PDF同图画多条ROC曲线 pdf(ROC.pdf, width 7, height 6) # 用第一个时间点建画布后续时间点叠加lines plot(res, time eval_times[1], col #C00000, lwd 2, title Multi-time ROC Curve) for (i in 2:length(eval_times)) { lines(res, time eval_times[i], col i 1, lwd 2) } # 添加对角线参考 abline(a 0, b 1, lty 2, col gray50) # 添加图例位置根据曲线分布选bottomright或topleft legend(bottomright, legend paste0(AUC at , eval_times, months , round(res$AUC, 3)), col c(#C00000, 3:(length(eval_times) 1)), lwd 2, bty n) dev.off()逻辑说明脚本的骨架是timeROC()函数这一次调用它同时算出了times向量里所有时间点的AUC、标准误、置信区间并保存了画ROC曲线所需的全部坐标数据。iidTRUE是计算渐近分布的标准做法没有它confint()算不出区间估计。weightedTRUE是生存ROC生效的关键等于启用IPCW加权这是和普通ROC最本质的区别。可视化部分先plot()第一个时间点的ROC作为底再用lines()叠加其余时间点的曲线不同颜色区分。参数说明cause1在存在竞争风险事件时指定“关注哪个事件”比如复发和死亡都是终点事件你想单独评估“复发”这个终点就设cause1如果只有单一终点保持默认1即可。eval_times的选择有讲究——不能等距拍脑袋设一般取随访时间的25%、50%、75%分位数或者按临床意义定比如肿瘤研究常取1年、3年、5年对应12/36/60个月。注意times最大值不能超过最大随访时间否则那个时间点的AUC会报错或者返回NA后面避坑章节会细说。3.2 输出物检查怎么判断ROC.pdf是靠谱的还是翻车了跑完代码之后别急着关RStudio先看两个东西。第一是控制台打印的res$AUC第二个是confint输出的置信区间。# 输出多时间点AUC汇总表 auc_summary - data.frame( time_point eval_times, AUC round(res$AUC, 3), lower_95 round(confint(res, level 0.95)$CI_vect[, 1], 3), upper_95 round(confint(res, level 0.95)$CI_vect[, 2], 3) ) print(auc_summary) write.csv(auc_summary, multi-time_ROC_AUC_summary.csv, row.names FALSE)逻辑说明这段是论文方法的补充材料常客——AUC_summary.csv会作为补充表格放进SCI论文附件。confint()返回的是一个矩阵CI_vect是置信区间矩阵行列对应每个时间点取第1列为下限第2列为上限。除了数字还要打开ROC.pdf看图。正常的图应该满足三个条件每条曲线都从左上角高灵敏度高特异度区域出发、整体在对角线上方、不同时间点的曲线能明显区分开且AUC大小关系和临床预期一致。如果看到某条曲线贴着对角线走说明该时间点模型预测能力接近随机猜这是正常的分析结果如实报告就好如果曲线大面积在对角线下方那不是模型差而是代码里status编码反了回到第2章改回来重跑。判断完了再做一步进阶验证画时间- AUC曲线。这是评估模型随时间变化稳定性的常用补充图timeROC包也直接支持。# 连续评估多个整月时间点观察AUC的时间趋势 continuous_times - seq(from 6, to max(eval_times) * 1.2, by 6) res_con - timeROC( T dat$time, delta dat$status, marker dat$marker, cause 1, times continuous_times, iid TRUE, weighted TRUE ) pdf(AUC_over_time.pdf, width 7, height 5) plot(res_con$times, res_con$AUC, type b, pch 19, col #C00000, xlab Time (months), ylab AUC, main AUC over Time, ylim c(0.4, 1)) dev.off()逻辑说明这段把离散的几个时间点扩展成每6个月一个的连续序列得到一条“AUC随时间变化”的折线能看出模型在早期、中期、晚期的区分能力是上升、下降还是平稳。临床上如果AUC随随访时间延长而上升说明模型对远期事件的预测更好如果早期高后期掉得厉害说明模型对近期高风险人群敏感远期能力弱。这个图可以作为正文图2或补充材料提交。参数说明seq(from6, tomax(eval_times)*1.2, by6)里的by6代表每隔6个月计算一次时间单位是月如果你的数据单位是天要改成by180这类对应值。ylim范围视实际情况调整别硬套c(0.4, 1)。4. 避坑指南多时间点生存ROC的六个常见问题4.1 报错“times必须小于最大随访时间”或AUC显示NA现象代码运行到timeROC()时报错提示时间点向量里有值大于最大观察时间或者没报错但某个时间点的AUC返回NaN/NA。原因timeROC要求评估时间点必须落在随访时间范围内。随访时间最长的样本是60个月你设了72个月作为评估点在72个月这个时点全队列里没有任何人达到“观察满72个月且未删失”的状态权重估计失去支撑结果自然无法计算。解决先跑summary(dat$time)看最大随访时间以最大随访时间tmax为上限。时间点不要正好卡在最大值上留出一点余量。# 自查最大随访时间以此为上限设置评估点 t_max - max(dat$time, na.rm TRUE) # 合理做法取t_max的70%-90%作为最晚评估点 eval_times - quantile(dat$time, probs c(0.25, 0.5, 0.75), na.rm TRUE) eval_times - round(eval_times, 1) # 取整方便报告 eval_times - eval_times[eval_times t_max * 0.9] # 剔除接近最大值的点 print(eval_times)逻辑说明以分位数作为评估点是多个R包教程里的通行做法能避免人为拍脑袋选的随意性也保证每个时间点都有足够样本支撑。剔除接近tmax的点是为了防止删失比例过高时AUC估计方差过大——时间点越靠后剩余风险集越小置信区间越宽画出来的图说服力越弱。4.2 报错“delta必须为数值向量”或曲线全部在对角线下方现象timeROC()报错提示delta类型不对或者图能出来但曲线都在对角线下面。原因status列在R里被读成了字符型或因子型比如1、0带引号或者用了Dead/Alive。字符被强制转数值时会出错。曲线翻转到对角线下方这事原因基本只有一种status编码记反了——事件0删失1。解决读数据时用stringsAsFactorsFALSE再做一步显式转换加检查。# 显式把status转成数值型0/1 dat$status - as.numeric(as.character(dat$status)) stopifnot(all(dat$status %in% c(0, 1))) # 不是0/1组合直接中断 # 交叉验证事件样本的time应该更短才算合理 aggregate(time ~ status, data dat, FUN median)逻辑说明stopifnot是安全的强制检查口条件不满足直接抛出异常避免带着脏数据往下跑。aggregate那行是验证编码方向是否符合临床直觉——按常理事件组的生存时间中位数应该显著短于删失组如果反了八成是status定义反了。4.3 加载timeROC包报错“程辑包不存在”或依赖缺失现象library(timeROC)提示不存在这个包或者install.packages(timeROC)报了“无法连接到CRAN仓库”“依赖包不存在”之类的问题。原因要么是R版本太老或镜像源没配置好要么是缺少底层依赖包它依赖survival和一些基础工具包。R包的安装报错信息经常不直接说明缺的是谁只给一个像“不存在叫‘xx程辑包’这个名字”的提示排查起来确实玄学。解决先设CRAN镜像再装依赖最后装timeROC。别直接裸敲install.packages(timeROC)。# 设置国内镜像兼容部分网络环境下CRAN连接不稳定的情况 options(repos c(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN/)) # 先确保依赖包就位 install.packages(c(survival, pROC)) # pROC非必需但画图时常用 # 再装timeROC install.packages(timeROC) # 装完验证能否正常加载 library(timeROC)逻辑说明R包安装报错90%是镜像源问题和依赖缺失先设镜像再装依赖是通用解法。pROC顺手装上是储备后续如果要画传统二分类ROC对比不用再装一次。4.4 图例把ROC曲线压住了或不显示AUC数值现象图明显比例失衡图例框太大把右下角的ROC曲线遮住或者图例的AUC数字远到看不见。原因legend的位置参数设成了固定的bottomright但实际曲线恰好也聚集在右下角。生存ROC曲线常见形态是陡峭的左上角加平缓的右下拐右下角通常有曲线穿过图例一放就压线。解决看图例位置改成topleft如果还挡手动用inset微调配合segments加AUC标签。# 动态判断图例放哪边以第一个时间点的曲线坐标为参考 library(prodlim) # timeROC依赖包内含plot方法 legend_pos - ifelse(mean(res$TP[1, ]) 0.5, bottomright, topleft) legend(legend_pos, legend paste0(t , eval_times, m, AUC , round(res$AUC, 2)), col c(#C00000, 3:(length(eval_times) 1)), lwd 2, bty n, cex 0.9)逻辑说明res$TP是timeROC结果里的灵敏度矩阵第一行对应第一个时间点。用它在图上的平均高度判断曲线偏向图的上半区还是下半区决定图例落在哪个角。cex0.9缩小字号也能腾出空间。btyn去掉图例边框整体更干净这在SCI图上更常见。4.5 多时间点曲线颜色太接近打印后分不清现象颜色是默认的2,3,4号色在屏幕上能看清但论文打印成黑白稿后3号绿色和4号蓝色几乎无法区分。解决自主定义一组色彩对比度高、色盲友好的配色并在画图时用col_vec统一管理。# 色盲友好的配色方案按时间点顺序排列 col_vec - c(#C00000, #0072B2, #E69F00, #009E73, #999999) # 使用时截取前n个颜色 plot_cols - col_vec[1:length(eval_times)]逻辑说明#C00000是深红、#0072B2是中蓝、#E69F00是橙色这组方案在多本可视化书籍里被验证过在色盲人群中和黑白打印下都能区分。4.6 变量列名带空格或特殊字符脚本运行到一半报错现象自己导入Excel导出的CSV列名是OS time月、Event(死亡1)这种带空格、中文、括号的格式脚本里用dat$OS time取列报语法错误或者timeROC()提示找不到列。原因R中列名含空格必须用反引号包裹访问脚本是按input.txt的简洁列名写的拿到带复杂列名的表直接套用必然翻车。解决读取后统一改列名。dat - read.csv(your_data.csv, check.names FALSE) # 保留原列名以防中文乱码 names(dat) - c(time, status, marker) # 按位置重命名前三列逻辑说明check.namesFALSE防止R自动把列名里的空格改成语点但用完后立刻手动重命名一举两得既避免原列名读取时的错位又让后续代码能直接访问dat$time。5. 进阶玩法把这份源代码改造成批量分析工具5.1 批量处理多个模型ROC全家桶一键生成你手里大概率不止一个模型要验证。比如你做了COX回归风险评分、Lasso筛选后的signature、还有单个biomarker表达量理论上每组都要出一套ROC图。目标是把bioR43.multiTimeROC.R从“单个模型专用”改造成“多模型通吃”。# 封装成函数一个模型一套结果 run_multitime_roc - function(data, marker_col, time_col time, status_col status, eval_times c(12, 24, 36), pdf_prefix ROC) { res - timeROC( T data[[time_col]], delta data[[status_col]], marker data[[marker_col]], cause 1, times eval_times, iid TRUE, weighted TRUE ) # 输出单模型AUC表 ci - confint(res, level 0.95)$CI_vect out - data.frame( marker marker_col, time eval_times, auc round(res$AUC, 3), lower round(ci[, 1], 3), upper round(ci[, 2], 3) ) # 画图 pdf(paste0(pdf_prefix, _, marker_col, .pdf), width 6.5, height 5.5) plot(res, time eval_times[1], col #C00000, lwd 2, title paste0(ROC: , marker_col)) for (i in 2:length(eval_times)) { lines(res, time eval_times[i], col c(#0072B2, #E69F00)[i - 1], lwd 2) } abline(a 0, b 1, lty 2, col gray50) legend(bottomright, legend paste0(eval_times, m AUC, round(res$AUC, 2)), col c(#C00000, #0072B2, #E69F00)[1:length(eval_times)], lwd 2, bty n) dev.off() return(out) } # 对三个候选模型批量执行 models - c(risk_score, lncRNA_sig, TNMB_stage) results - do.call(rbind, lapply(models, function(m) { run_multitime_roc(dat, marker_col m) })) print(results)逻辑说明函数里data[[marker_col]]用双中括号取动态列名这是批量处理的关键写法。do.call(rbind, lapply(...))是R里批量处理返回数据框的标准组合拳能把多个模型的AUC汇总成一张表。pdf_prefix参数控制输出文件名避免多个模型覆盖同一个PDF。参数说明marker_col可以传入任意列名包括基因表达量、风险评分、分期分组只要转成数值。组别数据如TNM分期可以转成有序数值1/2/3/4当marker用逻辑上是把分期当作等级变量评估其预测能力。5.2 把AUC对比结果写进论文参数和方法描述的参考写法做完分析总要写到论文里方法学部分和结果部分的写法是有套路的。我拆过的项目里评审对ROC分析最看重的就是三点时间点怎么选的、删失怎么处理的、置信区间有没有。下面是一个可以直接改的模板。多时间点ROC曲线用于评估各模型在不同随访时间点的预测区分度。评估时间点设定为随访时间的25%12个月、50%24个月和75%36个月分位数。采用timeROC包基于逆删失概率加权IPCW估计各时间点的AUC及其95%置信区间。结果显示风险评分在12、24、36个月的AUC分别为0.8295%CI: 0.75-0.89、0.7895%CI: 0.71-0.85和0.7495%CI: 0.66-0.82区分度随时间略有下降。这里AUC数值是我编的示意你替换成自己的真实结果即可。注意写法要点时间点的选择依据必须写清楚分位数临床节点一句话点出IPCW加权AUC结果按时间顺序报同时给置信区间最后加一句时间趋势描述。这个模板能应对大多数临床SCI期刊的方法学描述要求。再补充一个技巧如果三个模型在同一队列上比较光看各模型的AUC不够最好把三条曲线画在同一张图上做视觉比较。做法是改一下第3章的plot部分多循环一层模型。# 多模型在同一图上的比较 pdf(multi_model_compare.pdf, width 7, height 6) plot(res_1, time 24, col #C00000, lwd 2, title Model Comparison at 24 Months) lines(res_2, time 24, col #0072B2, lwd 2) lines(res_3, time 24, col #E69F00, lwd 2) legend(bottomright, legend c(Risk score (AUC0.78), lncRNA sig (AUC0.83), TNM stage (AUC0.65)), col c(#C00000, #0072B2, #E69F00), lwd 2, bty n) dev.off()逻辑说明这个图的信息密度比多时间点曲线更高——同一时间点下多模型的区分能力一目了然审稿人不需要自己去对比三个子图。实际操作中做多模型比较时固定时间点为中位随访时间或临床最关注的节点如肿瘤研究的3年或5年生存率节点更有说服力。参数说明比较图里res_1、res_2、res_3是分别用timeROC()对三个模型在同一个24个月时间点算出的结果对象。运行前务必确保三个模型的marker列和time/status列都对齐数据框排序一致避免错配。5.3 一个经验教训AUC是点估计置信区间才是真相拆这个资源包的过程中我用模拟数据做过一次测试把样本量砍到100例分组做多时间点ROC某个时间点的AUC点估计是0.72但95%置信区间是0.51到0.93——一个区间横跨了0.5这条“随机猜测”线。如果不看置信区间单报0.72论文写出来是有隐患的。从那以后我每次跑多时间点ROC都会强制走一遍先summary确认时间范围再用confint核对各时间点的置信区间最后才画图写结果。这份源代码帮你省掉了从零写包的工作量但统计报告的底线仍然在自己手里希望帮到你。本文还有配套的精品资源点击获取
