R语言：生存分析竞争风险模型，从理论到实战总结

临床研究中结局事件可能有多个, 例如结局B将阻止感兴趣事件A的出现或影响A发生的概率，B事件形成"竞争”关系,B为A竞争风险事件。常见场景是研究治疗方案A和复发的关系，如果患者在死亡就观察不到复发，“死亡”和“复发”存在竞争。传统生存分析回归前提假设是删失时间与失效时间是独立的且删失对回归系数β是无信息的，即结局不存在竞争风险。

临床生存数据常常伴有多个结局，结局间若存在竞争关联关系，简单地将其他原因别死亡与截止依然存活视为相同删失是错误的，违背了删失独立假设，此时若使用传统生存分析会高估累积发生率，所以，需要引入竞争风险cox分析来正确估计累积发生率。

竞争风险理论

竞争风险模型是多状态模型中的一种标准状态结构，示意如图1。竞争风险模型适用于多个终点生存数据，关心终点A与不关心终点B非相互独立且还存在竞争关系，A发生导致B不会发生，例如慢性肾患者死亡与透析、心梗导致的死亡与其他死因，生殖细胞癌死亡与继发恶性肿瘤，先心病术后死亡与随访终点肺静脉梗阻存在竞争风险。总的来看都是，死亡与关心终点存在竞争风险关系。竞争风险分析单因素分析常为估计终点A发生率，多因素分析常为估计预后影响因素及效应值。

图1 竞争风险模型示意图。λ为竞争风险函数

1.1 单因素KM > CIF

KM法无论A、B、合计事件发生率均高于累积发生函数(cumulative incidence function, CIF)法，如图2。KM对应生存曲线，差异性检验主要对应log-rank 检验。累积发生函数F(t)，F (t)=1-S(t)=Pr(T<t)，其中生存函数S(t)，又称生存率，即KM。CIF又称累积发生率，CIF对应曲线为Nelson-Aalen累积风险曲线，差异性检验主要对应Gray’s检验。当存在竞争风险时，结局不再是二分类生存、死亡，此时CIF≠F(t)，而CIF意为各自的关心事件累积发生函数、竞争事件累积发生函数。CIF总= CIFa + CIFb。

1.2 多因素净危险率λ (t) ≠ 竞争粗危险率λk(t)

若无竞争风险，即满足假设“删失独立”的前提下，Cox比例风险模型为 log⁡(λ(t))=log⁡(λ_0 (t))+Χβ，λ(t)称净危险率，λ0(t)为基线风险函数，即协变量向量为0时的风险函数，可写成λ(t)=λ_0 (t)exp(Χβ)。存在竞争风险，此时“删失独立”条件不满足，存在两种模型：原因别风险函数(cause-specific hazard function, CS)、部分分布风险函数(subdistribution hazard function, SD )。

CS公式为λ_k^cs (t)=(lim)┬(∆t→0)⁡〖(P(t≤T<t+∆t,D=k|T≥t))/∆t〗，CS表示t时刻未发生任何事件的观察个体发生第k类事件的瞬时概率强度。

SD公式为λ_k^sd (t)=(lim)┬(∆t→0)⁡〖(P(t≤T<t+∆t,D=k|T>t∪(T<t∩K≠k)))/∆t〗，SD表示t时刻未发生第k类事件的观察个体发生第k类事件的瞬时概率强度，又称CIF回归模型、Fine-Gray模型。

CS适合回答流行病病因学问题，解决了比例风险模型中不能同时较准确地考虑多个终点事件,适合病因学研究，CS回归系数反映了协变量对无事件风险集对象中主要终点A发生率增加的相对作用，缺点是结果的解释不是很直观；SD适合临床问题、建立临床预测模型及风险评分，仅对终点A绝对发生率感兴趣。SD获取Nelson-Aalen累积风险曲线，χ2，Gray检验参数。

图2 累积发生函数CIF与Kaplan-Meier估计发生率示意图

