生存分析进阶三板斧：时依Cox、竞争风险与机器学习

KM曲线和Cox回归，几乎是生存分析的两张默认名片。很多分析流程都是先画KM曲线做Log-rank检验，再跑一个Cox回归报告HR和95%CI。这个套路在简单场景下够用，但一旦遇到PH假定不满足、存在竞争风险、特征多且关系非线性，结果就可能失真甚至方向相反。

这次我们来看三类更进阶的生存分析手段，我习惯叫它“生存分析进阶三板斧”：

• 第一板斧：时依协变量Cox模型，专门处理PH假定不满足的情况。
• 第二板斧：竞争风险模型，处理“患者先发生其他事件，导致目标事件无法被观测到”的问题。
• 第三板斧：机器学习生存分析，用随机生存森林、生存SVM、DeepSurv等模型处理高维特征和非线性关系。

文章会从适用场景、R/Python实现、结果解读、常见报错和论文报告习惯五个方面展开。适合正在做临床数据分析、生信挖掘、用户流失预测、设备故障时间预测的读者。如果你现在还在“无脑KM + 无脑Cox”，这篇文章值得直接收藏。

1. 三板斧核心能力速览

三板斧不是互相替代的关系。

时依Cox是Cox模型的扩展，解决“比例风险假定不成立”的统计诊断问题；竞争风险模型解决的是“结局定义和删失机制”的问题；机器学习生存分析则偏向预测建模，适合做风险分层和筛选变量。实际项目里，通常是先用KM和基础Cox做探索，再根据数据诊断结果选择进阶方案。

2. 适用场景与使用边界

2.1 适合谁用

• 临床科研人员：肿瘤随访数据、术后复发、心脑血管事件、死亡与疾病进展并存。
• 生信分析人群：组学特征筛选、预后模型构建、风险评分。
• 金融风控：贷款逾期、客户流失、设备故障时间预测。
• 可靠性工程：产品寿命、维修时间、故障间隔。

2.2 能解决什么问题

基础Cox模型假设协变量对风险的影响是恒定的，也就是PH假定。但真实数据里，药物早期效果好、晚期效果减弱，或者某种生物标志物的预测能力只在确诊后前两年有意义。这时Cox模型的HR是一个平均效应，会掩盖时间趋势。

竞争风险场景更常见。比如研究“肿瘤复发”，但部分患者先死亡了。如果不处理死亡，只把死亡当作删失，累计复发率会被高估，因为真实世界里有相当一部分人根本没机会复发。

机器学习的引入，则是为了处理Cox模型不擅长的非线性、交互作用和高维变量。特别是基因表达谱、影像组学这类特征数量远大于样本量的数据，传统回归很难稳定估计。

2.3 不适合什么场景

• 没有事件时间、没有删失结构的普通二分类数据，不要硬套生存分析。
• 样本量极小、事件数极少时，机器学习生存模型容易过拟合，优先用朴素Cox或带惩罚的Cox。
• 因果推断问题，比如“这个药是否真的有效”，需要方案设计、随机化和因果推断框架，不是单纯跑一个回归就能回答。
• 不使用代码、只依赖SPSS菜单操作的场景，时依协变量和随机生存森林的实现会比较麻烦。

2.4 使用边界与合规提醒

涉及真实患者数据时，必须确认伦理审批、数据脱敏和隐私保护措施到位。临床数据集通常不能直接公开，更不能上传到未授权的第三方平台。生存分析模型用于发表论文时，要保留数据字典、分析脚本和版本记录。涉及商业数据、用户行为数据，也同样要注意授权范围，不能拿未授权数据做模型训练。

3. 环境准备与数据格式

3.1 R环境准备

R语言是生存分析最成熟的生态。核心包建议一次装好：

install.packages(c( "survival", "cmprsk", "randomForestSRC", "riskRegression", "timeROC", "ggplot2" ))

