做临床随访或者真实世界数据分析的人大概率都经历过这样的流程拿到一张带着time列和status列的表格先画 KM 曲线看 log-rank 的 p 值然后跑 Cox 回归报告 HR 和 95% 置信区间最后写进论文。这套组合拳确实是生存分析的入门标配熟练到几乎不用思考。但数据不会永远这么听话。随访时间一长患者会死于各种原因你研究的是肿瘤特异性生存可很多患者最后死于心脑血管事件你研究的是手术 A 和手术 B 的长期疗效可术后早期 A 组手术风险高、远期复发风险反而低HR 并不是一个常数你辛辛苦苦跑出 HR0.72审稿人却反问一句“这个相对风险在临床上到底意味着患者多活了几个月”这三个问题正好对应进阶生存分析最常见的三板斧竞争风险模型、时变效应处理、RMST 限制平均生存时间。它们不是炫技而是基础 KM 曲线和 Cox 回归覆盖不了的真实场景。这篇文章会用可复现代码把这三板斧讲清楚并给出模型选择、结果报告和审稿应对的思路。1. 基础KM曲线与Cox回归的边界在哪里1.1 基础方法适合什么场景先别急着否定 KM 曲线和 Cox 回归。它们仍然是生存分析的基石在很多场景下完全够用随访期相对短研究期间事件类型单一患者结局只有“发生”和“删失”两种状态不同治疗组的风险比在整个随访期内大致恒定研究问题停留在“有没有差异”“有没有关联”层面。比如一个随访 12 个月的术后并发症研究终点是明确的“并发症发生”患者即使死亡也会被当作删失处理医药干预对事件的影响在短时间内不会剧烈变化。这种情况下KM 曲线 log-rank 检验 Cox 回归就是标准答案没有必要堆砌复杂模型。1.2 基础方法撑不住的三种情况基础流程的核心假设是“删失机制与结局无关且终点事件唯一”。一旦以下三种情况出现结果就可能出现系统性偏差情况一存在竞争风险。随访期内患者可能死于与目标结局无关的原因。例如研究“肿瘤特异性死亡”但大量患者死于心血管疾病。普通 KM 曲线会把心血管死亡当作删失相当于假设这些患者“仍然有机会在未来发生目标事件”但现实中他们已经不可能再发生肿瘤死亡了。这会让目标事件的累积发生率被高估。情况二HR 不恒定。Cox 回归输出的是一个平均效应。如果治疗在早期有效、晚期无效或者某种治疗前期风险高、后期收益大单个 HR 就会掩盖这种时间动态变化。很多人一看到cox.zph()的 p 值小于 0.05 就慌了不知下一步该怎么处理。情况三只有相对效应缺少绝对效应。HR0.72 在统计上是显著的但临床上医生更想问的是“用新方案患者 3 年内平均能多活几个月”这一句话 HR 回答不了需要 RMST。一句话总结基础与进阶方法的分工基础方法回答“有没有关系”进阶方法回答“这个关系是否恒定、是否被其他事件干扰、临床上到底有多大”。2. 生存分析核心概念先花三分钟统一语言2.1 一条生存数据由什么构成进入进阶方法之前先统一概念。一条标准的生存数据有三个要素要素含义例子时间 T从起点到事件发生或删失的时间手术日期到死亡日期事件状态 status是否发生事件以及事件类型0删失1目标事件2竞争事件协变量 X可能影响生存的特征年龄、性别、治疗组、分期普通 Cox 回归要求status是 0/1 二值变量。到竞争风险模型时status会变成 0/1/2 多分类这是很多人第一次进阶时卡住的点数据结构里根本没有为“第二种事件”预留位置。2.2 生存函数、风险函数与CIF几个必须分辨的指标生存函数 S(t)到时间 t 为止仍未发生事件的概率KM 曲线估计的就是它。风险函数 h(t)在 t 时刻存活的个体中瞬时发生事件的概率。累积发生率函数 CIF(t)到时间 t 为止某个特定原因事件累积发生的概率。在只有单一终点时CIF 和 1 - KM 是等价的。但存在竞争风险时两者不再相等。普通 KM 把竞争事件当作删失等于把已经从风险集中移除的患者又“保留”在分母里导致目标事件发生率被高估。这也是为什么竞争风险场景下必须画 CIF而不是画 1 - KM。2.3 删失与竞争事件的区别删失censor和竞争事件competing event在数学上都表现为“观察不到目标事件”但含义完全不同删失患者在随访截止时仍未发生事件或者失访。理论上如果他继续被观察未来仍可能发生目标事件。竞争事件患者先发生了另一种事件使目标事件在真实世界中无法再发生。比如患者死亡后不可能再“复发”。正是因为竞争事件不是删失所以不能用普通 Cox 的删失逻辑去处理。下一章进入第一板斧。3. 第一板斧竞争风险模型3.1 什么时候必须考虑竞争风险判断标准很简单研究期间内是否有一种事件会发生并且它的发生会阻止目标事件被观察到。典型场景包括肿瘤长期随访中非肿瘤死亡与肿瘤死亡竞争心血管研究中非心血管死亡与心血管事件竞争骨科植入物研究中假体失效与患者死亡竞争老年队列研究里任何“死亡”都可能与“入院”“跌倒”等事件竞争。如果竞争事件占比很低比如不足 5%用普通方法影响不大但如果随访时间长、患者年龄大竞争事件比例往往达到 20% 甚至更高这时候不做竞争风险分析结论很容易被质疑。3.2 原因别风险与子分布风险别再混淆竞争风险的建模有两套主流框架很多人一直分不清模型估计对象风险集定义适合回答的问题原因别风险模型Cause-specific每种事件的瞬时风险仍然存活、尚未发生任何事件的人群某种干预对目标事件的病因学效应Fine-Gray 子分布风险模型事件对应的 CIF 累计发生率包含已发生竞争事件的患者预测、预后分层、临床决策通俗理解原因别模型关心“在还活着的人里这个因子的致病作用有多大”Fine-Gray 模型关心“到某个时间点这个人群里有多少比例的人发生了目标事件并且竞争事件会拉低这个比例”。因此两者结果不一致很正常。原因别 HR 显著而 Fine-Gray HR 不显著说明该因子确实增加了目标事件的生物学风险但患者更容易先死于竞争事件反过来Fine-Gray 显著而原因别不显著常见于竞争事件本身受该因子影响较大的情况。报告时最好把两种模型都给出并说明研究定位是病因探索还是预后预测。3.3 R代码CIF、原因别Cox与Fine-Gray进阶生存分析的标准实现目前以 R 生态最成熟。这里用经典cmprsk包演示。假设数据有三列time为随访时间status取 0删失、1目标事件、2竞争事件以及协变量age、sex、treat。# 安装并加载包 install.packages(cmprsk) library(cmprsk) # 读取数据确认status编码 dat - read.csv(survival_data.csv) table(dat$status) # 0删失, 1目标事件, 2竞争事件 # 第一步绘制两组的目标事件CIF cif - cuminc( ftime dat$time, fstatus dat$status, group dat$treat ) # 绘图手动指定曲线图例 plot(cif, curvlab c(A组-删失, A组-目标事件, A组-竞争事件, B组-删失, B组-目标事件, B组-竞争事件))注意cuminc()默认把status0当作删失编码也就是参数cencode0。如果你项目里删失编码不是 0需要显式指定。接下来分别拟合两种模型# 第二步原因别Cox模型对目标事件建模竞争事件视为删失 cs_model - coxph( Surv(time, as.numeric(status 1)) ~ age sex treat, data dat ) summary(cs_model) # 第三步Fine-Gray子分布风险模型 fg_model - crr( ftime dat$time, fstatus dat$status, cov1 model.matrix(~ age sex treat, data dat)[, -1] ) summary(fg_model)crr()的cov1必须是矩阵形式所以用model.matrix()转换。默认failcode1表示把status1当作目标事件cencode0表示把status0当作删失与前面的数据结构一致。3.4 Fine-Gray与原因别模型的结果取舍实际项目中我更推荐这样的报告策略先用cuminc()画出各组 CIF并在关键随访时间点如 12、24、36 个月报告 CIF 估计值组间 CIF 差异用 Grays 检验cuminc()返回结果里的Tests字段就是各组事件的 p 值论文正文同时给出原因别 Cox 和 Fine-Gray 的结果避免审稿人质疑方法选择如果研究目的是预测个体预后以 Fine-Gray 为主如果是探索病因机制以原因别模型为主。需要提醒的是Fine-Gray 模型的“风险集”包含了已发生竞争事件的患者这在概率上是一个有意的构造用来直接建模 CIF。正因如此它的结果不能直接解释为“活着的患者中某个因素导致的瞬时风险”这是它被质疑最多的地方。4. 第二板斧时变效应解决HR不恒定4.1 PH假设到底意味着什么Cox 回归的前提是比例风险假设Proportional HazardsPH任意两个个体的风险函数之比在整个随访期内保持不变。也就是说治疗组的 HR 在随访第 1 个月、第 12 个月、第 36 个月都是同一个数值。这个假设在现实中经常被违背。例如免疫治疗起效慢前期可能看不出差异后期优势才显现手术干预早期风险高于保守治疗但长期预后更好某种药物的获益随时间衰减HR 从 0.6 逐渐趋近于 1。如果 PH 不成立直接报告一个 HR 会严重误导读者因为这一现象背后的真实效应是时变的。4.2 用Schoenfeld残差检验PH假设我们先用标准 Cox 模型拟合再检验 PH 假设。R 中使用cox.zph()对 Schoenfeld 残差做检验同时画出残差随时间的变化library(survival) cox_model - coxph(Surv(time, status) ~ age sex treat, data dat) summary(cox_model) # 检验PH假设 zph_result - cox.zph(cox_model) print(zph_result) # 查看某个变量残差随时间的趋势 plot(zph_result[3])cox.zph()输出中每个协变量对应一个 p 值p 值小于 0.05 说明该变量的效应随时间变化PH 假设不成立。这里的截断值可以放宽到 0.05但如果样本量很大即使很小的偏差也会显著因此还要结合残差图来判断是否存在实际意义上的趋势。4.3 时变系数模型tt()实战确认 PH 假设不满足后最直接的解决办法之一是对该变量引入时间交互项即“时变系数”模型。在survival包中通过tt()参数实现# 对treat做log(t1)交互模拟HR随时间对数衰减 cox_tt - coxph( Surv(time, status) ~ age sex treat tt(treat), data dat, tt function(x, t, ...) { x * log(t 1) } ) summary(cox_tt)这时treat在时间 t 的效应不再是常数而是HR(t) exp(β_treat β_tt * log(t 1))其中β_tt显著就说明治疗效应确实随时间变化。tt()里的变换函数可以根据数据形态替换常见的有线性t、对数log(t1)、以及分段函数。选择变换函数时最好先看 Schoenfeld 残差图是线性趋势、对数趋势还是折线趋势再决定用哪种形式不要盲目套用。另一条常见路径是分层。当某个变量本身不是关注焦点但该变量违反 PH 假设时用strata()把它分层可以得到基线风险函数不同、但其他变量仍有 HR 的模型cox_stratified - coxph( Surv(time, status) ~ age sex strata(treat), data dat ) summary(cox_stratified)分层后treat不再输出 HR但其他协变量估计更稳健。还可以用时间分段模型把随访期切成几段分别估计不同阶段的 HR。R 中使用survSplit()把数据转换成含(tstart, time]区间的长格式dat_split - survSplit( Surv(time, status) ~ ., data dat, cut c(12, 24), episode period ) cox_period - coxph( Surv(tstart, time, status) ~ treat treat:strata(period) age sex, data dat_split ) summary(cox_period)这种方法直观但切点选择要谨慎最好在研究方案里预先定义而不是根据 p 值来回试。4.4 时变协变量与常见误区上一节解决的是“协变量效应随时间变化”还有一类问题是“协变量本身随时间变化”比如患者的血压、用药剂量、疾病状态。这类问题的标准做法是把数据改成长格式用(start, stop]区间记录每个患者在不同时间段的值。Python 用户可以使用lifelines的CoxTimeVaryingFitterfrom lifelines import CoxTimeVaryingFitter # 长格式数据id_col标识患者start/stop为区间status为区间内是否事件 df_long ... ctv CoxTimeVaryingFitter() ctv.fit( df_long, id_colid, start_colstart, stop_colstop, event_colstatus, show_progressTrue ) ctv.print_summary()这里有一个必须反复强调的坑内生性时变协变量会引入严重偏倚。如果某个协变量受治疗本身影响又作为时变协变量放进模型就会形成因果倒置。比如把“患者是否发生并发症”当作时变协变量预测死亡但治疗组更容易发生某种并发症模型中就会错误地吸收治疗效应。另外时变协变量的引入还要警惕永恒时间偏倚immortal time bias。如果患者在某个事件发生之前不可能被暴露比如“术后接受了二期治疗”这个变量患者在治疗开始前的一段时间天然不可能死亡如果把这部分时间错误归入治疗组会得到偏高的保护效应。解决思路是设计阶段就用 landmark 分析或恰当的区间切分而不是建模时补救。5. 第三板斧RMST把HR翻译成人话5.1 RMST解决什么问题HR 是一个相对风险临床医生很难直接感知道“HR0.72”意味着什么。而且当 PH 假设不成立时Cox 模型输出的单个 HR 本质上是整个随访期效应的某种加权平均解释起来更加危险。RMSTRestricted Mean Survival Time限制平均生存时间的概念非常简单从 0 到某个指定时间点 tau 之间生存曲线下的面积。它表示患者在 tau 时间内平均存活的月数。举个例子两组患者随访 36 个月A 组 RMST28.5 个月B 组 RMST24.3 个月差值是 4.2 个月。这句话临床医生一听就懂A 组患者在 3 年内平均比 B 组多活了 4.2 个月。RMST 不需要 PH 假设对生存曲线尾部不稳定的问题也更稳健因此在长期随访和肿瘤临床试验中越来越受重视。5.2 tau怎么选才不被审稿人挑战RMST 的估计强烈依赖 tau 的选择因为生存曲线在 tau 之后的部分被截断掉了。推荐做法在统计分析计划SAP中提前定义 tau通常是研究中位随访时间、试验设计的最短随访时间或临床上有意义的观察窗口不要看了数据后再挑一个让结果最漂亮的 tau这属于 p-hacking报告时应做敏感性分析选择 tau±6 个月或一组合理的 tau 值看结论是否稳定。5.3 R代码survRM2的两组RMST比较R 中survRM2包是最常用的实现支持两组比较和基于伪观测法的协变量调整install.packages(survRM2) library(survRM2) # 两组RMST比较tau36个月 rmst_fit - rmst2( time dat$time, status dat$status, arm dat$treat, tau 36 ) print(rmst_fit) # 调整协变量的RMST比较 rmst_adj - rmst2( time dat$time, status dat$status, arm dat$treat, tau 36, covariates dat[, c(age, sex)] ) print(rmst_adj)这里status仍使用 0/1 编码RMST 关心的是“事件发生前平均生存时间”所以竞争风险场景下要把目标事件单独编码。如果竞争风险占比很高更严谨的做法是在竞争风险框架下比较 CIF再结合 RMST 做辅助解释。5.4 Python实现RMSTlifelines目前没有直接提供 RMST 函数但根据定义我们可以用 KM 曲线做数值积分。这条代码可以用在任何 Python 项目中import numpy as np from lifelines import KaplanMeierFitter def rmst_from_km(time, event, tau): 基于KM生存曲线的RMST计算梯形法近似面积 kmf KaplanMeierFitter() kmf.fit(time, event_observedevent) times kmf.survival_function_.index.values surv kmf.survival_function_[KM_estimate].values # 截取到tau mask times tau tt np.concatenate([[0], times[mask], [tau]]) ss np.concatenate([[1], surv[mask], [surv[mask][-1] if np.any(mask) else 1]]) return np.trapz(ss, tt) # 用法示例 rmst_a rmst_from_km(df_a[time], df_a[status], tau36) rmst_b rmst_from_km(df_b[time], df_b[status], tau36) print(fA组RMST{rmst_a:.2f}月, B组RMST{rmst_b:.2f}月) print(fRMST差值{rmst_a - rmst_b:.2f}月)手动实现的好处是逻辑透明方便集成到自己的 Python 分析流水线中缺点是没有直接给出置信区间和 p 值。正式分析时建议两两对照 R 的survRM2结果确认数值一致后再用于论文。5.5 与HR的关系和报告当 PH 假设成立时RMST 的结论通常与 HR 方向一致。但 RMST 提供的是绝对差异HR 提供的是相对风险两者是互补关系。最佳报告方式是“HR RMST差值”同时给出表格里列出中位随访时间、各组事件数、HR95% CI、RMST95% CI和 RMST 差值正文里重点解释 RMST 差值的临床含义如果 PH 不成立明确说明 HR 仅作参考主要结论以 RMST 差异或时变效应结果为准。6. 三板斧到底怎么选一个决策框架讲完三板斧很多人会问是不是以后所有分析都要全部上不是。方法论必须服务于研究问题。下面是我在项目中常用的判断流程数据场景首选方法补充方法单终点、无竞争事件、PH成立KM Cox可加 RMST 辅助随访期长、存在非目标死亡CIF 竞争风险模型RMST治疗效果随时间变化时变系数模型 / 分层 / 时间分段分阶段报告 HR需要解释临床绝对获益RMSTHR KM预后预测、个体风险分层Fine-Gray机器学习生存模型更具体的判断逻辑先看研究问题。是探究病因还是预测预后预测优先 Fine-Gray病因优先原因别模型。再画图。CIF 曲线之间是否分离、是否交叉Schoenfeld 残差图是否有明显趋势该上 RMST 就上。当随访时间确定、临床意义是关键 生存分析进阶三板斧：竞争风险、时变效应与RMST
做临床随访或者真实世界数据分析的人大概率都经历过这样的流程拿到一张带着time列和status列的表格先画 KM 曲线看 log-rank 的 p 值然后跑 Cox 回归报告 HR 和 95% 置信区间最后写进论文。这套组合拳确实是生存分析的入门标配熟练到几乎不用思考。但数据不会永远这么听话。随访时间一长患者会死于各种原因你研究的是肿瘤特异性生存可很多患者最后死于心脑血管事件你研究的是手术 A 和手术 B 的长期疗效可术后早期 A 组手术风险高、远期复发风险反而低HR 并不是一个常数你辛辛苦苦跑出 HR0.72审稿人却反问一句“这个相对风险在临床上到底意味着患者多活了几个月”这三个问题正好对应进阶生存分析最常见的三板斧竞争风险模型、时变效应处理、RMST 限制平均生存时间。它们不是炫技而是基础 KM 曲线和 Cox 回归覆盖不了的真实场景。这篇文章会用可复现代码把这三板斧讲清楚并给出模型选择、结果报告和审稿应对的思路。1. 基础KM曲线与Cox回归的边界在哪里1.1 基础方法适合什么场景先别急着否定 KM 曲线和 Cox 回归。它们仍然是生存分析的基石在很多场景下完全够用随访期相对短研究期间事件类型单一患者结局只有“发生”和“删失”两种状态不同治疗组的风险比在整个随访期内大致恒定研究问题停留在“有没有差异”“有没有关联”层面。比如一个随访 12 个月的术后并发症研究终点是明确的“并发症发生”患者即使死亡也会被当作删失处理医药干预对事件的影响在短时间内不会剧烈变化。这种情况下KM 曲线 log-rank 检验 Cox 回归就是标准答案没有必要堆砌复杂模型。1.2 基础方法撑不住的三种情况基础流程的核心假设是“删失机制与结局无关且终点事件唯一”。一旦以下三种情况出现结果就可能出现系统性偏差情况一存在竞争风险。随访期内患者可能死于与目标结局无关的原因。例如研究“肿瘤特异性死亡”但大量患者死于心血管疾病。普通 KM 曲线会把心血管死亡当作删失相当于假设这些患者“仍然有机会在未来发生目标事件”但现实中他们已经不可能再发生肿瘤死亡了。这会让目标事件的累积发生率被高估。情况二HR 不恒定。Cox 回归输出的是一个平均效应。如果治疗在早期有效、晚期无效或者某种治疗前期风险高、后期收益大单个 HR 就会掩盖这种时间动态变化。很多人一看到cox.zph()的 p 值小于 0.05 就慌了不知下一步该怎么处理。情况三只有相对效应缺少绝对效应。HR0.72 在统计上是显著的但临床上医生更想问的是“用新方案患者 3 年内平均能多活几个月”这一句话 HR 回答不了需要 RMST。一句话总结基础与进阶方法的分工基础方法回答“有没有关系”进阶方法回答“这个关系是否恒定、是否被其他事件干扰、临床上到底有多大”。2. 生存分析核心概念先花三分钟统一语言2.1 一条生存数据由什么构成进入进阶方法之前先统一概念。一条标准的生存数据有三个要素要素含义例子时间 T从起点到事件发生或删失的时间手术日期到死亡日期事件状态 status是否发生事件以及事件类型0删失1目标事件2竞争事件协变量 X可能影响生存的特征年龄、性别、治疗组、分期普通 Cox 回归要求status是 0/1 二值变量。到竞争风险模型时status会变成 0/1/2 多分类这是很多人第一次进阶时卡住的点数据结构里根本没有为“第二种事件”预留位置。2.2 生存函数、风险函数与CIF几个必须分辨的指标生存函数 S(t)到时间 t 为止仍未发生事件的概率KM 曲线估计的就是它。风险函数 h(t)在 t 时刻存活的个体中瞬时发生事件的概率。累积发生率函数 CIF(t)到时间 t 为止某个特定原因事件累积发生的概率。在只有单一终点时CIF 和 1 - KM 是等价的。但存在竞争风险时两者不再相等。普通 KM 把竞争事件当作删失等于把已经从风险集中移除的患者又“保留”在分母里导致目标事件发生率被高估。这也是为什么竞争风险场景下必须画 CIF而不是画 1 - KM。2.3 删失与竞争事件的区别删失censor和竞争事件competing event在数学上都表现为“观察不到目标事件”但含义完全不同删失患者在随访截止时仍未发生事件或者失访。理论上如果他继续被观察未来仍可能发生目标事件。竞争事件患者先发生了另一种事件使目标事件在真实世界中无法再发生。比如患者死亡后不可能再“复发”。正是因为竞争事件不是删失所以不能用普通 Cox 的删失逻辑去处理。下一章进入第一板斧。3. 第一板斧竞争风险模型3.1 什么时候必须考虑竞争风险判断标准很简单研究期间内是否有一种事件会发生并且它的发生会阻止目标事件被观察到。典型场景包括肿瘤长期随访中非肿瘤死亡与肿瘤死亡竞争心血管研究中非心血管死亡与心血管事件竞争骨科植入物研究中假体失效与患者死亡竞争老年队列研究里任何“死亡”都可能与“入院”“跌倒”等事件竞争。如果竞争事件占比很低比如不足 5%用普通方法影响不大但如果随访时间长、患者年龄大竞争事件比例往往达到 20% 甚至更高这时候不做竞争风险分析结论很容易被质疑。3.2 原因别风险与子分布风险别再混淆竞争风险的建模有两套主流框架很多人一直分不清模型估计对象风险集定义适合回答的问题原因别风险模型Cause-specific每种事件的瞬时风险仍然存活、尚未发生任何事件的人群某种干预对目标事件的病因学效应Fine-Gray 子分布风险模型事件对应的 CIF 累计发生率包含已发生竞争事件的患者预测、预后分层、临床决策通俗理解原因别模型关心“在还活着的人里这个因子的致病作用有多大”Fine-Gray 模型关心“到某个时间点这个人群里有多少比例的人发生了目标事件并且竞争事件会拉低这个比例”。因此两者结果不一致很正常。原因别 HR 显著而 Fine-Gray HR 不显著说明该因子确实增加了目标事件的生物学风险但患者更容易先死于竞争事件反过来Fine-Gray 显著而原因别不显著常见于竞争事件本身受该因子影响较大的情况。报告时最好把两种模型都给出并说明研究定位是病因探索还是预后预测。3.3 R代码CIF、原因别Cox与Fine-Gray进阶生存分析的标准实现目前以 R 生态最成熟。这里用经典cmprsk包演示。假设数据有三列time为随访时间status取 0删失、1目标事件、2竞争事件以及协变量age、sex、treat。# 安装并加载包 install.packages(cmprsk) library(cmprsk) # 读取数据确认status编码 dat - read.csv(survival_data.csv) table(dat$status) # 0删失, 1目标事件, 2竞争事件 # 第一步绘制两组的目标事件CIF cif - cuminc( ftime dat$time, fstatus dat$status, group dat$treat ) # 绘图手动指定曲线图例 plot(cif, curvlab c(A组-删失, A组-目标事件, A组-竞争事件, B组-删失, B组-目标事件, B组-竞争事件))注意cuminc()默认把status0当作删失编码也就是参数cencode0。如果你项目里删失编码不是 0需要显式指定。接下来分别拟合两种模型# 第二步原因别Cox模型对目标事件建模竞争事件视为删失 cs_model - coxph( Surv(time, as.numeric(status 1)) ~ age sex treat, data dat ) summary(cs_model) # 第三步Fine-Gray子分布风险模型 fg_model - crr( ftime dat$time, fstatus dat$status, cov1 model.matrix(~ age sex treat, data dat)[, -1] ) summary(fg_model)crr()的cov1必须是矩阵形式所以用model.matrix()转换。默认failcode1表示把status1当作目标事件cencode0表示把status0当作删失与前面的数据结构一致。3.4 Fine-Gray与原因别模型的结果取舍实际项目中我更推荐这样的报告策略先用cuminc()画出各组 CIF并在关键随访时间点如 12、24、36 个月报告 CIF 估计值组间 CIF 差异用 Grays 检验cuminc()返回结果里的Tests字段就是各组事件的 p 值论文正文同时给出原因别 Cox 和 Fine-Gray 的结果避免审稿人质疑方法选择如果研究目的是预测个体预后以 Fine-Gray 为主如果是探索病因机制以原因别模型为主。需要提醒的是Fine-Gray 模型的“风险集”包含了已发生竞争事件的患者这在概率上是一个有意的构造用来直接建模 CIF。正因如此它的结果不能直接解释为“活着的患者中某个因素导致的瞬时风险”这是它被质疑最多的地方。4. 第二板斧时变效应解决HR不恒定4.1 PH假设到底意味着什么Cox 回归的前提是比例风险假设Proportional HazardsPH任意两个个体的风险函数之比在整个随访期内保持不变。也就是说治疗组的 HR 在随访第 1 个月、第 12 个月、第 36 个月都是同一个数值。这个假设在现实中经常被违背。例如免疫治疗起效慢前期可能看不出差异后期优势才显现手术干预早期风险高于保守治疗但长期预后更好某种药物的获益随时间衰减HR 从 0.6 逐渐趋近于 1。如果 PH 不成立直接报告一个 HR 会严重误导读者因为这一现象背后的真实效应是时变的。4.2 用Schoenfeld残差检验PH假设我们先用标准 Cox 模型拟合再检验 PH 假设。R 中使用cox.zph()对 Schoenfeld 残差做检验同时画出残差随时间的变化library(survival) cox_model - coxph(Surv(time, status) ~ age sex treat, data dat) summary(cox_model) # 检验PH假设 zph_result - cox.zph(cox_model) print(zph_result) # 查看某个变量残差随时间的趋势 plot(zph_result[3])cox.zph()输出中每个协变量对应一个 p 值p 值小于 0.05 说明该变量的效应随时间变化PH 假设不成立。这里的截断值可以放宽到 0.05但如果样本量很大即使很小的偏差也会显著因此还要结合残差图来判断是否存在实际意义上的趋势。4.3 时变系数模型tt()实战确认 PH 假设不满足后最直接的解决办法之一是对该变量引入时间交互项即“时变系数”模型。在survival包中通过tt()参数实现# 对treat做log(t1)交互模拟HR随时间对数衰减 cox_tt - coxph( Surv(time, status) ~ age sex treat tt(treat), data dat, tt function(x, t, ...) { x * log(t 1) } ) summary(cox_tt)这时treat在时间 t 的效应不再是常数而是HR(t) exp(β_treat β_tt * log(t 1))其中β_tt显著就说明治疗效应确实随时间变化。tt()里的变换函数可以根据数据形态替换常见的有线性t、对数log(t1)、以及分段函数。选择变换函数时最好先看 Schoenfeld 残差图是线性趋势、对数趋势还是折线趋势再决定用哪种形式不要盲目套用。另一条常见路径是分层。当某个变量本身不是关注焦点但该变量违反 PH 假设时用strata()把它分层可以得到基线风险函数不同、但其他变量仍有 HR 的模型cox_stratified - coxph( Surv(time, status) ~ age sex strata(treat), data dat ) summary(cox_stratified)分层后treat不再输出 HR但其他协变量估计更稳健。还可以用时间分段模型把随访期切成几段分别估计不同阶段的 HR。R 中使用survSplit()把数据转换成含(tstart, time]区间的长格式dat_split - survSplit( Surv(time, status) ~ ., data dat, cut c(12, 24), episode period ) cox_period - coxph( Surv(tstart, time, status) ~ treat treat:strata(period) age sex, data dat_split ) summary(cox_period)这种方法直观但切点选择要谨慎最好在研究方案里预先定义而不是根据 p 值来回试。4.4 时变协变量与常见误区上一节解决的是“协变量效应随时间变化”还有一类问题是“协变量本身随时间变化”比如患者的血压、用药剂量、疾病状态。这类问题的标准做法是把数据改成长格式用(start, stop]区间记录每个患者在不同时间段的值。Python 用户可以使用lifelines的CoxTimeVaryingFitterfrom lifelines import CoxTimeVaryingFitter # 长格式数据id_col标识患者start/stop为区间status为区间内是否事件 df_long ... ctv CoxTimeVaryingFitter() ctv.fit( df_long, id_colid, start_colstart, stop_colstop, event_colstatus, show_progressTrue ) ctv.print_summary()这里有一个必须反复强调的坑内生性时变协变量会引入严重偏倚。如果某个协变量受治疗本身影响又作为时变协变量放进模型就会形成因果倒置。比如把“患者是否发生并发症”当作时变协变量预测死亡但治疗组更容易发生某种并发症模型中就会错误地吸收治疗效应。另外时变协变量的引入还要警惕永恒时间偏倚immortal time bias。如果患者在某个事件发生之前不可能被暴露比如“术后接受了二期治疗”这个变量患者在治疗开始前的一段时间天然不可能死亡如果把这部分时间错误归入治疗组会得到偏高的保护效应。解决思路是设计阶段就用 landmark 分析或恰当的区间切分而不是建模时补救。5. 第三板斧RMST把HR翻译成人话5.1 RMST解决什么问题HR 是一个相对风险临床医生很难直接感知道“HR0.72”意味着什么。而且当 PH 假设不成立时Cox 模型输出的单个 HR 本质上是整个随访期效应的某种加权平均解释起来更加危险。RMSTRestricted Mean Survival Time限制平均生存时间的概念非常简单从 0 到某个指定时间点 tau 之间生存曲线下的面积。它表示患者在 tau 时间内平均存活的月数。举个例子两组患者随访 36 个月A 组 RMST28.5 个月B 组 RMST24.3 个月差值是 4.2 个月。这句话临床医生一听就懂A 组患者在 3 年内平均比 B 组多活了 4.2 个月。RMST 不需要 PH 假设对生存曲线尾部不稳定的问题也更稳健因此在长期随访和肿瘤临床试验中越来越受重视。5.2 tau怎么选才不被审稿人挑战RMST 的估计强烈依赖 tau 的选择因为生存曲线在 tau 之后的部分被截断掉了。推荐做法在统计分析计划SAP中提前定义 tau通常是研究中位随访时间、试验设计的最短随访时间或临床上有意义的观察窗口不要看了数据后再挑一个让结果最漂亮的 tau这属于 p-hacking报告时应做敏感性分析选择 tau±6 个月或一组合理的 tau 值看结论是否稳定。5.3 R代码survRM2的两组RMST比较R 中survRM2包是最常用的实现支持两组比较和基于伪观测法的协变量调整install.packages(survRM2) library(survRM2) # 两组RMST比较tau36个月 rmst_fit - rmst2( time dat$time, status dat$status, arm dat$treat, tau 36 ) print(rmst_fit) # 调整协变量的RMST比较 rmst_adj - rmst2( time dat$time, status dat$status, arm dat$treat, tau 36, covariates dat[, c(age, sex)] ) print(rmst_adj)这里status仍使用 0/1 编码RMST 关心的是“事件发生前平均生存时间”所以竞争风险场景下要把目标事件单独编码。如果竞争风险占比很高更严谨的做法是在竞争风险框架下比较 CIF再结合 RMST 做辅助解释。5.4 Python实现RMSTlifelines目前没有直接提供 RMST 函数但根据定义我们可以用 KM 曲线做数值积分。这条代码可以用在任何 Python 项目中import numpy as np from lifelines import KaplanMeierFitter def rmst_from_km(time, event, tau): 基于KM生存曲线的RMST计算梯形法近似面积 kmf KaplanMeierFitter() kmf.fit(time, event_observedevent) times kmf.survival_function_.index.values surv kmf.survival_function_[KM_estimate].values # 截取到tau mask times tau tt np.concatenate([[0], times[mask], [tau]]) ss np.concatenate([[1], surv[mask], [surv[mask][-1] if np.any(mask) else 1]]) return np.trapz(ss, tt) # 用法示例 rmst_a rmst_from_km(df_a[time], df_a[status], tau36) rmst_b rmst_from_km(df_b[time], df_b[status], tau36) print(fA组RMST{rmst_a:.2f}月, B组RMST{rmst_b:.2f}月) print(fRMST差值{rmst_a - rmst_b:.2f}月)手动实现的好处是逻辑透明方便集成到自己的 Python 分析流水线中缺点是没有直接给出置信区间和 p 值。正式分析时建议两两对照 R 的survRM2结果确认数值一致后再用于论文。5.5 与HR的关系和报告当 PH 假设成立时RMST 的结论通常与 HR 方向一致。但 RMST 提供的是绝对差异HR 提供的是相对风险两者是互补关系。最佳报告方式是“HR RMST差值”同时给出表格里列出中位随访时间、各组事件数、HR95% CI、RMST95% CI和 RMST 差值正文里重点解释 RMST 差值的临床含义如果 PH 不成立明确说明 HR 仅作参考主要结论以 RMST 差异或时变效应结果为准。6. 三板斧到底怎么选一个决策框架讲完三板斧很多人会问是不是以后所有分析都要全部上不是。方法论必须服务于研究问题。下面是我在项目中常用的判断流程数据场景首选方法补充方法单终点、无竞争事件、PH成立KM Cox可加 RMST 辅助随访期长、存在非目标死亡CIF 竞争风险模型RMST治疗效果随时间变化时变系数模型 / 分层 / 时间分段分阶段报告 HR需要解释临床绝对获益RMSTHR KM预后预测、个体风险分层Fine-Gray机器学习生存模型更具体的判断逻辑先看研究问题。是探究病因还是预测预后预测优先 Fine-Gray病因优先原因别模型。再画图。CIF 曲线之间是否分离、是否交叉Schoenfeld 残差图是否有明显趋势该上 RMST 就上。当随访时间确定、临床意义是关键