校准曲线这东西凡是做过临床预测模型的人应该都不陌生。不管是用Cox回归做生存预测还是用Logistic回归做二分类诊断模型建完之后审稿人几乎必问一句你们的模型校准度怎么样尤其如果你准备按TRIPOD声明来报告模型校准曲线和区分度指标就是绕不开的两块内容。我在之前的系列文章里讲过Cox模型的构建、验证、列线图绘制这篇单独把校准曲线拎出来专门聊清楚它在Cox模型里到底怎么画、怎么看、怎么解读顺便把能直接跑的R语言代码全部贴出来。这篇文章适合正在做临床预测模型、被校准曲线折腾过的医生和科研人员也适合R语言刚入门但想搞懂校准概念的朋友。1. 校准曲线到底在做什么1.1 模型评估的两把尺子区分度与校准度很多人一提预测模型性能第一反应就是C-index或者AUC。C-index确实重要它回答的问题是在随机抽取的一对患者中模型能不能把更早发生事件的那位排在风险更高的位置。换句话说它衡量的是模型的排序能力叫区分度discrimination。但区分度高不代表模型的预测概率就对。举个例子某个模型对所有患者预测的生存率都系统性偏高10个百分点它依然可能拥有很高的C-index因为患者之间的相对排序没有变高的还是高低的还是低。可如果你用这个模型去告诉一位患者“你两年生存率有70%”而真实情况可能只有60%这就不是排序的问题了这是预测绝对值准不准的问题。校准度calibration衡量的正是这个模型预测出的概率和实际观察到的概率到底有多大差距。校准曲线就是把这种差距画出来给人看。理想状态下模型预测生存率是70%那么在实际随访数据中同样一批预测值为70%的患者两年后应该大约有70%的人还活着。所以建立预测模型区分度和校准度是两条腿缺一条都站不稳。C-index看的是“谁能排前面”校准曲线看的是“预测的概率靠不靠谱”。1.2 Cox模型校准曲线的独特之处同样是校准曲线在Logistic回归和Cox回归里画的逻辑完全不一样这是很多人第一次上手就懵掉的地方。Logistic回归预测的是一个固定时间点的结局概率比如“患者是否患病”。你只要把预测概率和实际患病比例放在一起比较就行操作简单直观。可Cox模型输出的是生存函数也就是随时间变化的生存概率。同一名患者1年生存率是一个值3年生存率又是另一个值。因此Cox模型的校准曲线必须指定一个时间点比如“2年生存率校准曲线”然后在这个时间点上比较预测生存率和观测生存率。这里还有个隐藏的坑生存数据里有删失censoring。随访结束时还没发生事件的患者不代表他永远不会发生事件只是还没观察到而已。所以“观测生存率”不能简单地用“某组内还活着的人数除以总人数”来计算必须用Kaplan-Meier估计或者伪观测值、逆删失加权这类方法来处理删失。这也是为什么Cox校准曲线比Logistic校准曲线难画、也更容易画错的原因。搞懂了这两点再看下面的代码和分析思路就会清晰很多。2. 画校准曲线的两套主流方案选型2.1 实战中最常用的rms路线建模校准一条龙在R语言里面画Cox校准曲线医学统计领域用得最多的方案是rms包。Frank Harrell开发的这个包几乎是临床预测模型的标准工具从数据描述、模型拟合、列线图绘制到模型验证全部在同一个生态里完成。rms方案的核心思路是先用cph()函数拟合Cox模型然后直接调用calibrate()函数对模型做校准。calibrate()内部采用Bootstrap重抽样得到三条曲线理想曲线Ideal、表观曲线Apparent和偏差校正曲线Bias-corrected。这三条曲线的含义后面会详细讲简单说偏差校正曲线才是你论文里最该展示的结果因为它校正了模型自身过拟合带来的乐观偏差。用rms方案有几个明显的好处。第一代码量少几行就能出图适合大多数人快速出结果。第二和rms生态下的nomogram()列线图、scores()评分表衔接流畅你前一步刚建好列线图下一步直接校准数据格式不用来回转换。第三calibrate()的输出不仅有图还有量化的校准统计指标方便写进论文结果部分。不过这个方案也有局限。它的参数设置比较讲究比如必须提前设定好time.inc模型拟合时必须加上x TRUE, y TRUE很多新手在这些细节上栽跟头。另外它默认的图形样式偏基础要达到发表级别的美观程度通常还需要用ggplot2重构。2.2 更灵活的riskRegression路线多模型、多时间点一条龙如果你的需求不是简单画一条曲线而是想在一个图中比较多个模型的校准情况或者同时查看多个时间点的校准曲线那我更推荐riskRegression包。riskRegression的核心函数是Score()功能很强大。你传入一个已经拟合好的Cox模型可以是survival包的coxph()不要求必须是rms的cph()它就能计算出Brier分数、时间依赖性AUC、校准曲线等多种评估指标。配合plotCalibration()函数可以快速画出校准图。riskRegression还有一个优势它对时间点非常友好。比如你想同时看1年、3年、5年的校准曲线一条代码就能在同一个坐标系里叠加显示这在rms方案里操作起来要麻烦得多。我在做多时间点生存模型验证时基本都会用riskRegression来出一版对比图。缺点是riskRegression包的API风格和rms差别较大它的Hist()对象和Surv()对象经常让初学者头晕而且它的输出内容相对更“统计风味”一些需要花点时间理解各项指标的含义。2.3 方案对比小结维度rms方案riskRegression方案支持模型类型主要支持cph()对象支持coxph()等多数生存模型建模生态统一性高可与列线图、验证集等直接衔接一般需要额外建模多时间点校准需要分别调用或循环处理原生支持一次出多个时间点多模型对比不直观原生支持一张图比较多个模型图形美观度基础图形较朴素相对美观但仍有美化空间上手难度参数细节多容易踩坑需要适应Hist()等API风格到底选哪套我的建议很直接如果你只是为当前这个Cox模型补一张校准曲线且已经用了rms建模那就用rms方案如果你要做多模型对比、多时间点验证或者用的是survival包建模我推荐riskRegression。两套代码我下面都会贴出来你可以都跑一遍看哪个更顺手。3. 全套实操从数据到发表级校准曲线3.1 准备工作与核心陷阱在写正式代码前先把演示环境准备好。下面的代码全部用R语言实现测试数据是survival包自带的lung数据集。这是一份晚期肺癌随机临床试验数据包含228名患者的生存时间time单位是天、生存状态status1为删失2为死亡以及年龄、性别、ECOG评分等变量。用它演示Cox校准曲线非常合适。如果你电脑上还没装R环境建议先去R语言官网下载最新版R再装一个RStudio作为IDE。安装好之后执行下面的代码安装本次需要用到的包install.packages(rms) install.packages(survival) install.packages(riskRegression) install.packages(ggplot2)数据准备和预处理的时候我先提醒两个高频错误几乎所有刚接触这份数据的人都会踩。第一个坑是lung数据里status变量的编码。在这个数据集里status 1表示删失还没观察到死亡status 2表示事件发生死亡。而在R的Surv()函数里内部约定通常是1表示事件发生0表示删失。因此你需要把原来的status变量做一次映射转换否则后面的分析会得出完全错误的结果。第二个坑是缺失值。lung数据里的ph.ecog有几个缺失值如果不处理cph()函数在拟合时会默认使用na.action机制删除缺失行但后续datadist()和calibrate()可能因为变量不同步而报错。稳妥起见演示代码里直接用na.omit()把缺失行删掉。library(survival) library(rms) library(riskRegression) library(ggplot2) # 加载数据并剔除缺失值 data(lung) lung2 - na.omit(lung) # 转换生存状态编码原编码2死亡这里转换为1事件 lung2$status2 - ifelse(lung2$status 2, 1, 0) head(lung2)转换完之后status2就可以直接放进Surv()里了。接下来正式建模。3.2 rms方案代码cph calibrate plotrms这套流程有几个必须记住的“仪式”少一个都可能出问题。第一步是用datadist()设置数据分布并把它指定为全局选项第二步是在cph()里设置x TRUE, y TRUE这样calibrate()才能拿到原始数据做重抽样第三步是设置time.inc这是校准的时间点单位必须和time变量一致。假设我们要画两年生存率校准曲线time是以天为单位那么两年就是365.25 * 2 730.5天。这个换算容易出错尤其当原始时间是“月”时千万别习惯性地写24。下面是可以直接跑的完整代码# 1. 设置datadistrms生态必需 dd - datadist(lung2) options(datadist dd) # 2. 设定校准时间点2年 time_inc - 365.25 * 2 # 3. 拟合Cox模型 fit_cph - cph(Surv(time, status2) ~ age sex ph.ecog, data lung2, x TRUE, y TRUE, time.inc time_inc, surv TRUE) # 4. 查看模型概要 print(fit_cph)代码里我加了surv TRUE这个参数让模型额外保存生存函数估计保证后续绘图时能顺利计算每个个体的预测生存概率。缺了这个参数在某些版本的rms里calibrate()会表现出莫名其妙的问题。接下来就是核心一步# 5. Bootstrap法校准200次重抽样分20组比较 set.seed(2024) # 设随机种子保证结果可复现 cal - calibrate(fit_cph, u time_inc, # 校准时间点 B 200, # 重抽样次数 pred.group 20, # 将预测概率分为20组 cmethod hare) # 采用HARE方法估计观测生存率 # 6. 绘图 plot(cal, xlim c(0, 1), ylim c(0, 1), xlab Predicted 2-Year Survival Probability, ylab Observed 2-Year Survival Probability, subtitles TRUE)跑完这段代码你会看到一张经典的校准曲线图。图中通常有三条线一条是对角虚线代表理想状态下的完美校准线一条是表观校准线显示模型在训练数据上的预测和观测一致性还有一条是偏差校正后的校准线它更真实地反映了模型在新患者群体中可能的表现。B 200这个参数是Bootstrap重抽样次数。200次是比较常规的默认选择既能保证结果稳定速度也不会太慢。如果样本量大、事件数多可以提高到500次但代价是计算时间明显增加。我自己的经验是在228例这种小样本上200次足够再往上提升对最终图形的影响肉眼几乎看不出来。pred.group 20表示把患者按预测概率分成20个小组每组内分别比较平均预测生存率和KM估计的观测生存率。分组越多曲线上的点越密看起来越细腻但如果样本量不足每组内的人数太少KM估计就会变得很不稳定。小样本数据建议用15到20组大样本可以适当增加到30组甚至更多。cmethod hare是观测生存率的估计方法。HARE是一种基于风险回归的平滑估计方法比直接使用Kaplan-Meier更稳。一般保持默认就行如果希望结果更传统也可以改成cmethod KM两者结果差距不大。3.3 读懂calibrate的输出与统计表格图画出来之后很多人只看一眼曲线形态其实calibrate()还输出了非常有价值的量化信息。直接打印cal对象可以看到图例中的统计量比如样本量n、事件数d、分组数p以及模型在第pred.group百分位点上的C统计量。用summary()还能得到更详细的表格显示在50%和90%分位点上的预测生存率与观测生存率# 查看校准汇总结果 summary(cal)输出的结果大致长这样分组百分位预测生存率Predicted观测生存率Observed50% 分位点0.3210.33290% 分位点0.6010.610这个表格很实用。比如在50%分位点模型预测的中位生存概率是0.321而实际观察到的同等患者两年生存率是0.332两者相差0.011说明在校准中位水平上模型的预测偏差很小。90%分位点类似。我写论文时通常会把这个表整理成文字放进结果部分比如“模型在50%和90%分位点上的预测生存率与实际观测生存率相差均在0.02以内校准良好”。判断校准好坏不要只看整条曲线贴不贴对角线这两个分位点的数值差其实更直观。差异在0.05以内一般算可接受超过0.1就要警惕模型在该风险区间的预测能力可能存在问题。3.4 用ggplot把rms校准曲线重绘成论文风格rms自带的plot()画出的图准确是准确但样式比较老旧。很多期刊对图片清晰度、配色、字体都有要求直接在原图基础上用还是会显得敷衍。我通常的做法是提取calibrate()对象内部的数据转成数据框后用ggplot2重新绘制这样想怎么调就怎么调。calibrate()返回的对象里apparent存储的是表观校准数据calibrated存储的是偏差校正数据两者都是两列矩阵第一列是预测生存率第二列是观测生存率。提取出来就能画# 提取校准数据 cal_app - as.data.frame(cal$apparent) cal_corrected - as.data.frame(cal$calibrated) # 添加数据来源标记方便合并绘图 cal_app$type - Apparent cal_corrected$type - Bias-corrected cal_all - rbind(cal_app, cal_corrected) colnames(cal_all)[1:2] - c(Predicted, Observed) # ggplot2重绘 ggplot(cal_all, aes(x Predicted, y Observed, color type)) geom_abline(intercept 0, slope 1, linetype dashed, color black) geom_line(aes(group type), linewidth 0.8) geom_point(size 2) coord_equal() scale_color_manual(values c(Apparent gray50, Bias-corrected firebrick)) labs(x Predicted 2-Year Survival Probability, y Observed 2-Year Survival Probability, color ) theme_bw(base_size 12) theme(legend.position right)这种图放在论文里比默认的base R图形美观很多。需要特别说明的是cal对象里是否总是同时包含apparent和calibrated跨版本可能略有差异但在我测试过的rms 6.x版本里都是这样的结构。如果你运行时报列名找不到可以用str(cal)查看一下对象内部结构再提取。另外提醒一句保存图片时务必设置合适的尺寸和分辨率。如果走期刊投稿一般要求300 dpi以上可以用ggsave(..., dpi 300, width 5, height 4)保存字体和线条粗细也要根据期刊要求调整。3.5 riskRegression方案代码coxph Score plotCalibration如果你的模型是用survival::coxph()拟合的不一定非要用rms重写一遍直接用riskRegression就能完成校准。而且这套方案在绘制多个时间点的校准曲线时优势明显。先拟合一个普通的Cox模型# 用survival包拟合Cox模型 fit_cox - coxph(Surv(time, status2) ~ age sex ph.ecog, data lung2, x TRUE, y TRUE)然后调用Score()# 计算校准指标 score_obj - Score( list(Cox model fit_cox), formula Hist(time, status2) ~ 1, data lung2, times 730.5, # 2年时间点 metrics brier, # 同时计算Brier分数 plots calibration ) # 绘制校准曲线 plotCalibration(score_obj, times 730.5, type survival)这里有两个容易忽略的细节。第一Score()的formula参数必须用Hist(time, status2) ~ 1而不是Surv(time, status2) ~ 1这是riskRegression包的固定风格。第二times参数指定想要校准的时间点单位依然要和原始数据一致用天就填730.5。plotCalibration()画出来的图是分组散点加理想对角线每个点代表预测生存率相近的一组患者。如果同时传入多个模型还可以在图中叠加显示多条校准曲线用于模型间的横向比较# 多模型校准曲线示例 fit_cox2 - coxph(Surv(time, status2) ~ age sex, data lung2, x TRUE, y TRUE) score_obj2 - Score( list(full model fit_cox, reduced model fit_cox2), formula Hist(time, status2) ~ 1, data lung2, times 730.5, metrics brier, plots calibration ) plotCalibration(score_obj2, times 730.5, type survival)这样的图特别适合在模型比较章节展示完整模型和简化模型的校准线都靠近理想线但完整模型的曲线更贴近对角线说明增加变量确实改善了校准度。审稿人看到这种图往往会觉得你的验证工作做得很扎实。如果你想同时看1年、2年、3年三个时间点的校准曲线times c(365.25, 730.5, 1095.75)即可plotCalibration()会按时间点分别出图非常方便。4. 校准曲线常见问题与排查实录4.1 时间点选择的经验与单位转换校准曲线的“校准时间点”选择是我见过出错频率最高的地方。绝大多数人的困惑是我应该选择哪个时间点我的经验是优先选择临床上有决策意义的时间点而不是机械地使用中位随访时间。比如做肺癌术后预测模型医生和患者最关心的是3年或5年生存率那校准时间点就选3年或5年别因为随访时间只够2年就强行画5年节点会因为没有足够的风险集样本而导致曲线末端剧烈波动。时间单位这个坑也一样常见。如果你的time变量是以天为单位时间点就只能写天数如果是月就写月份数如果是年就写年数。经常有人time.inc默认没设或者设成1结果模型拟合和校准时间对不上画出来的曲线完全不知所云。交叉检查的方法很简单summary(fit_cph)会输出预先设定的time.inc值先确认这一项是不是你想校准的时间节点再往下走。如果你的数据是以年为单位比如time记录为2.3、5.1这种那time.inc 2就代表2年。我建议无论原始单位是什么最好在数据准备阶段把time统一成“年”或者“天”并在代码注释中标明免得隔一个月回来看自己都忘了。4.2 内部校准与外部校准的区别很多人的校准曲线是在建模数据本身上画的这在严格意义上是“内部校准”。calibrate()里用Bootstrap重抽样的目的就是通过反复在重抽样样本上拟合模型、再回到原始数据验证来估算模型因为过拟合产生的乐观度最后把这种乐观度从表观校准里减掉得到偏差校正校准曲线。所以rms画出来的Bias-corrected线比直接在原始数据上比较预测和观测的Apparent线更有参考价值。如果你的数据量足够大更好的做法是预留独立验证集或者使用外部队列在验证集上直接绘制校准曲线。这种情况下不需要Bootstrap校正直接把模型预测值应用到新数据按预测概率分组比较各组预测生存率与KM估计的观测生存率即可。代码实现也不难把训练集上拟合好的模型拿到验证集上用predict()算出验证集每位患者的预测生存概率然后按分位数分组每组内用survfit()估计观测生存率最后画散点。这种方式得到的是真正的外部校准证据临床说服力远高于内部校准。4.3 常见报错速查表下面把我这些年跑校准曲线时遇到的典型报错整理成表每一条都是实测踩过的坑。报错信息原因解决办法fitmust havexTRUEandyTRUE拟合cph()或coxph()时没保存原始数据在模型函数中加上x TRUE, y TRUEThe time.inc variable was not defined调用calibrate()时没有指定时间点且cph()也没预设time.inc在cph()里加上time.inc或在calibrate()中指定u参数object dd not found忘记设置datadist执行dd - datadist(data); options(datadist dd)cannot use~for formula或Hist相关报错在riskRegression里误用Surv()把formula参数改为Hist(time, status) ~ 1There were missing values in...数据里有缺失值建模前用na.omit()或在模型中显式处理缺失推荐先做多重插补再建模校准曲线末端点突然翘起或断裂时间点太靠后该时间点还在风险集内的样本太少检查summary(survfit(...), times ...)确认时间点处仍有足够风险集或改选更早的时间点4.4 几个容易忽视的小细节细节一随机种子的设置。calibrate()依赖Bootstrap重抽样不设种子的话每次运行结果会在小数点后有微小波动。虽然不影响整体判断但投稿时如果审稿人要求复现不同次运行结果不一致会给人很不严谨的印象。所以在calibrate()之前一定加set.seed()。细节二样本量与分组数的匹配。pred.group不是越大越好。我见过有人把预测概率分成50组画出来的曲线密密麻麻全是锯齿反而看不出趋势。样本量在200例左右时20组是安全的选择样本量上千的时候可以考虑30到40组。如果你的模型目的只是看整体趋势分组少一点反而更稳健。细节三一定要结合Brier分数。校准曲线是“看图说话”Brier分数则把预测误差变成了一个数值结合使用更利于论文量化报告。riskRegression的Score()函数里metrics brier可以直接算出这个值Brier分数越低说明预测越准确生存模型中它同时考虑到了校准度和区分度。细节四校准曲线和列线图的关系。如果你之前已经画过列线图那么这两者是配套的。列线图展示的是“怎么利用模型去预测”校准曲线回答的是“这个预测到底有多靠谱”。做临床预测模型报告时先列线图再校准曲线再加C-index这一套组合拳是标准配置。关于校准曲线最后再分享一点个人体会。我在做模型验证时发现校准曲线的形态其实比C-index更容易暴露模型的真实问题。有时候模型C-index看着还行但校准曲线明显偏离对角线说明模型在某些风险区间存在系统性偏差。这类问题光靠调阈值、改截断值掩盖不了得从模型结构、变量筛选或者数据质量上找原因。画校准曲线不是流程上的走过场它是真正帮你审视模型质量的一盏灯。希望这篇的代码和避坑经验能让你少走一些我当时走过的弯路。