survival提供KM、Cox、cox.zph等基础能力，cmprsk提供竞争风险的CIF和Fine-Gray回归，randomForestSRC提供随机生存森林，riskRegression提供时间依赖AUC和Brier Score。

3.2 Python环境准备

Python生态用scikit-survival和lifelines比较顺手。

pip install scikit-survival lifelines pycox

如果要用Fine-Gray模型，Python原生支持不如R完整，实际项目里可以R做统计建模，Python做生产线部署。

3.3 数据格式要求

生存分析数据至少有三类信息：

• 生存时间：必须是正数，单位统一，比如天、月、年。
• 事件状态：普通生存分析用0/1，0表示删失，1表示目标事件。
• 协变量：可以是连续变量、分类变量或高维特征矩阵。

竞争风险模型的事件状态要多一个编码：

用R读入数据后，可以先检查一下结构：

str(dat) summary(dat$time) table(dat$status)

数据清洗时重点看时间是否有0或负数，状态编码是不是只有0、1，分类变量是否被正确设为factor。

4. 第一板斧：时依协变量Cox模型

4.1 什么时候该用

基础Cox回归有一个强假定：任意协变量的风险比不随时间变化。这个假定不满足时，会出现几种典型信号：

• KM曲线交叉，Log-rank检验P值不显著但趋势明显。
• Schoenfeld残差图显示残差随时间有趋势。
• 药物在随访早期降低风险，晚期反而无效。

一旦出现这些信号，不一定要放弃Cox模型。用“时依协变量Cox模型”可以显式建模“效应随时间变化”的过程。

4.2 先跑PH假定检验

先拟合一个基础Cox模型，再用cox.zph检验PH假定：

library(survival) fit_cox <- coxph( Surv(time, status) ~ age + sex + treatment, data = dat ) zph <- cox.zph(fit_cox) print(zph) # 可视化Schoenfeld残差 plot(zph)

输出结果里会给出每个协变量的P值和全局P值。如果某个变量的P < 0.05，或者残差曲线明显不水平，说明该变量的效应可能随时间变化。

注意：cox.zph是诊断工具，P值受样本量影响很大。大样本里轻微偏离也会得到小P值，所以还要结合残差图判断偏离程度。

4.3 用时变系数函数构建模型

survival包里的tt()函数可以构造时间交互项：

fit_tt <- coxph( Surv(time, status) ~ age + sex + treatment + tt(treatment), data = dat, tt = function(x, t, ...) x * log(t + 1) ) summary(fit_tt)

这个模型里，treatment是主效应，tt(treatment)表示治疗效应与log(t + 1)的交互。如果交互项显著，说明治疗的风险比随时间变化，不能再用单一HR表达。

常用的时间变换有：

• x * log(t + 1)
• x * t
• x * sqrt(t)

选择哪个函数，可以用AIC或看残差图来比较。

4.4 结果怎么解释

假设模型结果：

• treatment的系数为负，说明治疗早期有保护作用。
• treatment:log(t+1)的系数为正，说明保护作用随时间减弱。

实际报告时，不能只报一个HR。更合理的做法是给出HR随时间变化的曲线，或者列出几个代表时间点的HR，比如6个月、12个月、24个月时的HR。

用timereg包可以估计时变系数并画置信带，适合做稳健展示。

4.5 替代方案：时间分段Cox

如果不想用复杂的tt()函数，可以用时间分段模型。假设随访24个月，按12个月切分成两段：

dat_split <- survSplit( Surv(time, status) ~ ., data = dat, cut = c(12), episode = "tgroup" ) fit_seg <- coxph( Surv(tstart, time, status) ~ age + sex + treatment:strata(tgroup), data = dat_split ) summary(fit_seg)

这样会得到两段的HR：前12个月的HR，以及12个月后的HR。切分点需要根据临床意义和随访分布确定，不要盲目用中位数。

5. 第二板斧：竞争风险模型

5.1 什么是竞争风险