1.3 CS、SD模型示例

CS与SD两模型差异风险集是否包含竞争事件B，SD分母仅剔除之前的A不剔除B，仅对终点A绝对发生率感兴趣。例如，时间点0风险集为10个观察对象，时间点2风险集CS模型为8人=10-1A-1B，SD模型为9人=10-1A，事件A风险分别为0/9、0/8，时间点3风险集CS模型为6人=8-2B，SD模型为9人=9-0A，事件A风险分别为2/9、2/6。

图3  风险集起始为10个观察对象，关心事件A为正方形，竞争风险事件B为三角形。

1.4 有无竞争风险模型区别

二、顶刊案例

2.1 中文范文

小编2017年于中华流行病杂志利用全国多中心10年连续入组先心病数据库，展示了术后死亡与PVO(关心结局)存在竞争风险关系，计算了KM、CIF率的差异；CS、SD模型的HR差异。发生率KM > CIF。SD对应发生率，CS对应原因别风险。例如，SD模型提示某var因素增加/减少术后PVO相应的发生率%，CS模型提示var因素增加/减少术后PVO原因别风险%。

Table3 术后PVO多水平cox分析、部分分布风险模型、原因别风险模型

2.2 SCI范文

三、R、SAS实战

本章节方法学与思维导图如下。

3.1 R代码

R中支持竞争风险的包主要有cmprsk、rms、QHScrnomo、mstate、tidycmprsk、riskRegression、tidycmprs包等等。cmprsk 包作者网站https://luca-scr.github.io/software/，里面附带的数据目前无法下了。

本文案例演示来自cmprsk自带数据bmtcrr.rdata。

本文R代码主要来自文献3、文献7，需要代码请移步文献3、7获取。

3.2 SAS代码

/*cohort为数据名，time为生存时间，结局y多类别间存在竞争风险，赋值编码0=生存、1=主要关心结局A、2=竞争风险事件B*/
/*cov1-covn为协变量名*/
ods graphics on; ods html dpi=300;
%cif(data=cohort,out=cif1_data,time=time,status=y,event=1,censored=0); /*不分组A总的CIF*/
%cif(data=cohort,out=cif1_data,time=time,status=y,event=1,censored=0,group=treat); /*手术组treat组间gray's比较A的CIF*/
%cif(data=cohort,out=cif2_data,time=time,status=y,event=2,censored=0,group=treat); /*手术组treat组间gray's比较B的CIF*/
/*A校正CIF，SD模型*/
proc phreg data=cohort plots(overlay=stratum)=cif; 
model time*y(0)=cov1~covn /eventcode = 1 rl; run;
/*A的CS模型，y(0,2)表示删失事件，zph为Schoenfeld残差图、assess为鞅式残差图，均可验证PH比例风险假设*/
proc phreg data= cohort zph;
model time*y(0,2) = cov1~covn /ties=efron rl; assess ph/resample seed=1234; run;
proc phreg data=cohort plots(overlay=stratum)=cif; /*B校正CIF，SD模型*/
model time*y(0)= cov1~covn /eventcode = 2 rl; run;
proc phreg data= cohort zph; /*B的 CS模型*/
model time*y(0,1) = cov1~covn /ties=efron rl; assess ph/resample seed=1234; run;

四、发散思维

4.1 多状态模型

传统竞争风险模型虽然能调整其他结局事件对主要结局事件发生风险的竞争影响，但无法直接比较个体发生不同结局事件之间的风险差异。多状态竞争风险模型为此提供了相对 合适的解决方案，后续会讲解。

4.2 竞争注意问题

竞争事件(死亡)比例>10%时候，若采用传统km/cox方法可能造成严重bias；当竞争比例<10%，依然可能出现假阳性或假阴性。CS回答上游的流行病病因学问题，SD回答下游的绝对发生率问题，SD模型越来越多的被用来临床预测模型和风险评分。虽然SD模型越来越多，但是SD模型并非一定比Cox模型更优，因为解释起来会费劲点，最佳策略是把SD当做sensitivity analysis与经典cox同时展示结果的稳定可靠。

