如果你已经掌握了生存分析的入门流程——先画 KM 曲线再做 log-rank 检验然后跑一个 Cox 回归最后从输出里读出 HR 和 P 值——那么再往前走一步时大概率会遇到下面三个尴尬问题。第一个问题研究终点是“疾病进展”但有一部分患者在进展前先去世了。直接用 KM 估计“进展概率”曲线会偏高因为死亡被当成删失处理了。第二个问题Cox 模型的结果很漂亮但用cox.zph()做比例风险假设检验时 P 值小于 0.05模型可能站不住脚。第三个问题导师或业务方问的不是“HR 是多少”而是“这位患者 1 年生存率是多少”你给出的个体化概率如何计算又如何证明它靠谱。这三个问题正好对应本文要讲的“生存分析进阶三板斧”第一斧竞争风险模型第二斧Cox 回归的假设诊断与修正第三斧从 HR 到个体概率的预测与模型验证。读完这篇文章你可以用 R 把一条完整的进阶分析流程跑通并且理解每一步背后的统计逻辑而不是只停留在“能用函数、能出图”的阶段。1. 为什么说 KM 曲线和 Cox 回归只是入口很多初学者把 KM 曲线和 Cox 回归当成生存分析的全部这可以理解因为它们确实覆盖了从“可视化”到“影响因素分析”的核心环节。但真实研究里的数据往往比教程例子复杂得多有删失、有竞争事件、有随时间变化的效应还要面对“这个模型到底准不准”的追问。KM 和 Cox 解决的是基础问题进阶三板斧解决的才是实战问题。1.1 生存分析回答的是“时间到事件”问题生存分析的标准三个要素是生存时间、终点事件、删失状态。以临床研究为例生存时间通常是从入组到事件发生的时间终点事件可以是死亡、复发、疾病进展删失表示观察结束时事件尚未发生或者患者失访、退出研究。很多新手只关注“事件是否发生”却忽略了“时间长度”才是生存分析的因变量这是一个很容易犯的错。从工程和数据视角看生存分析也不只属于医学。用户流失预测、设备故障时间、风控中的违约时间、营销活动中的响应时间本质上都是“时间到事件”数据。理解这一点后你会发现本文提到的竞争风险、PH 假设、模型验证在业务分析中同样成立。1.2 KM 曲线与 Cox 回归的边界在哪里KM 曲线的价值是描述性、探索性的它能给出不同组别的生存概率曲线并用 log-rank 检验比较组间差异。但它不能同时纳入多个变量也不能给出风险大小。Cox 回归的价值是解释性、推断性的它通过比例风险假设把多个协变量与生存风险联系起来输出的是 HR 和置信区间。但是这两个工具的边界也很明显问题KM 曲线Cox 回归能否处理竞争风险不能会把竞争事件当删失不能直接回答累积发生率是否依赖比例风险假设不涉及依赖且需要检验能否给出个体化概率只能给组别概率可以但需要额外转换和验证能否回答“预测准不准”不能不能直接回答所以进阶三板斧的本质是在 KM 和 Cox 之上补上三个能力处理多终点、校验模型假设、做个体预测与验证。它们不是替代 KM 和 Cox而是让整个分析更接近真实问题。2. 第一板斧竞争风险模型竞争风险Competing Risks是生存分析里最容易被忽略、又最容易造成结论错误的问题。如果你的目标事件不是“最终事件”那么你很可能已经遇到了竞争风险只是还没有意识到。2.1 什么是竞争风险竞争风险指的是在目标事件发生之前个体发生了另一事件导致目标事件无法再被观察到。比如研究“疾病进展”患者死亡后就不会再发生进展研究“首次入 ICU”患者转院或死亡后就不再处于风险状态研究“用户再次购买”用户注销账号后就无法再购买。关键点在于竞争事件不是随机删失而是改变了事件结局的“竞争终点”。如果忽略它等于假设竞争事件永远不会发生这显然不符合真实世界。2.2 为什么 KM 和普通 Cox 会失真KM 估计处理删失时默认删失个体在未来仍然有机会发生目标事件。但在竞争风险场景下已经发生竞争事件的个体理论上永远不会再发生目标事件。把竞争事件当作删失会把“不可能再发生事件的人”继续留在风险集里从而高估目标事件的累积发生率。普通 Cox 回归也有类似问题。当我们写coxph(Surv(time, status 1) ~ ...)竞争事件会被当作删失处理。模型估计的是 cause-specific hazard也就是“在尚未发生任何事件且仍然存活的前提下单位时间内发生目标事件的瞬时概率”。这个解释在病因学研究中很有意义但它不能直接回答“患者最终有多大可能发生进展”这类预后问题。回答这个问题需要用累积发生率函数 CIF。2.3 Fine-Gray 模型与 cause-specific 模型怎么选竞争风险建模有两种常见思路一是 cause-specific hazard 模型基于风险集做 Cox 回归关注的是特定事件的原因别风险。它适合回答“某个因素是否影响某条事件路径的瞬时发生率”。二是 Fine-Gray 子分布风险模型直接对累积发生率函数建模关注的是目标事件最终发生的累积概率。它更适合回答“某个因素是否增加或降低目标事件的最终风险”常用于预测和预后研究。维度Cause-specific hazardFine-Gray 子分布风险建模对象特定事件的瞬时风险目标事件的累积发生率风险集仍未发生任何事件的人未发生目标事件的人包括已发生竞争事件的人临床问题这个因素是否影响该事件发生速度这个因素是否影响最终发生概率常用方法Cox 分层 竞争事件按删失cmprsk 包中的 crr()适用场景病因学研究预后预测、临床决策这二者的结果可能方向不同所以论文或报告中一定要写清楚用的是哪一种模型、回答的是哪一个问题。2.4 R 实现用 mgus2 数据集跑竞争风险模型以survival包自带的mgus2数据为例。这是“意义未明的单克隆丙种球蛋白病”随访数据fstat的取值含义是0 表示删失1 表示疾病进展2 表示死亡。如果把“疾病进展”作为目标事件死亡就是典型的竞争事件。先准备数据library(survival) library(cmprsk) data(mgus2, package survival) # 筛选完整变量避免 crr() 因缺失值报错 mgus_clean - mgus2[complete.cases(mgus2[, c(age, sex, hgb, creat)]), ] table(mgus_clean$fstat)这里没有对fstat做额外编码因为crr()要求事件状态是整数并且通过failcode和cencode来区分目标和删失代码。先画累积发生率函数 CIF并做 Gray 检验cif_fit - cuminc(ftime mgus_clean$futime, fstatus mgus_clean$fstat, group mgus_clean$sex) print(cif_fit) plot(cif_fit)print(cif_fit)会输出分组的 Gray 检验结果重点关注pv列。plot(cif_fit)画出不同事件、不同组别的累积发生率曲线注意它和 KM 曲线的区别CIF 曲线会同时考虑竞争事件的影响终点不会超过竞争事件的累积风险。接着建立 Fine-Gray 回归模型cov_mat - model.matrix(~ age sex hgb creat, data mgus_clean)[, -1] fgr_fit - crr(ftime mgus_clean$futime, fstatus mgus_clean$fstat, cov1 cov_mat, failcode 1, cencode 0) summary(fgr_fit)failcode 1表示把疾病进展作为目标事件cencode 0表示把 0 作为删失代码那么 fstatus 中值为 2 的死亡事件就自动作为竞争事件处理。cov_mat是协变量矩阵需要手动构造这也是crr()与coxph()在接口上的最大区别。2.5 结果怎么看summary(fgr_fit)的输出包含系数、exp(coef)、标准误、z 值和 P 值。这里的exp(coef)是子分布风险比解释时要避免写成普通 Cox 的 HR。更准确的说法是在其他变量不变时该变量的值每增加一个单位目标事件累积发生率对应的子分布风险变化多少倍。一个典型的错误是把竞争风险模型的结果直接等同于普通 Cox 结果。它们估算的数学对象不同解读语境也不同。如果研究目的是描述病因可以同时报告 cause-specific 模型和 Fine-Gray 模型如果研究目的是预测Fine-Gray 和 CIF 通常更直观。3. 第二板斧Cox 回归的 PH 假设检验与修正如果说竞争风险解决的是“终点事件定义”问题那么比例风险假设解决的是“模型是否合理”的问题。Cox 回归的魅力在于它不要求指定基线风险函数但这个灵活性是有代价的代价就是 PH 假设。3.1 比例风险假设想表达什么Cox 模型的核心表达式是h(t|X) h0(t) * exp(Xβ)其中 h0(t) 是基线风险函数Xβ 是协变量的线性组合。这个公式的意思是所有个体的风险函数形状相同协变量只是按比例放大或缩小风险。换句话说某个变量的 HR 在整个随访期间保持不变。但真实数据经常不满足这个假设。比如年龄对早期死亡的影响可能很大对晚期死亡的影响可能减小新药的短期效果明显但长期效果被耐药性抵消。这时一个单一的 HR 无法准确描述协变量的作用。3.2 用 cox.zph 做 PH 检验R 中常用cox.zph()基于 Schoenfeld 残差进行检验data(lung, package survival) # 将 status 从 2/1 编码转换为 1/0便于后续建模 lung2 - lung lung2$status - ifelse(lung2$status 2, 1, 0) lung2$sex - factor(lung2$sex, levels 1:2, labels c(male, female)) fit_cox - coxph(Surv(time, status) ~ age sex ph.ecog, data lung2) summary(fit_cox) zph_fit - cox.zph(fit_cox) print(zph_fit) plot(zph_fit)print(zph_fit)会对每个协变量输出卡方值和 P 值。如果某个变量或整体模型的 P 值小于 0.05就说明该变量可能存在对时间的依赖PH 假设存疑。plot(zph_fit)画出残差随时间变化的曲线如果曲线有明显的趋势而不是围绕水平线波动说明效应确实随时间变化。这里要提醒一点PH 检验不是“过不过”的问题而是“偏离程度是否影响解释”的问题。样本量大时微小的偏离也可能得到显著 P 值样本量小时明显的偏离也可能检验不出来。所以不要只盯 P 值还要结合残差图和专业判断。3.3 不满足 PH 怎么办分层和时变系数处理 PH 假设不满足有两条常见路径。第一种是分层。比如性别不满足 PH 假设将sex作为分层变量fit_strata - coxph(Surv(time, status) ~ age ph.ecog strata(sex), data lung2) summary(fit_strata)strata(sex)允许男性和女性拥有不同的基线风险函数但不估计sex的主效应。这样既保证了 PH 假设的基本成立又承认了性别对生存曲线的整体影响。代价是你不能再直接解读“性别的 HR”。第二种是时变系数。如果连续变量如年龄随时间变化可以用tt()在模型中加入交互项fit_tt - coxph(Surv(time, status) ~ age tt(age) sex ph.ecog, data lung2, tt function(x, t, ...) x * t) summary(fit_tt)这里的tt()告诉coxph()对 age 进行时间变换函数function(x, t, ...) x * t表示把 age 的效应建模为随 t 线性变化。实际项目中t的变换形式可以是log(t)、sqrt(t)甚至更灵活的样条函数。选择哪种形式需要结合专业背景和模型拟合指标判断而不是机械地套用。3.4 模型比较的参考指标加入时变项之后模型不一定更好还需要比较AIC(fit_cox, fit_tt)AIC 越小说明模型在拟合优度和复杂度之间取得了更好的平衡。如果加入tt(age)后 AIC 明显下降说明时变效应确实存在如果 AIC 基本不变甚至上升则说明复杂化没有必要。不要只用 P 值决定模型形态这是进阶分析和基础分析的一个重要区别。4. 第三板斧从 HR 到个体预测与模型验证很多分析做到 Cox 回归就结束了但 Cox 回归的输出是变量效应不是个体概率。实际业务或临床中需要回答的问题是一个具体的人在某个时间点发生事件的概率是多少。4.1 HR 与个体预测之间的距离HR 是协变量每增加一个单位时风险的比例变化它描述的是“平均效应”。但两个病人即使拥有相同的 HR 相关变量他们的绝对风险也可能不同因为绝对风险还受基线风险函数和随访时间影响。所以从 Cox 模型到个体预测需要把回归系数和累积基线风险结合起来计算出个体在给定时间点的生存概率。这一步在 R 中通常由rms包完成。4.2 用 rms 构建列线图列线图Nomogram是把 Cox 模型预测可视化的一种方式。它将回归系数转换为可读的分数用简单的刻度直接读出个体在某时间点的生存概率。library(rms) # 为 rms 设置数据分布对象 ddist - datadist(lung2) options(datadist ddist) # 用 cph 重新拟合 Cox 模型xTRUE、yTRUE、survTRUE 必不可少 fit_cph - cph(Surv(time, status) ~ age sex ph.ecog, data lung2, x TRUE, y TRUE, surv TRUE) # 构造生存概率函数 surv_prob - Survival(fit_cph) nom - nomogram(fit_cph, fun list( function(x) surv_prob(365, x), function(x) surv_prob(730, x) ), funlabel c(1-year survival probability, 2-year survival probability)) plot(nom)使用rms时最容易踩的坑有两个。第一必须设置options(datadist ddist)否则nomogram()会报错。第二cph()必须带上x TRUE, y TRUE, surv TRUE否则后面的Survival()和calibrate()无法工作。4.3 C 指数与校准曲线预测模型不能只看是否显著还要看区分度和校准度。区分度衡量模型能否把高风险和低风险个体分开常用 C 指数实现。concordance(fit_cox)C 指数取值范围是 0.5 到 10.5 表示完全随机1 表示完全一致。需要注意的是C 指数对样本组成敏感不同队列之间的 C 指数不能简单直接比较。校准度衡量“预测概率”和“实际概率”是否一致。用calibrate()可以观察set.seed(2024) cal_res - calibrate(fit_cph, u 365, B 200) plot(cal_res)u 365表示评估 1 年生存概率的校准B 200表示自举重抽样次数。理想情况下校准曲线应该贴近对角线。如果曲线偏离严重说明预测概率存在系统性的高估或低估。4.4 决策曲线分析C 指数和校准曲线回答的是“模型准不准”决策曲线分析DCA回答的是“模型有没有用”。DCA 把一个模型在不同阈值下的净获益与“所有人都不治疗”“所有人都治疗”两个极端策略进行比较适合评估临床或业务决策场景中的增量价值。R 中可以通过dcurves或rmda包实现生存数据的 DCA。这两个包在不同版本中的函数参数有所调整使用前建议先查看包文档。DCA 的结果一般不单独看曲线高不高而是看新模型相比旧模型是否在某些阈值区间内带来净获益。4.5 验证策略个体预测模型如果没有验证很容易过拟合。常见验证策略分为内部验证和外部验证。内部验证是在当前数据集内评估常用的方法包括数据分割和自举法。自举法优于简单分割因为它能在不损失样本量的情况下给出更稳定的乐观度估计。前面calibrate()中的B 200就是自举法的一种应用。外部验证是在另一个中心、另一个时间段或另一批人群上验证模型这是预测模型证据等级中最可靠的验证方式。如果条件允许任何要进入实际决策的模型都应该最少完成一次外部验证。5. 把三板斧接成一条完整分析流水线在实际项目中三板斧并不是彼此独立的三个步骤而是一条完整的决策链。下面是一个参考工作流明确研究问题目标事件是什么竞争事件是什么时间原点从哪里开始数据预处理检查删失编码、缺失值、变量类型。探索性分析先画 KM 曲线但若存在竞争风险必须补充 CIF 和 Gray 检验。主分析用 Cox 回归并做 PH 假设检验若不满足使用分层或时变系数。预测建模用cph()拟合模型构建列线图。模型评估计算 C 指数、绘制校准曲线、做决策曲线分析。验证与报告记录随机种子、样本量、事件数、验证方式并按规范报告结果。这套流水线中每一步的输出都会影响下一步而不是机械地“跑一个模型然后截图”。很多数据分析报告之所以不可信不是因为代码报错而是因为事件定义不清、竞争事件被忽略、PH 假设没有检查、预测模型没有验证。如果你是 Python 技术栈也想在业务项目中完成类似分析可以这样对应lifelines支持 KM、Cox、PH 假设检验核心类是 CoxPH