简单理解就是：一个人还没等到目标事件发生，先发生了其他事件，导致目标事件无法继续被观测。

典型例子：

• 研究“血液肿瘤复发”，但患者先发生非复发死亡。
• 研究“心梗复发”，但患者先死于车祸。
• 研究“首次骨折”，但患者先死亡。

如果把这些竞争事件直接当删失处理，普通KM估计的累积发生率会偏高，因为删失在KM里被默认成“以后还可能发生事件”，而竞争事件发生后，这个可能性已经不存在了。

5.2 先画CIF曲线

竞争风险下，描述结局的曲线不用KM，而是累计发生函数CIF。

用R的cmprsk包：

library(cmprsk) dat$status <- ifelse(dat$event == "death", 1, ifelse(dat$event == "relapse", 2, 0)) cif <- cuminc( ftime = dat$time, fstatus = dat$status, group = dat$treatment ) print(cif) plot(cif)

fstatus必须是一个数字向量，0是删失，1是目标事件，2是竞争事件。group可以是分组变量。cuminc会自动输出Gray检验的P值。

看到CIF曲线后，如果两条曲线明显分开，说明不同组的目标事件累积发生率有差异。

5.3 Fine-Gray亚分布风险模型

CIF曲线只能做分组比较，回归分析要用Fine-Gray模型，得到的是“亚分布风险比”sHR。

fg_fit <- crr( ftime = dat$time, fstatus = dat$status, cov1 = data.frame( age = dat$age, sex = dat$sex, treatment = dat$treatment ), failcode = 1, cencode = 0 ) summary(fg_fit)

解释注意：

• 普通Cox的HR，解释为“暴露组相对非暴露组，事件发生风险增加多少倍”。
• Fine-Gray的sHR，解释为“暴露组相对非暴露组，在存活动力学中目标事件累积发生率的风险比”。

sHR的数值不能直接当作HR来写。两个模型回答的问题不同。

5.4 cause-specific Cox 与 Fine-Gray 怎么选

竞争风险模型有两个流派：

• cause-specific Cox：将竞争事件当作删失处理，估计“原因别风险”。
• Fine-Gray模型：直接建模目标事件的累积发生率，考虑竞争事件对风险集的影响。

粗略选择原则：

• 想探索某个因素的病因学机制，用cause-specific Cox。
• 想预测某个患者个体发生目标事件的概率，用Fine-Gray。
• 如果研究目标是“复发概率”，临床决策更关心累积发生率，Fine-Gray更直接。

稳妥做法是两种都跑，放在敏感性分析里比较结论是否一致。

6. 第三板斧：机器学习生存分析

6.1 为什么还需要机器学习

Cox模型有两个限制：

• 假设协变量对对数风险是线性加成。
• 依赖PH假定。

实际数据里，年龄和风险可能是U型关系，基因之间可能存在复杂交互。用Cox硬建模，要么变量被丢弃，要么预测效果差。

机器学习生存模型，比如随机生存森林、梯度提升生存模型、DeepSurv，可以自动处理非线性、交互和高维特征。代价是可解释性变差，且需要更严格的验证。

6.2 R语言随机生存森林

randomForestSRC是R里比较成熟的实现：

library(randomForestSRC) rsf_fit <- rfsrc( Surv(time, status) ~ ., data = dat, ntree = 500, nodesize = 20, seed = 42 ) print(rsf_fit) # 查看OOB误差 rsf_fit$err.rate[rsf_fit$ntree]

rfsrc输出里自带OOB的C-index、Brier Score和变量重要性。ntree建议先500，nodesize设为10-20。如果变量数特别多，可以让mtry保持默认，再根据OOB误差微调。

预测新样本的生存概率：

pred_rsf <- predict(rsf_fit, newdata = newdat) plot(pred_rsf$time.interest, pred_rsf$survival[1, ], type = "l")

6.3 Python实现随机生存森林

