简介这份R语言源代码面向需要绘制SCI科研生存ROC曲线的研究者尤其适用于具备一定R基础、正在准备医学论文或课题展示的科研人员。资源围绕生存数据的时间依赖ROC分析展开使用timeROC等包完成ROC计算、AUC提取、最佳阈值选择以及图形美化最终输出可直接用于论文发表的生存ROC曲线图。压缩包共包含3个文件涵盖R脚本、txt数据文件和PDF说明文档整体仅11KB结构简洁、注释清晰便于对照学习和按需调整。目前已有433人学习/下载适合在诊断效能评估、预后模型验证等场景中快速应用。通过运行和解读这份代码用户不仅能理解生存ROC与普通ROC的差异还能掌握敏感度、特异性等核心指标的计算与可视化方法且其中涉及的阈值优化和图形定制思路也可迁移到其他分类模型评价任务中是一份实用价值较高的R语言绘图参考资料。1. 生存ROC曲线不是普通ROC加一条时间轴先把这份源代码的定位看清楚拿到「R语言绘制SCI科研生存ROC曲线源代码.zip」这个资源时多数人第一反应是我已经会用pROC画ROC了生存ROC不过是把结局换成生存状态。实际跑一遍你会发现pROC包在带删失的生存数据上根本算不对——它要求每个样本都有一个固定的二分类标签而生存数据里有一部分病人在随访结束时仍未发生事件这部分“删失”样本既不能算阳性也不能算阴性但你又不能直接丢掉。这正是时间依赖ROCtime-dependent ROC存在的意义。这个zip里的bioR42.timeROC.R、input.txt和ROC.pdf构成一个最小可运行闭环核心是用timeROC包计算不同随访时间点上的AUC并把结果输出为符合SCI投稿基本要求的曲线图。适合已经会用R读数据、装包但对生存分析中的时间依赖评价指标还不够熟的科研人员。下面我会把原理、代码、数据格式和出图参数逐一拆开照着改就能用在你的队列数据上。2. 时间依赖ROC的原理与timeROC包选型AUC随随访时间变化2.1 为什么pROC和ROCR画不了生存数据普通ROC曲线要求每个样本有一个确定的状态标签事件发生或未发生。roc()函数内部会把状态因子化计算所有可能阈值下的真阳性率和假阳性率。但在随访数据中假设随访截止时间是第42个月有的人在第30个月死亡有的人第30个月失访——失访样本在第42个月时是“未发生”还是“已经发生”你只知道他最后被观察到时还活着。直接标记为“未发生”会低估事件率标记为“发生”则是错误。这个信息缺失是删失censoring而生存ROC需要把时间维度和删失分布同时放进去不能简单用pROC::roc(status, predicted)。很多初学者会过滤掉删失样本再画ROC这是最危险的误用。删失不是随机丢失它常与病人的基础风险相关过滤后AUC被系统性地高估或低估。timeROC的核心思路是在每个目标时间点t对仍然活着且未删失的样本用Kaplan-Meier或条件概率方法估计该时间点的灵敏度和特异度这样既用了删失样本的随访信息又避免了硬性分类。2.2 timeROC如何定义灵敏度和特异度假设我们要看第42个月文件里的“42”很可能指42个月随访点的预测区分能力。记病人的观测时间time事件指示delta1发生目标事件0删失或未发生预测标记marker例如风险评分或生物标志物。对于时间点ttimeROC包按以下定义计算灵敏度 P(marker c | 在t时刻前发生事件) 特异度 P(marker ≤ c | 在t时刻时尚未发生事件)注意第二个概率的估计要处理删失。timeROC包默认使用“marginal”加权法即用Kaplan-Meier估计删失分布对每个样本赋权重。如果你设weightingmarginal代码会用事件时间分布去校正那些随访时间不足t但未发生事件的样本。另一种是weightingcox利用Cox模型估计删失分布当标记与删失相关时更稳健但计算量也更大。2.3 选timeROC而不是survivalROC的理由做时间依赖ROCR里还有survivalROC包、timeROC包以及SurvROC等。我选择timeROC是因为它在这类需求中的覆盖更完整输出AUC的逐点置信区间并且支持iid自助法直接算标准误能在同一个plot里叠加多个时间点的ROC曲线方便表达“模型在早期预测好还是晚期预测好”对cause1的多竞争风险场景也能处理计算速度在几千样本量下完全可以接受且依赖包少。survivalROC只能给出单个时间点的结果画多个时间点要写循环而且置信区间需要额外嵌入。对于一篇需要展示35个随访时间点的SCI论文timeROC的封装更直接。下面的第4章会展示这份zip里的代码到底怎么落地。3. input.txt的数据组织与清洗把临床随访表变成R可读的输入3.1 列字段定义time、status、marker打开input.txt你大概率看到的是一个以制表符分隔的文本文件。为了配合bioR42.timeROC.R数据文件至少要包含三列观测时间、删失状态、预测标记。源码中对应的变量名通常是time、status和marker但如果你拿到的是自己项目的随访表列名可能是OS.time、OS.event、riskScore需要在读入后重新映射。时间单位要统一。如果文件里写的是天数而你要看第42个月的ROC就需要先dat$time - dat$time / 30.44。我建议在清洗阶段就把时间转换为“月”或“年”避免后面在times参数里搞混。3.2 用read.table载入并检查数据先用基础命令读入并做结构检查# 读入txt分隔符用\t如果出现乱码检查fileEncoding dat - read.table(input.txt, header TRUE, sep \t, stringsAsFactors FALSE) # 检查列名和类型确认没有把数值读成字符 str(dat) # 对关键三列做缺失值和范围检查 summary(dat$time) table(dat$status, useNA ifany) range(dat$marker, na.rm TRUE)str(dat)会显示每一列的类类型如果time列被读成了chr大概率是文件里有空格或引号需要检查原始数据。table(dat$status)这一步尤其重要统计1和0的数量确认事件发生比例是否过低。生存ROC在事件数少于30个时曲线会非常不稳定即使程序能跑出AUC审稿人也会质疑。3.3 状态变量的编码约定timeROC要求delta是0/1数值1表示发生了目标事件cause0表示删失。如果你的原始编码是“Dead/Alive”或“1/2”要转换成0/1。这里有个常见误用把“未发生但随访满WHOLE随访期”的样本也设为1。例如第42个月评估一个人在随访第50个月死亡在第42个月时他还没发生事件应该算作“尚未发生”而不是事件。正确的做法是在构造时间时取“死亡时间与评估时间的最小值”状态同步设置dat$eval_time - 42 dat$time - pmin(dat$time, dat$eval_time) dat$status - ifelse(dat$time dat$eval_time dat$event 0, 0, dat$event)只有明确在42个月之前发生事件的样本才应在第42个月的ROC中作为阳性参与计算。timeROC包的T参数直接用原始时间它会内部处理“时间大于t但未删失”的情况所以只要你传入的是真实随访时间而不是截断后的时间包会正确处理。但如果你自己先截断了反而会破坏删失结构。4. bioR42.timeROC.R运行拆解从timeROC()到ROC.pdf的完整流程4.1 核心计算代码注解项目里的bioR42.timeROC.R核心计算部分大致等价于下面的代码。实际运行时你只需要把dat换成自己的数据library(timeROC) # time-dependent ROC的主包 library(survival) # 用于后续KM估计和其他生存函数 # 读入项目自带的数据 dat - read.table(input.txt, header TRUE, sep \t) # 计算多个时间点的ROC这里以12,24,36,42个月为例 result - timeROC( T dat$time, # 观测时间单位需要与times一致 delta dat$status, # 事件指示1事件0删失 marker dat$marker, # 预测标记越高表示风险越大 cause 1, # 感兴趣的事件编码 times c(12, 24, 36, 42), # 需要评估的随访时间点 iid TRUE, # 计算iid加权输出置信区间 weighting marginal # 删失分布加权方法 ) # 查看各时间点的AUC print(result$AUC)timeROC()的返回值是一个列表其中AUC是各时间点的AUC向量。如果iidTRUE还会包含inference和conf.int字段。关于weighting参数marginal假设删失独立于标记cox则建模标记对删失的影响。如果你的数据里删失比例高40%且标记与治疗方案或随访强度相关我一般会切换到coxresult_cox - timeROC( T dat$time, delta dat$status, marker dat$marker, cause 1, times c(12, 24, 36, 42), iid TRUE, weighting cox )注意cox要求survival包已加载并且计算时间会明显增加。如果两次加权结果差异很大说明删失机制不可忽略需要用cox方案并且报告中说明。4.2 关键参数times、cause、iid的边界times的取值范围必须在随访时间内。比如随访最长是60个月你写timesc(12,24,36,42)没问题但如果某个时间点超过最大随访时间timeROC会报错或给出无效AUC。建议先用quantile(dat$time[dat$status1], c(0.25,0.5,0.75))看事件时间的四分位数然后再决定展示哪些时间点。cause参数在多竞争风险场景下很重要。如果状态变量里1癌症相关死亡2其他原因死亡你想评估“癌症相关死亡”的预测能力就设cause1另一事件自动作为删失处理。但你必须在状态变量中保持数字编码不能用因子。iid参数控制是否计算影响函数只有设为TRUE时后面的conf.int才有效。在你只是快速预览AUC趋势时可以iidFALSE速度提升不少写论文时则必须iidTRUE。4.3 绘图与PDF输出plot.timeROC的参数细节项目里的ROC.pdf是由下面的绘图逻辑生成的。timeROC自带plot方法直接传入result和目标时间点# 输出PDF文件尺寸一般设定为7x7英寸符合SCI投稿要求 pdf(ROC.pdf, width 7, height 7) # 绘制第42个月的ROC曲线 plot( result, time 42, col #E6194B, # 曲线颜色建议用高对比色 lwd 2.5, # 线宽 xlab 1 - Specificity, ylab Sensitivity, main Time-dependent ROC at Month 42, cex.main 1.5, cex.lab 1.2 ) # 画对角线 abline(a 0, b 1, lty 2, col gray50, lwd 1.5) # 将AUC值添加到图例 legend( bottomright, legend sprintf(AUC %.3f, result$AUC[t42]), col #E6194B, lwd 2.5, bty n ) dev.off()plot.timeROC的time参数必须是times中出现过的值。绘制多个时间点时你可以用同一个result多次调用plot用addTRUE叠加但要注意plot方法默认会新建画布需要先用plot(...)初始化一个时间点再依次调用addTRUE。如果你希望所有时间点用不同颜色可以这样plot(result, time 12, col #440154, lwd 2, ...) plot(result, time 24, col #21918c, lwd 2, add TRUE) plot(result, time 36, col #fde725, lwd 2, add TRUE)这里有一个容易忽略的坑lwd和col在plot.timeROC中不能通过贝塞尔曲线外的列表参数统一设置每次addTRUE都需要重新指定否则会沿用第一次设置。另外xlab和ylab只在第一次调用时生效后续叠加调用中即便写上也会被忽略。5. 把生存ROC图发到SCI曲线美化、置信区间与常见报错5.1 用ggplot2重绘生存ROC曲线timeROC的base R绘图风格比较朴素如果期刊要求矢量图且主题统一我会选择把计算出的坐标点抓出来重新用ggplot2绘制。核心是先手动提取特异度和灵敏度向量library(ggplot2) # 提取第42个月ROC曲线的坐标点 roc_df - data.frame( fpr 1 - result$Sp$42, # 假阳性率 tpr result$Se$42 # 真阳性率 ) ggplot(roc_df, aes(x fpr, y tpr)) geom_abline(intercept 0, slope 1, linetype dashed, color gray50) geom_line(color #E6194B, linewidth 1.2) coord_equal() labs( title Time-dependent ROC at Month 42, x 1 - Specificity, y Sensitivity ) annotate(text, x 0.7, y 0.2, label paste0(AUC , round(result$AUC[4], 3))) theme_minimal(base_size 14)这里result$Sp和result$Se是timeROC返回值中的两个列表键名是时间点对应的字符串必须用42反引号或as.character(42)索引。有些版本把键名存为数字你需要先names(result$Se)查看再决定怎么写。用coord_equal()保证横纵轴比例一致否则曲线看起来会失真。5.2 用自助法计算95%置信区间timeROC包内置的confint方法基于iid影响函数比较快但不是百分位自助法。如果你想用更贴近临床论文风格的结果可以在外围写一个简单的自助循环set.seed(123) boot_auc - numeric(500) n - nrow(dat) for (b in 1:500) { idx - sample(n, n, replace TRUE) boot_dat - dat[idx, ] res_boot - timeROC( T boot_dat$time, delta boot_dat$status, marker boot_dat$marker, cause 1, times 42, iid FALSE ) boot_auc[b] - res_boot$AUC[1] } # 输出2.5%和97.5%分位数 quantile(boot_auc, c(0.025, 0.975))自助法每次都要重算timeROC500次在中等数据量下约一二十秒。注意自助样本中如果某次抽样恰好没有发生事件timeROC会报错需要在循环内加tryCatch跳过无效迭代res_boot - tryCatch( timeROC(...), error function(e) return(NULL) ) if (is.null(res_boot)) next最后把分位数写到图片的图例或表格中。对于SCI文章审稿人更看重置信区间而非仅仅一个点估计。5.3 我在实际项目中踩过的三个坑第一时间单位不一致。如果input.txt里的时间是天数而times里写的是30、60AUC还算得出来但曲线对应的临床解释会变得很奇怪。我习惯在脚本开头统一除以30.44转为月再在times里写月数。第二status列包含2。多状态数据里如果没排除Competing RisktimeROC会把2当作删失处理但如果2是0以外的数字最好先转化成二分类。即使目标是“全因死亡”也要确保目录里没有其他编码。第三marker的方向搞反。timeROC默认marker越大风险越高如果你的marker是保护因子越大越安全得到的AUC会小于0.5。解决办法是传入-marker或调整代码。验证方法很简单在R中运行cor(dat$marker, dat$time[dat$status1], methodspearman)如果负相关且AUC小于0.5就取负号。最后在输出PDF前用dev.off()收尾否则ROC.pdf文件中看不到图形如果你在RStudio里反复运行脚本建议每次绘制前先graphics.off()清空画布避免叠加残留。本文还有配套的精品资源点击获取
R语言生存ROC曲线绘制：timeROC包原理、代码与SCI论文实战
简介这份R语言源代码面向需要绘制SCI科研生存ROC曲线的研究者尤其适用于具备一定R基础、正在准备医学论文或课题展示的科研人员。资源围绕生存数据的时间依赖ROC分析展开使用timeROC等包完成ROC计算、AUC提取、最佳阈值选择以及图形美化最终输出可直接用于论文发表的生存ROC曲线图。压缩包共包含3个文件涵盖R脚本、txt数据文件和PDF说明文档整体仅11KB结构简洁、注释清晰便于对照学习和按需调整。目前已有433人学习/下载适合在诊断效能评估、预后模型验证等场景中快速应用。通过运行和解读这份代码用户不仅能理解生存ROC与普通ROC的差异还能掌握敏感度、特异性等核心指标的计算与可视化方法且其中涉及的阈值优化和图形定制思路也可迁移到其他分类模型评价任务中是一份实用价值较高的R语言绘图参考资料。1. 生存ROC曲线不是普通ROC加一条时间轴先把这份源代码的定位看清楚拿到「R语言绘制SCI科研生存ROC曲线源代码.zip」这个资源时多数人第一反应是我已经会用pROC画ROC了生存ROC不过是把结局换成生存状态。实际跑一遍你会发现pROC包在带删失的生存数据上根本算不对——它要求每个样本都有一个固定的二分类标签而生存数据里有一部分病人在随访结束时仍未发生事件这部分“删失”样本既不能算阳性也不能算阴性但你又不能直接丢掉。这正是时间依赖ROCtime-dependent ROC存在的意义。这个zip里的bioR42.timeROC.R、input.txt和ROC.pdf构成一个最小可运行闭环核心是用timeROC包计算不同随访时间点上的AUC并把结果输出为符合SCI投稿基本要求的曲线图。适合已经会用R读数据、装包但对生存分析中的时间依赖评价指标还不够熟的科研人员。下面我会把原理、代码、数据格式和出图参数逐一拆开照着改就能用在你的队列数据上。2. 时间依赖ROC的原理与timeROC包选型AUC随随访时间变化2.1 为什么pROC和ROCR画不了生存数据普通ROC曲线要求每个样本有一个确定的状态标签事件发生或未发生。roc()函数内部会把状态因子化计算所有可能阈值下的真阳性率和假阳性率。但在随访数据中假设随访截止时间是第42个月有的人在第30个月死亡有的人第30个月失访——失访样本在第42个月时是“未发生”还是“已经发生”你只知道他最后被观察到时还活着。直接标记为“未发生”会低估事件率标记为“发生”则是错误。这个信息缺失是删失censoring而生存ROC需要把时间维度和删失分布同时放进去不能简单用pROC::roc(status, predicted)。很多初学者会过滤掉删失样本再画ROC这是最危险的误用。删失不是随机丢失它常与病人的基础风险相关过滤后AUC被系统性地高估或低估。timeROC的核心思路是在每个目标时间点t对仍然活着且未删失的样本用Kaplan-Meier或条件概率方法估计该时间点的灵敏度和特异度这样既用了删失样本的随访信息又避免了硬性分类。2.2 timeROC如何定义灵敏度和特异度假设我们要看第42个月文件里的“42”很可能指42个月随访点的预测区分能力。记病人的观测时间time事件指示delta1发生目标事件0删失或未发生预测标记marker例如风险评分或生物标志物。对于时间点ttimeROC包按以下定义计算灵敏度 P(marker c | 在t时刻前发生事件) 特异度 P(marker ≤ c | 在t时刻时尚未发生事件)注意第二个概率的估计要处理删失。timeROC包默认使用“marginal”加权法即用Kaplan-Meier估计删失分布对每个样本赋权重。如果你设weightingmarginal代码会用事件时间分布去校正那些随访时间不足t但未发生事件的样本。另一种是weightingcox利用Cox模型估计删失分布当标记与删失相关时更稳健但计算量也更大。2.3 选timeROC而不是survivalROC的理由做时间依赖ROCR里还有survivalROC包、timeROC包以及SurvROC等。我选择timeROC是因为它在这类需求中的覆盖更完整输出AUC的逐点置信区间并且支持iid自助法直接算标准误能在同一个plot里叠加多个时间点的ROC曲线方便表达“模型在早期预测好还是晚期预测好”对cause1的多竞争风险场景也能处理计算速度在几千样本量下完全可以接受且依赖包少。survivalROC只能给出单个时间点的结果画多个时间点要写循环而且置信区间需要额外嵌入。对于一篇需要展示35个随访时间点的SCI论文timeROC的封装更直接。下面的第4章会展示这份zip里的代码到底怎么落地。3. input.txt的数据组织与清洗把临床随访表变成R可读的输入3.1 列字段定义time、status、marker打开input.txt你大概率看到的是一个以制表符分隔的文本文件。为了配合bioR42.timeROC.R数据文件至少要包含三列观测时间、删失状态、预测标记。源码中对应的变量名通常是time、status和marker但如果你拿到的是自己项目的随访表列名可能是OS.time、OS.event、riskScore需要在读入后重新映射。时间单位要统一。如果文件里写的是天数而你要看第42个月的ROC就需要先dat$time - dat$time / 30.44。我建议在清洗阶段就把时间转换为“月”或“年”避免后面在times参数里搞混。3.2 用read.table载入并检查数据先用基础命令读入并做结构检查# 读入txt分隔符用\t如果出现乱码检查fileEncoding dat - read.table(input.txt, header TRUE, sep \t, stringsAsFactors FALSE) # 检查列名和类型确认没有把数值读成字符 str(dat) # 对关键三列做缺失值和范围检查 summary(dat$time) table(dat$status, useNA ifany) range(dat$marker, na.rm TRUE)str(dat)会显示每一列的类类型如果time列被读成了chr大概率是文件里有空格或引号需要检查原始数据。table(dat$status)这一步尤其重要统计1和0的数量确认事件发生比例是否过低。生存ROC在事件数少于30个时曲线会非常不稳定即使程序能跑出AUC审稿人也会质疑。3.3 状态变量的编码约定timeROC要求delta是0/1数值1表示发生了目标事件cause0表示删失。如果你的原始编码是“Dead/Alive”或“1/2”要转换成0/1。这里有个常见误用把“未发生但随访满WHOLE随访期”的样本也设为1。例如第42个月评估一个人在随访第50个月死亡在第42个月时他还没发生事件应该算作“尚未发生”而不是事件。正确的做法是在构造时间时取“死亡时间与评估时间的最小值”状态同步设置dat$eval_time - 42 dat$time - pmin(dat$time, dat$eval_time) dat$status - ifelse(dat$time dat$eval_time dat$event 0, 0, dat$event)只有明确在42个月之前发生事件的样本才应在第42个月的ROC中作为阳性参与计算。timeROC包的T参数直接用原始时间它会内部处理“时间大于t但未删失”的情况所以只要你传入的是真实随访时间而不是截断后的时间包会正确处理。但如果你自己先截断了反而会破坏删失结构。4. bioR42.timeROC.R运行拆解从timeROC()到ROC.pdf的完整流程4.1 核心计算代码注解项目里的bioR42.timeROC.R核心计算部分大致等价于下面的代码。实际运行时你只需要把dat换成自己的数据library(timeROC) # time-dependent ROC的主包 library(survival) # 用于后续KM估计和其他生存函数 # 读入项目自带的数据 dat - read.table(input.txt, header TRUE, sep \t) # 计算多个时间点的ROC这里以12,24,36,42个月为例 result - timeROC( T dat$time, # 观测时间单位需要与times一致 delta dat$status, # 事件指示1事件0删失 marker dat$marker, # 预测标记越高表示风险越大 cause 1, # 感兴趣的事件编码 times c(12, 24, 36, 42), # 需要评估的随访时间点 iid TRUE, # 计算iid加权输出置信区间 weighting marginal # 删失分布加权方法 ) # 查看各时间点的AUC print(result$AUC)timeROC()的返回值是一个列表其中AUC是各时间点的AUC向量。如果iidTRUE还会包含inference和conf.int字段。关于weighting参数marginal假设删失独立于标记cox则建模标记对删失的影响。如果你的数据里删失比例高40%且标记与治疗方案或随访强度相关我一般会切换到coxresult_cox - timeROC( T dat$time, delta dat$status, marker dat$marker, cause 1, times c(12, 24, 36, 42), iid TRUE, weighting cox )注意cox要求survival包已加载并且计算时间会明显增加。如果两次加权结果差异很大说明删失机制不可忽略需要用cox方案并且报告中说明。4.2 关键参数times、cause、iid的边界times的取值范围必须在随访时间内。比如随访最长是60个月你写timesc(12,24,36,42)没问题但如果某个时间点超过最大随访时间timeROC会报错或给出无效AUC。建议先用quantile(dat$time[dat$status1], c(0.25,0.5,0.75))看事件时间的四分位数然后再决定展示哪些时间点。cause参数在多竞争风险场景下很重要。如果状态变量里1癌症相关死亡2其他原因死亡你想评估“癌症相关死亡”的预测能力就设cause1另一事件自动作为删失处理。但你必须在状态变量中保持数字编码不能用因子。iid参数控制是否计算影响函数只有设为TRUE时后面的conf.int才有效。在你只是快速预览AUC趋势时可以iidFALSE速度提升不少写论文时则必须iidTRUE。4.3 绘图与PDF输出plot.timeROC的参数细节项目里的ROC.pdf是由下面的绘图逻辑生成的。timeROC自带plot方法直接传入result和目标时间点# 输出PDF文件尺寸一般设定为7x7英寸符合SCI投稿要求 pdf(ROC.pdf, width 7, height 7) # 绘制第42个月的ROC曲线 plot( result, time 42, col #E6194B, # 曲线颜色建议用高对比色 lwd 2.5, # 线宽 xlab 1 - Specificity, ylab Sensitivity, main Time-dependent ROC at Month 42, cex.main 1.5, cex.lab 1.2 ) # 画对角线 abline(a 0, b 1, lty 2, col gray50, lwd 1.5) # 将AUC值添加到图例 legend( bottomright, legend sprintf(AUC %.3f, result$AUC[t42]), col #E6194B, lwd 2.5, bty n ) dev.off()plot.timeROC的time参数必须是times中出现过的值。绘制多个时间点时你可以用同一个result多次调用plot用addTRUE叠加但要注意plot方法默认会新建画布需要先用plot(...)初始化一个时间点再依次调用addTRUE。如果你希望所有时间点用不同颜色可以这样plot(result, time 12, col #440154, lwd 2, ...) plot(result, time 24, col #21918c, lwd 2, add TRUE) plot(result, time 36, col #fde725, lwd 2, add TRUE)这里有一个容易忽略的坑lwd和col在plot.timeROC中不能通过贝塞尔曲线外的列表参数统一设置每次addTRUE都需要重新指定否则会沿用第一次设置。另外xlab和ylab只在第一次调用时生效后续叠加调用中即便写上也会被忽略。5. 把生存ROC图发到SCI曲线美化、置信区间与常见报错5.1 用ggplot2重绘生存ROC曲线timeROC的base R绘图风格比较朴素如果期刊要求矢量图且主题统一我会选择把计算出的坐标点抓出来重新用ggplot2绘制。核心是先手动提取特异度和灵敏度向量library(ggplot2) # 提取第42个月ROC曲线的坐标点 roc_df - data.frame( fpr 1 - result$Sp$42, # 假阳性率 tpr result$Se$42 # 真阳性率 ) ggplot(roc_df, aes(x fpr, y tpr)) geom_abline(intercept 0, slope 1, linetype dashed, color gray50) geom_line(color #E6194B, linewidth 1.2) coord_equal() labs( title Time-dependent ROC at Month 42, x 1 - Specificity, y Sensitivity ) annotate(text, x 0.7, y 0.2, label paste0(AUC , round(result$AUC[4], 3))) theme_minimal(base_size 14)这里result$Sp和result$Se是timeROC返回值中的两个列表键名是时间点对应的字符串必须用42反引号或as.character(42)索引。有些版本把键名存为数字你需要先names(result$Se)查看再决定怎么写。用coord_equal()保证横纵轴比例一致否则曲线看起来会失真。5.2 用自助法计算95%置信区间timeROC包内置的confint方法基于iid影响函数比较快但不是百分位自助法。如果你想用更贴近临床论文风格的结果可以在外围写一个简单的自助循环set.seed(123) boot_auc - numeric(500) n - nrow(dat) for (b in 1:500) { idx - sample(n, n, replace TRUE) boot_dat - dat[idx, ] res_boot - timeROC( T boot_dat$time, delta boot_dat$status, marker boot_dat$marker, cause 1, times 42, iid FALSE ) boot_auc[b] - res_boot$AUC[1] } # 输出2.5%和97.5%分位数 quantile(boot_auc, c(0.025, 0.975))自助法每次都要重算timeROC500次在中等数据量下约一二十秒。注意自助样本中如果某次抽样恰好没有发生事件timeROC会报错需要在循环内加tryCatch跳过无效迭代res_boot - tryCatch( timeROC(...), error function(e) return(NULL) ) if (is.null(res_boot)) next最后把分位数写到图片的图例或表格中。对于SCI文章审稿人更看重置信区间而非仅仅一个点估计。5.3 我在实际项目中踩过的三个坑第一时间单位不一致。如果input.txt里的时间是天数而times里写的是30、60AUC还算得出来但曲线对应的临床解释会变得很奇怪。我习惯在脚本开头统一除以30.44转为月再在times里写月数。第二status列包含2。多状态数据里如果没排除Competing RisktimeROC会把2当作删失处理但如果2是0以外的数字最好先转化成二分类。即使目标是“全因死亡”也要确保目录里没有其他编码。第三marker的方向搞反。timeROC默认marker越大风险越高如果你的marker是保护因子越大越安全得到的AUC会小于0.5。解决办法是传入-marker或调整代码。验证方法很简单在R中运行cor(dat$marker, dat$time[dat$status1], methodspearman)如果负相关且AUC小于0.5就取负号。最后在输出PDF前用dev.off()收尾否则ROC.pdf文件中看不到图形如果你在RStudio里反复运行脚本建议每次绘制前先graphics.off()清空画布避免叠加残留。本文还有配套的精品资源点击获取
