生存分析进阶：竞争风险模型、时依协变量与RMST实战指南

从最基础的 Kaplan-MeierKM曲线和 Cox 比例风险回归入手几乎成了每个入门生存分析的人必经之路。KM 曲线画出来简单直观Cox 回归一行代码也能跑出 HR风险比和 p 值。但一旦进入真实项目尤其是临床随访数据、用户流失预测、设备故障分析这类场景你很快会碰到几个基础方法解释不了的问题有人还没等到目标事件就先发生了别的结局某个关键因素在随访过程中变了或者 HR 虽然显著却很难告诉业务方“到底多活了多少个月”。这篇文章就来拆解这三个进阶方向分别用可运行代码演示竞争风险模型、时依协变量 Cox 模型和 RMST限制性平均生存时间的计算与解读。读完你既能补上基础方法之外的短板也能直接把这些代码改改用到自己的数据上。1. 为什么说 KM 曲线 Cox 回归只是起点1.1 基础方法的定位与优势KM 曲线是估计生存概率的非参数方法。它把随访时间切成一堆事件发生时刻在每个时刻用“当前风险集人数”和“事件发生数”计算条件生存概率再连乘得到整个生存曲线。它不需要假设生存时间服从什么分布所以应用范围很广。对于两组生存曲线的比较通常配合 log-rank 检验。Cox 回归则是半参数模型它的核心假设是比例风险PH 假设也就是不同组别的风险函数成比例比值不随时间变化。Cox 模型的优点是可以同时纳入年龄、治疗、临床指标等多个协变量输出每个因素的风险比 HR 和置信区间是生存分析中最常用的多因素回归工具。这两件套在入门阶段性价比极高。数据规整、终点单一、删失机制独立、PH 假设成立时它们已经够用。很多医学论文、用户留存分析、设备可靠性分析都用这套流程拿到核心结论。1.2 基础方法解决不了的三类现实问题第一类问题终点不只有一个。比如研究“疾病复发”但一部分患者先死亡了。死亡之后再也不可能观察到复发两种终点之间存在竞争关系。如果像 KM 曲线那样把死亡当作普通删失处理会系统性高估复发的累积发生率因为删失并不是随机的它和“复发”这个终点有生物学上的关联。这种场景需要用竞争风险模型。第二类问题协变量会随时间变化或者 PH 假设不成立。基线 Cox 模型默认一个病人的协变量值是入组时测量的之后不再改变。但现实中用药方案会调整、血压会波动、用户是否活跃会变化。这些时变信息如果被忽略就会丢失重要信号。另一种情况是某个分组变量对风险的影响随时间衰减或增强导致风险比不再是常数这时需要更灵活的建模方式。第三类问题HR 不容易直接转化为临床或业务收益。HR 0.7 只说明处理组风险比对照组低 30%但业务方更想知道“平均能多存活多少个月”“在 3 年窗口期内能多争取多久的无事件时间”。生存曲线下面积也就是 RMST能给出更直观的答案。1.3 三板斧分别解决什么这篇文章要讲的进阶三板斧正好对应上面三类问题板斧核心方法解决了什么问题第一板斧竞争风险模型多种终点互相竞争时正确估计累积发生率并回归分析影响因素第二板斧时依协变量 Cox 模型协变量随时间变化或 PH 假设不成立时充分利用动态信息第三板斧RMST用受限平均生存时间作为效应指标直观回答“获益多少时间”这三板斧不是替代 KM 和 Cox而是在它们的基础上补足盲区。下面我们用同一份模拟数据把每一步走通。2. 环境准备与数据说明2.1 Python 环境与依赖本文代码基于 Python 3.9推荐在 Jupyter Notebook 或任意 Python 脚本中运行。需要安装以下库pip install lifelines scikit-survival numpy pandas matplotlib版本方面建议使用当前各库的较新稳定版本。lifelines 负责 KM 曲线、Cox 回归和时依协变量 Cox 模型scikit-survivalsksurv负责竞争风险模型。安装完成后可以先确认导入是否正常import numpy as np import pandas as pd import matplotlib.pyplot as plt import lifelines from lifelines import KaplanMeierFitter, CoxPHFitter, CoxTimeVaryingFitter import sksurv from sksurv.linear_model import CompetingRiskSurvivalAnalysis print(lifelines:, lifelines.__version__) print(sksurv:, sksurv.__version__)不同版本的 API 可能会有细微调整如果某个方法名在当前版本中提示不存在请以官方文档为准。2.2 模拟一份带竞争终点的随访数据为了方便演示我们构造一份 300 条记录的模拟数据。每条记录对应一个研究对象随访终点状态定义为0删失也就是随访结束或失访时目标事件没有发生1发生“复发”事件这是我们要研究的主要终点2发生“死亡”事件这是和复发竞争的终点。模拟数据中包含两个协变量年龄age和是否接受新药治疗treatment。复发时间受治疗影响治疗组复发风险更低死亡时间受年龄影响年龄越大死亡风险越高。删失时间随机生成。最终观测时间是复发时间、死亡时间、删失时间三者中的最小值同时为了模拟真实随访场景所有超过 60 个月的观测都按 60 个月截断删失处理。np.random.seed(2024) n 300 age np.random.normal(55, 10, n).clip(30, 85) treatment np.random.binomial(1, 0.5, n) # 复发风险治疗降低复发风险年龄增大复发风险略升 rec_scale 30 / np.exp(0.7 * treatment 0.02 * (age - 50)) # 死亡风险年龄越大风险越高 death_scale 45 / np.exp(0.04 * (age - 50)) # 删失时间 cens_scale 50 t_rec np.random.exponential(scalerec_scale) t_death np.random.exponential(scaledeath_scale) t_cens np.random.exponential(scalecens_scale) T_comp np.minimum(np.minimum(t_rec, t_death), t_cens) E_comp np.where(T_comp t_rec, 1, np.where(T_comp t_death, 2, 0)) # 超过60个月统一按删失处理 T np.minimum(T_comp, 60) E np.where(T_comp 60, 0, E_comp) df_competing pd.DataFrame({ patient_id: range(1, n 1), months: T, event: E, age: age, treatment: treatment }) print(df_competing[event].value_counts()) print(df_competing.groupby(event)[months].describe())运行后可以观察一下各状态的数量。事件 1 和事件 2 都会占到一定比例删失占比通常在三成左右。这份数据就是我们后续所有演示的基础。3. 第一板斧竞争风险模型3.1 什么是竞争风险竞争风险指的是在研究主要终点时存在另一个事件会阻碍主要终点被观察。经典的例子是研究肿瘤患者的复发死亡就是竞争事件研究心脏病发作非心血管死亡就是竞争事件。传统做法是把死亡当作删失用 KM 法估计“复发概率”。但这样做的隐含假设是删失和复发独立。死于其他原因的人和仍然在风险集中的人未来复发概率相同。这个假设在临床场景中往往不成立因为年纪大、身体差的人既更容易死亡也更容易复发。把死亡当成删失会让风险集里“不太容易复发”的人被剔除最终高估累计复发率。3.2 用 KM 高估复发率的直观演示我们先按传统方式用 KM 曲线估计复发累积发生率。做法是把死亡事件全部转成删失df_competing[event_rec_km] (df_competing[event] 1).astype(int) kmf_rec KaplanMeierFitter() kmf_rec.fit( df_competing[months], df_competing[event_rec_km], labelKM 复发(死亡算删失) )然后估计正确的累积发生率函数 CIF。CIF 的计算思路和 KM 不同它把竞争事件也纳入风险集计算得到的才是“到某个时间点复发发生的概率”。这里我们手动实现一个 CIF 估计函数方便看清楚每一步def cif_at_time(df, event_of_interest, t): 计算在时间 t 内指定事件 event_of_interest 的累积发生率 CIF(t)。 思路在每个事件时刻用总风险集和事件 k 的发生数 乘以“之前未发生任何事件”的总生存概率。 # 总体生存概率任意事件都不算删失 df_temp df.copy() df_temp[any_event] (df_temp[event] ! 0).astype(int) kmf_any KaplanMeierFitter() kmf_any.fit(df_temp[months], df_temp[any_event]) event_table kmf_any.event_table # 包含 at_risk, observed 等列 # 只收集目标事件的唯一发生时间 target_times np.sort(df.loc[df[event] event_of_interest, months].unique()) cif 0.0 for event_t in target_times: if event_t t: break at_risk (df[months] event_t).sum() d_k ((df[months] event_t) (df[event] event_of_interest)).sum() # 在 event_t 之前未发生任何事件的生存概率即 S(event_t - 0) prev_times event_table.index[event_table.index event_t] S_before 1.0 if len(prev_times) 0: S_before ( 1 - event_table.loc[prev_times, observed] / event_table.loc[prev_times, at_risk] ).prod() cif d_k / at_risk * S_before return cif有了 CIF 函数我们就可以在时间轴上计算复发累积发生率并和 KM 做对比time_grid np.arange(0, 61, 1) cif_values [cif_at_time(df_competing, 1, t) for t in time_grid] km_survival kmf_rec.survival_function_at_times(time_grid) km_cumulative 1 - km_survival.values.flatten() plt.figure(figsize(8, 5)) plt.step(time_grid, cif_values, wherepost, labelCIF 复发(竞争风险正确估计)) plt.step(time_grid, km_cumulative, wherepost, label1 - KM 生存(死亡当删失)) plt.xlabel(随访时间月) plt.ylabel(累积发生率) plt.title(竞争风险下KM 高估复发风险) plt.legend() plt.grid(alpha0.3) plt.show()从图中可以看到在随访后期两条曲线的差距会逐渐变大。这说明当死亡风险不可忽视时1 - KM 生存概率并不是复发概率的一致估计必须用 CIF 替换。3.3 Fine-Gray 子分布风险回归CIF 是非参数估计适合做描述。如果需要同时分析年龄、治疗等因素对复发累积发生率的影响可以使用 Fine-Gray 子分布风险回归模型也就是 subdistribution hazard model。它建模的对象不是瞬时风险而是累积发生函数对应的子分布风险回归系数可以解释为子分布风险比SHR。scikit-survival 中的CompetingRiskSurvivalAnalysis可以完成这一任务。在准备数据时需要把生存数据转成结构化数组event字段用整数编码0 表示删失1 和 2 分别表示两种事件。X df_competing[[age, treatment]].copy() y np.array( [(e, t) for e, t in zip(df_competing[event], df_competing[months])], dtype[(event, i4), (time, f8)] ) fg_model CompetingRiskSurvivalAnalysis() fg_model.fit(X, y) result pd.Series(fg_model.coef_, indexX.columns) result pd.DataFrame(result, columns[coef_]) result[SHR] np.exp(result[coef_]) print(result.round(3))以治疗因素为例如果treatment的系数是负值说明治疗组比对照组的复发子分布风险更低。这里的解释要特别谨慎它对应的是累积发生率尺度上的组间差异而不是传统 Cox 模型中“复发瞬时风险”的差异。两个事件类型之间如何取舍CompetingRiskSurvivalAnalysis会同时考虑所有编码为非 0 的事件。如果你想分别研究“复发”和“死亡”的 subdistribution hazard可以在完整数据上拟合后从业务角度选择解释对象。更多时候我们会针对主要终点单独解释回归系数把其他终点作为竞争事件纳入计算。3.4 竞争风险模型结果解读Fine-Gray 模型输出的 SHR 与 Cox 模型输出的 HR 含义不同。HR 描述的是瞬时风险比SHR 描述的是子分布风险比直接对应累积发生率。临床文献中常见的表述方式是调整年龄后治疗组 3 年累积复发风险比对照组低 XX%。这种表述更能回答患者和医生真正关心的问题。在模拟数据上治疗因素应该表现出明显的保护作用年龄也会显著影响死亡竞争风险。如果你的实际数据中主要终点和竞争终点之间没有相关性那么传统 KM 与 CIF 的结果差距不会太大差异越大说明竞争风险的影响越不可忽视。4. 第二板斧时依协变量 Cox 模型4.1 基线 Cox 的局限性CoxPHFitter 默认把每个个体的协变量当作基线值也就是随访开始时的测量值。但这种设定在处理“治疗过程中换药”“用户从活跃变为沉默”“环境暴露强度变化”等场景时会丢失重要的动态信息。另一方面如果某个协变量对风险的影响不满足 PH 假设比如治疗早期效果明显、后期效果减弱那么固定值的 Cox 模型也会有问题。虽然我们可以用分层、加时间交互项来修补但更通用的做法之一是构建时依协变量模型。当然严格来说“时依协变量”和“时依效应”是两个不同的问题前者是协变量取值随时间变化后者是协变量效应随时间变化但它们的建模载体都可以是时间分段数据。4.2 准备长格式数据时依协变量 Cox 模型要求数据按“区间”展开。每个个体可能有多行记录每行对应一个时间区间用start和stop标记区间的起止时间event标记该区间结束时是否发生事件。协变量取值在该区间内保持不变。下面构造一份示例数据。基线数据沿用上一节生成的df_competing但额外模拟一个“治疗状态切换”过程部分患者在随访中某一天切换了治疗状态比如从标准治疗切换到新药或者停止用药。我们记录每个人的切换时间switch_time然后按区间拆分成长格式。np.random.seed(88) switch_time np.random.gamma(shape2, scale7, sizen) long_rows [] for i, row in df_competing.iterrows(): pid row[patient_id] stop_time row[months] event row[event] sw switch_time[i] if sw stop_time: # 随访期内未切换只需要一行 long_rows.append([pid, 0.0, stop_time, event, row[treatment], row[age]]) else: # 切换前区间无事件 if sw 1e-6: long_rows.append([pid, 0.0, sw, 0, row[treatment], row[age]]) # 切换后区间事件可能发生在这里 long_rows.append([pid, sw, stop_time, event, 1 - row[treatment], row[age]]) df_long pd.DataFrame( long_rows, columns[patient_id, start, stop, event, current_treatment, age] ) print(df_long.head(10))这里有几个需要特别注意的地方start和stop必须是数值且每个区间内stop start。event只在最后一个区间能为 1 或 2中间区间的事件必须为 0。同一个体的多个区间靠patient_id关联区间不能重叠且要覆盖完整的随访时间。如果某个区间的长度几乎为 0模型会不稳定需要在构造数据时剔除或合并。4.3 拟合 CoxTimeVaryingFitterlifelines 中的CoxTimeVaryingFitter专门处理这种长格式数据。它不再要求 PH 假设覆盖整个随访期而是在每个短区间内近似满足比例风险。ctv CoxTimeVaryingFitter(penalizer0.01) ctv.fit( df_long, id_colpatient_id, start_colstart, stop_colstop, event_colevent, show_progressTrue ) ctv.print_summary()输出结果中会看到current_treatment和age的对数风险系数、exp(coef) 和 p 值。current_treatment的含义是在当前区间对应的治疗状态下单位变化对瞬时风险的影响。相比于基线 Cox 模型中的treatment它能捕捉“治疗状态变化后风险随之变化”这一动态过程。4.4 与基线 Cox 模型对比为了突出时依信息的价值我们用同样的数据拟合一个基线 Cox 模型但只使用基线治疗状态忽略切换信息df_base_cox df_competing.copy() df_base_cox[status_any_event] (df_base_cox[event] ! 0).astype(int) cph_baseline CoxPHFitter() cph_baseline.fit( df_base_cox[[months, status_any_event, treatment, age]], duration_colmonths, event_colstatus_any_event, show_progressTrue ) cph_baseline.print_summary()再与时依模型结果对比。你会发现如果切换治疗确实对复发或死亡风险有影响那么忽略切换的基线模型会低估治疗效果或产生偏倚。这里要特别提醒切莫把“时依协变量模型”当成万能工具如果协变量本身是事件结局的前兆或者由结局倒推而来会引入内生性偏倚。这在用药依从性分析中尤其常见所谓的“健康用户效应”会让药物看起来格外有效。5. 第三板斧RMST5.1 什么是 RMSTRMST全称 Restricted Mean Survival Time限制性平均生存时间。它的定义是生存曲线从 0 到某个时间点 tau 之间的曲线下面积。直觉上RMST(tau) 表示“在 tau 时间窗口内平均能够存活或无事件的时间”。之所以要加一个“限制性”是因为在实际随访中很难完整观察所有个体的生存时间。随访到头时很多人仍然存活或尚未发生事件。如果直接计算平均生存时间这些删失个体无法简单赋值。RMST 把关注窗口限制在 [0, tau] 内在这个范围内生存曲线可以被 KM 估计出来曲线下面积自然也可以计算。5.2 计算两组 RMST 差异我们继续用模拟数据研究“治疗组”与“对照组”在无事件生存时间上的差异。为了同时体现竞争风险的思路这里暂以“复发”作为目标事件把死亡当作删失来演示 KM 与 RMST 的匹配更严格的竞争风险场景下可以将生存函数替换为 CIF 对应的无事件概率。tau 48 def rmst_estimate(df, treatment_value, tau, event_colevent_rec_km): sub df[df[treatment] treatment_value] kmf KaplanMeierFitter() kmf.fit(sub[months], sub[event_col]) time_grid np.linspace(0, tau, 2000) surv kmf.survival_function_at_times(time_grid).values.flatten() return np.trapz(surv, time_grid) rmst_control rmst_estimate(df_competing, 0, tau) rmst_treat rmst_estimate(df_competing, 1, tau) print(f对照组 RMST({tau}个月) {rmst_control:.2f} 个月) print(f治疗组 RMST({tau}个月) {rmst_treat:.2f} 个月) print(fRMST 差异 {rmst_treat - rmst_control:.2f} 个月)结果可以这样解释在 48 个月的窗口内治疗组平均比对照组多了 X 个月的无事件生存时间。这个数字比 HR 更容易被业务方和临床医生理解。5.3 用 Bootstrap 计算置信区间要判断 RMST 差异是否统计显著需要计算置信区间。常见做法是用 Bootstrap 重采样重复估计多次取 2.5% 和 97.5% 分位数作为置信区间。def rmst_diff_bootstrap(df, tau, n_boot1000, event_colevent_rec_km): diffs [] for _ in range(n_boot): boot df.sample(frac1.0, replaceTrue) rmst_0 rmst_estimate(boot, 0, tau, event_col) rmst_1 rmst_estimate(boot, 1, tau, event_col) diffs.append(rmst_1 - rmst_0) diffs np.array(diffs) return np.percentile(diffs, 2.5), np.percentile(diffs, 97.5) ci_low, ci_high rmst_diff_bootstrap(df_competing, tau) print(fRMST 差异 95% Bootstrap 置信区间: ({ci_low:.2f}, {ci_high:.2f}))如果置信区间不包含 0说明两组在窗口期内的平均无事件生存时间差异显著。实际分析中Bootstrap 次数可以增加到 5000 或 10000结果会更稳定。注意每个 Bootstrap 样本中如果某一组的样本量过小KM 曲线会不稳定这时需要谨慎解释。5.4 RMST 与 HR 的关系HR 和 RMST 并不矛盾只是从两个角度回答问题。HR 关心“风险强度”RMST 关心“时间获益”。举例来说HR 0.7 说明治疗组风险比对照组低 30%但患者更想知道“我平均能多活几个月”后者只有在 RMST 的框架下才直观。更关键的是RMST 不依赖 PH 假设。当 PH 假设明显不成立时HR 是一个随时间变化的混合指标解读困难而 RMST 在任何时候都能计算。因此近年来很多临床研究把 RMST 作为次要效应指标甚至替代 HR 作为主要指标。分析实践中建议同时报告 HR 和 RMST 差异互为补充。6. 三板斧如何选型场景推荐方法关键输出指标单一终点PH 假设成立KM Cox生存概率HR存在竞争终点竞争风险模型CIFSHR协变量随时间变化时依协变量 Cox动态 HRPH 假设不成立但需要解释组间差异含时间交互的 Cox 或 RMST时依效应RMST 差异业务方要求“获益多少时间”RMSTRMST 差异及置信区间组合使用时一般先画 CIF 或 KM 曲线看整体趋势再做多因素回归找到核心影响因素最后用 RMST 补充效应量和时间获益的解读。7. 常见问题与排查思路问题现象常见原因解决思路sksurv 拟合报错event 字段类型不对或编码不是整数使用结构化数组 dtype[(event, i4), (time, f8)]0 表示删失CIF 曲线在后期震荡某个时间点风险集太小合并时间段或把 tau 限制在数据支撑范围内CoxTimeVaryingFitter 报“overlapping intervals”同一 id 的区间存在重叠或重复检查 start/stop 是否连续中间区间 event 是否为 0时依模型系数极不稳定区间过短或协变量变化过于频繁增加区间最小长度或对时变协变量做平滑分组RMST 在不同 tau 下结论不一致tau 选择对结果影响大选择有临床或业务意义的时间窗口做敏感性分析竞争风险模型下 HR 与 SHR 方向不一致两种模型定义不同明确报告模型类型避免混用排查时最好先画描述性图比如 CIF 曲线、生存曲线、Schoenfeld 残差图找到数据中的异常模式再进入建模环节。8. 最佳实践与工程建议8.1 数据清洗阶段生存分析的数据质量决定模型上限。首先要明确终点的定义特别是“删失”的判定口径。其次要确认随访时间的起点和终点避免 immortal time bias。所谓 immortal time是指研究开始后到治疗开始前的这段“必然存活”的时间如果处理不当会夸大治疗效果。时依协变量模型是处理这类问题的常用手段之一。8.2 模型验证竞争风险模型和时依协变量模型同样需要验证。对于 Cox 系列模型Schoenfeld 残差检验可以帮助判断 PH 假设是否成立。如果 PH 假设明显不成立可以考虑分层 Cox、加入时间交互项或切换为 RMST 框架。对于 Fine-Gray 模型可以通过校准曲线评估预测的 CIF 与实际观察 CIF 的一致性。样本量允许时使用交叉验证或 Bootstrap 评估模型稳定性。8.3 报告与生产落地学术报告或交付文档中应该同时包含各组事件状态的数量分布、主要终点 CIF 曲线、回归系数与置信区间、RMST 差异及 tau 的选择依据。如果模型会投入业务预测系统还要额外关注模型的输入特征是否在生产环境能够实时获取时依协变量的数据延迟会不会影响预测时效模型更新频率以及是否需要周期性重训练。安全方面涉及真实患者或用户数据时必须经过合法授权和脱敏处理。任何删除数据或修改随访状态的操作都要在测试环境中充分验证并保留备份。9. 总结与下一步本文围绕生存分析中的三个进阶方向展开竞争风险模型解决多终点竞争问题时依协变量 Cox 模型解决动态因素建模问题RMST 解决效应量解读问题。配套的三段可运行代码覆盖了从数据构造、模型拟合到结果解读的完整流程。直接改换成自己的数据时重点检查事件编码、区间拆分和窗口期选择这几个关键环节。下一步可以继续学习多状态模型Multi-State Model、动态预测Dynamic Prediction以及基于机器学习的生存模型例如随机生存森林和 DeepSurv。这些方法在复杂场景中各有优势但底层逻辑仍然离不开今天说的 CIF、时依风险和受限平均生存时间这些核心概念。建议先把手头数据用这三板斧完整跑一遍理解每一步的输出含义再去扩展到更复杂的模型。如果这篇文章对你有帮助可以收藏备用后续遇到生存分析中的竞争风险和时变因素问题时直接翻出来对照着排查。