Python推荐scikit-survival：

import pandas as pd from sksurv.ensemble import RandomSurvivalForest from sksurv.util import Surv from sksurv.metrics import concordance_index_censored # 构造结构化生存标签 y = Surv.from_dataframe("status", "time", dat) # X 是特征表，不要包含time和status X = dat.drop(columns=["time", "status"]) rsf = RandomSurvivalForest( n_estimators=500, min_samples_leaf=20, random_state=42, n_jobs=-1 ) rsf.fit(X, y) # 训练集C-index，正确评估需要交叉验证 c_index = concordance_index_censored( y["status"], y["time"], rsf.predict(X) ) print(c_index)

注意：Surv.from_dataframe的status列必须是布尔值或0/1，time必须是浮点数。预测值越大代表风险越高，不是生存概率。

6.4 模型评估不能只看C-index

生存模型不能像普通分类模型那样只算AUC。需要用时间依赖指标：

• C-index：全局区分度，适合快速比较。
• 时间依赖AUC：反映某个时间点的区分能力。
• Brier Score：预测概率与真实结局的校准程度。
• 校准曲线：预测生存概率和实际生存概率是否一致。

R语言的riskRegression可以同时算多个指标：

library(riskRegression) Score( list("RSF" = rsf_fit, "Cox" = fit_cox), formula = Surv(time, status) ~ 1, data = dat, metrics = "auc", times = c(12, 24) )

实际使用前需要确认Score对自定义模型对象的兼容性，具体以包文档为准。

6.5 可解释性

随机生存森林可以输出变量重要性VIMP，画偏依赖图：

plot.variable(rsf_fit, xvar.names = "age", partial = TRUE)

偏依赖图能显示某个变量取值变化对预测生存概率的影响趋势，对发现非线性关系很有帮助。

在Python里，也可以用置换重要性观察变量对C-index的影响。需要说明的是，机器学习模型的“特征重要性”不等于临床意义上的因果效应，不能直接写成“该变量是独立预后因素”。

7. 三板斧联用实战流程

7.1 决策流程参考

不要跳过基础分析直接上进阶模型。推荐按下面的思路走：

1. 清洗数据，确定结局、时间、竞争事件定义。
2. 做描述统计，计算事件率、删失比例。
3. 画KM曲线做分层探索，做Log-rank检验。
4. 拟合基础Cox，做cox.zphPH诊断。
5. 如果PH不满足，尝试时依协变量Cox或时间分段模型。
6. 如果存在竞争事件，画CIF，跑Fine-Gray回归。
7. 如果特征多、关系复杂，建模目标偏预测，则上随机生存森林。
8. 用Bootstrap或交叉验证评估区分度和校准度。
9. 报告时明确说明模型选择依据，不隐藏基础模型的不足。

7.2 批量建模：多个终点循环

实际项目里经常要同时分析多个结局，比如“疾病进展”“死亡”“复合终点”。可以用循环批量建模，把结果汇总成表格：

endpoints <- c("status_progression", "status_death", "status_composite") result_list <- list() for (ep in endpoints) { dat$ep <- dat[[ep]] fit <- coxph( Surv(time, ep) ~ age + sex + treatment, data = dat ) result_list[[ep]] <- summary(fit)$coefficients } result_list

批量建模时，建议每个模型都输出样本量、事件数、删失数。事件数太少的结果要谨慎解读。

7.3 接口服务与自动化

生存分析本身通常不是在线API服务。如果需要把训练好的模型集成到业务系统，更常见的做法是：

• R里的模型用rds保存，配合plumber封装成HTTP接口。
• Python的sksurv模型用joblib保存，配合FastAPI提供预测服务。
• 定期用脚本重训模型，输出指标到Excel或数据库中。

接口健康检查可以参考常见AI服务的做法：固定端口、限制访问IP、记录日志。

8. 计算效率与性能观察

生存分析不是典型的深度学习任务，但仍然有资源消耗问题。