4.3 模型评价

R中riskRegression包可以对基于竞争风险模型构建的预测模型进行进一步评价，比如计算C-index，绘制校准曲线等。

4.4 PH假设

无论经典cox、CS/SD竞争风险都需要符合PH假设，当存在竞争风险又不符合PH假设的时候怎么办呢？方法学https://academic.oup.com/aje/article/170/2/244/111339，介绍了不满足ph假设采用Parametric Mixture Model，但是小编查阅资料几乎很少后续PMM文章。更推荐采用时依系数法、分段模型法和分层分析法，文献12也提供了很好的解决方案。

4.5 CS/SD可能不一致

SD与CS结论不一致的时候，例如Peter C研究中CS模型提示肿瘤仅增加非心血管病原因别死亡风险，SD模型提示肿瘤降低心血管死亡率增加非心血管病死亡率，换句话说肿瘤降低心血管病死亡率是由于对应的非心血管原因别风险增加间接导致的。

五、小结

“

1.存在竞争风险时候，关心结局发生率CIF < KM。
2.竞争事件率>10%，推荐CIF、SD模型直接反应出绝对率的改变；竞争事件率<10%，最佳策略是把SD当做sensitivity analysis与经典cox同时展示结果的稳定可靠。
3.CS、SD也都需要ph假设，不满足ph需要分层cox等其他方法。
4.单因素时推荐QHScrnomo、tidycmprsk包绘制CIF曲线及Fine-Gray检验，多因素时推荐riskRegression包制作table2 (sHR)+ ROC曲线，QHScrnomo计算cindex、calibration、nomogram、tidycmprs包展示森林图，需要注意哑变量ref问题。”

参考文献：

1.聂志强等.临床生存数据新视角：竞争风险模型[J].中华流行病学杂志, 2017, 38(8):5.DOI:10.3760/cma.j.issn.0254-6450.2017.08.026.
2. 临床研究与医学统计，Fine-Gray检验与竞争风险模型
3. 医学和生信笔记，Fine-Gray检验、竞争风险模型、列线图绘制
4. Scrucca L. (2010) Regression modeling of competing risk using R: an in depth guide for clinicians. Bone Marrow Transplantation, 45, 1388–1395.
5. Zhang Z; written on behalf of AME Big-Data Clinical Trial Collaborative Group. Overview of model validation for survival regression model with competing risks using melanoma study data. Ann Transl Med 2018;6(16):325. doi: 10.21037/atm.2018.07.38
6. Beyersmann J. Simulating competing risks data in survival analysis. Stat Med 2009 Jan 5;28(1):956-971
7. 一只勤奋的科研喵，竞争风险分析: 单变量和多变量回归
8. Austin PC. Introduction to the Analysis of Survival Data in the Presence of Competing Risks. Circulation. 2016 Feb 9;133(6):601-9. doi: 10.1161/CIRCULATIONAHA.115.017719. PMID: 26858290.
9. Lau B. Competing risk regression models for epidemiologic data. Am J Epidemiol. 2009 Jul 15;170(2):244-56. doi: 10.1093/aje/kwp107. Epub 2009 Jun 3. PMID: 19494242.
10. https://zhuanlan.zhihu.com/p/539801600?utm_id=0
11. 一只勤奋的科研喵，预测模型 | 13. 竞争风险Nomogram (基于mstate包) 12. 统计咨询，竞争风险模型与半式竞争风险模型介绍及R实现
13. 李卫芹, 杨雷, 王胜锋, 等.  生存数据的多状态竞争风险模型分析方法探索 [J] . 中华预防医学杂志, 2021, 55(12) : 1524-1529. DOI: 10.3760/cma.j.cn112150-20211103-01019.