8.1 不同模型的计算开销差异

• 基础Cox回归：一般秒级完成，几百样本甚至不到1秒。
• 时依协变量Cox：取决于是否时间分段，分段后数据行数变多，计算量增加。
• Fine-Gray回归：数据量在万级以内通常很快。
• 随机生存森林：几百棵树在小样本上很快，但高维组学数据会明显变慢。
• DeepSurv等神经网络生存模型：需要GPU或者较大内存，调参成本更高。

具体时间没有统一数字，建议运行前先在小数据子集上做计时。

8.2 怎么观察资源占用

R里可以用：

system.time({ fit <- rfsrc(Surv(time, status) ~ ., data = dat, ntree = 500) })

Python里用：

import time start = time.perf_counter() rsf.fit(X, y) end = time.perf_counter() print(end - start)

训练前可以打开系统任务管理器或top命令观察CPU和内存占用。如果出现内存不足，优先减小ntree、增大nodesize、减少特征数。

8.3 降低计算成本的实用方法

• 随机生存森林先用ntree = 200跑通，再根据误差曲线决定是否增加。
• 高维数据先用单变量筛选或Lasso-Cox降维，再进随机森林。
• 交叉验证时，尽量用重复5次的分层抽样，而不是1000次Bootstrap。
• Python设置n_jobs = -1利用多核。
• 时间依赖AUC计算比较耗时，先选2到3个有代表性的时间点。

8.4 性能指标和业务指标结合

统计指标不能只看C-index。C-index 0.7以上就算中等区分度，但在临床里不一定有决策价值。建议同时报告：

• 删失比例和事件数。
• 时间依赖AUC。
• 校准曲线。
• 不同风险分层的实际事件率。

9. 常见问题与排查方法

10. 最佳实践与使用建议

10.1 先跑通最小可运行流程

第一次分析不要直接上全部数据。取一个子集，跑通KM、COX、cox.zph、CIF、Fine-Gray、随机生存森林，确认输出格式符合预期后再上全量数据。

10.2 数据、代码、输出分目录管理

推荐这样的项目结构：

project/ ├── data/ │ ├── raw/ │ └── processed/ ├── code/ ├── output/ │ ├── figures/ │ └── tables/ └── report/

分析脚本要有编号，比如01_data_clean.R、02_km_cox.R、03_competing_risk.R、04_rsf.R。所有输出文件保存日期和版本。

10.3 模型报告要完整

写论文时至少报告以下内容：

• 随访时间的中位数和范围。
• 各结局的事件数和删失数。
• 模型选择依据，特别是PH检验结果。
• 竞争风险的定义和编码方式。
• C-index、时间依赖AUC、Brier Score及置信区间。
• 校准曲线或校准表。

10.4 合规使用数据和模型

真实患者数据必须脱敏，不能上传到无授权的在线工具。分析脚本和模型文件尽量不要包含原始身份证号、姓名等敏感信息。涉及商业决策时，模型输出需要复核，不能只依赖单一指标。

11. 总结与下一步

这三板斧最值得先试的是第一板斧和第二板斧。它们不改变你的整体分析框架，只在你已有的KM+Cox基础上加一步诊断和一步修正。先从cox.zph和cuminc入手，跑完你会很清楚地知道自己的数据到底适合哪种模型。

第三板斧最适合预测建模场景，不要在样本量只有几十例时强行使用。如果你手里已经有几百个样本、几十个以上特征，再考虑随机生存森林或DeepSurv。

最容易踩的坑是状态编码。普通生存分析只有0和1，竞争风险模型多一个2，很多报错和结果偏差都来自这里。建议在分析开始前先写一行table(survival_data$status)，确认编码没有问题再往下走。

下一步可以扩展的内容包括：时间依赖ROC曲线、校准曲线、外部验证队列、DeepSurv多模态数据、以及把最佳模型封装成REST API供业务系统调用。