科研猫小课堂：敲黑板！竞争风险模型应该如何分析？

作者：科研猫 | 西红柿

责编：科研猫 | 馋猫

背景

在观察某一事件是否发生时，如果该事件受到其他事件的阻碍，在这种所谓的竞争风险研究中可能会有多个结果事件，有些结果会阻止感兴趣的事件发生或影响其发生的概率。所有形成竞争关系的结果事件，互为竞争的风险事件。

例如，有研究人员收集了本市2007年诊断为轻度认知障碍(MCI)的518例老年患者的临床资料，包括基本人口学特征、生活方式、体检、疾病信息等，并在2010-2013年完成6次随访调查，主要结果事件是这些老年患者会不会发生阿尔茨海默病(AD)。随访期内共发生AD 78例，这其中包括迁移28例，失访31例，死亡25例。影响MCI向AD过渡的因素有哪些？

在这种情况下，如果MCI患者在没有发生AD的观察期内死于癌症、心血管疾病、车祸等原因，则不能观察到AD的发病，也就是死亡终点与AD的发生产生了“竞争”。根据传统生存分析方法，个体死亡发生在AD之前，结局是死亡，而不是AD，被认为是删失数据，可能会导致偏差。对于死亡率较高的老年人群，当存在竞争风险事件时，传统生存分析方法(Kaplan-Meier方法、logrank检验、Cox比例风险回归模型)会高估感兴趣疾病的风险，从而导致竞争风险偏差。一些研究发现，大约46%的文献可能存在这样的偏差。

在这种情况下，竞争风险模型是适用的。所谓竞争风险模型是对生存数据的多个潜在结果进行处理的一种分析方法。早在1999年，Fine和Gray就提出了部分分布的半参数比例风险模型，常用的终点指标是累积发生率函数(CIF)。在这种情况下，AD前死亡可视为AD的竞争风险事件，并采用竞争风险模型进行统计分析。竞争风险的单变量分析通常被用来估计感兴趣的终点事件的发生率，多变量分析经常被用来探索预后因素和影响大小。

案例分析

此案例数据是从

http://www.stat.unipg.it/luca/R/

上下载的，研究人员计划比较骨髓移植和血液移植治疗白血病的疗效。终点事件定义为“复发”，一些患者在移植后不幸因为不良反应出现死亡，所以会因移植相关的死亡无法观察到患者末期的“复发”。换句话说，“移植相关死亡”和“复发”是竞争性风险事件。因此，本文采用竞争风险模型进行统计分析。

首先，从当前工作路径导入数据文件‘bmtcrr.csv’。可以找我们的工作人员领取。

1library(foreign) 
2bmt <-read.csv(‘bmtcrr.csv’) 
3str(bmt)
4library(cmprsk)
5bmt$D <- as.factor(bmt$D) 
6attach(bmt)

这是一个数据表结构的数据，有7个变量，总共177个观测值。

$ Sex : 因子变量，2个水平：“F”,“M”。

$ D : 因子变量，2个水平：“ALL(急性淋巴细胞白血病)”,“AML(急性髓系细胞白血病)”。

$ Phase : 因子变量，4个水平：“CR1”, “CR2”, “CR3”, “Relapse”。

$ Age : 年龄。

$ Status: 结局，0=删失，1=复发，2=竞争风险事件。

$ Source: 因子变量，2个水平：“BM+PB(骨髓移植+血液移植)”,“PB(血液移植)”。

$ ftime : 时间。

Fine-Gray检验（单因素分析）

与比较两组生存结果数据的对数秩和检验类似，考虑竞争风险事件也可以进行单因素分析。接下来，我们可以使用cuminc()函数进行单变量Fine-Gray检验

1fit1 <- cuminc(ftime,Status,D) 
2fit1

结果解释：“1”表示定义的终点，“2”表示竞争风险事件。第一行统计量=2.8623325，P值=0.09067592，说明控制竞争风险事件后(即第二行统计量和P值)，“ALL”和“AML”的累积复发风险无统计学差异，P=0.09067592。

$ est：表示每个时间点“ALL”和“AML”组的估计累计重复率和累计竞争风险事件发生率(定义的端点和竞争风险事件分别以“1”和“2”区分，与上面第1行和第2行的结果一致)。

$ var：表示每个时间点“ALL”和“AML”组(定义的终点和竞争风险事件分别由“1”和“2”标识，与上面第1行和第2行的结果一致)的估计累积重复率和累积竞争风险事件发生率的差异。下面，我们绘制了累积复发率和累积竞争风险事件发生率的生存曲线，以直观地表示上述数字化结果(图31)。

1plot(fit1,xlab = ‘Month’, ylab = ‘CIF’,lwd=2,lty=1, 
2     col = c(‘red’,’blue’,’black’,’forestgreen’))

图解：纵轴代表累计发生率（CIF），横轴代表时间。我们查看了对应于ALL1的红色曲线和对应于AML1的蓝色曲线(1表示定义的终点，2表示竞争风险事件)。从图中可以看出，ALL组的复发风险均高于AML组，但差异无统计学意义，P=0.09067592。同样，如果在AML2对应的草绿色曲线下方观察ALL2对应的黑色曲线，可以得出结论：ALL组的竞争风险事件发生率均低于AML组，也未达到统计学意义，P=0.50322531。从曲线可以看出，曲线在前20个月是“纠缠”的，所以没有得到统计学上有意义的结果。简单来说，这个数字可以用一句话来概括：在控制竞争风险事件后，“ALL”和“AML”的累积复发风险没有统计学差异，P=0.09067592。

Fine-Gray检验（多因素分析）

以下是考虑竞争风险事件的生存数据的多变量分析。在cmprsk包中，crr()函数可以方便地进行多因素分析。函数用法如下：

1crr(ftime, fstatus, cov1, cov2, 
2    tf, cengroup, failcode =1, 
3    cencode =0, subset, 
4    na.action = na.omit, 
5    gtol =1e-06, maxiter =10, 
6    init, variance = TRUE)

您可以参考crr()函数帮助文档了解各个参数的详细说明。这里应该注意的是，函数必须指定时间变量和结果变量，然后传入协变量矩阵或数据框。首先，定义进入模型的协变量，并将其定义为数据框。

1cov <- data.frame(age = bmt$Age, 
2                 sex_F = ifelse(bmt$Sex==‘F’,1,0), 
3                 dis_AML = ifelse(bmt$D==‘AML’,1,0), 
4                 phase_cr1 = ifelse(bmt$Phase==‘CR1’,1,0), 
5                 phase_cr2 = ifelse(bmt$Phase==‘CR2’,1,0), 
6                 phase_cr3 = ifelse(bmt$Phase==‘CR3’,1,0), 
7                 source_PB = ifelse(bmt$Source==‘PB’,1,0)) 
8cov

构建多元竞争风险模型。在这里，需要指定failcode =1， cencode =0，分别表示：终点事件赋值“1”和截尾赋值“0”，并且其他竞争性风险事件被赋值“2”。

1fit2 <- crr(bmt$ftime, bmt$Status, cov , failcode=1, cencode=0) 
2summary(fit2)

结果解释：在控制竞争性风险分布事件后，phase 变量，即疾病阶段，是患者复发的独立危险因素。以“复发”患者为参照，CR1、CR2、CR3期患者与“复发”患者的累计复发率、危险度比、95%CI分别为0.332(0.159、0.695)、0.361(0.180、0.724)、0.481(0.155、1.490)，对应的P值分别为0.0034、0.0041、0.2000。

小结

本部分详细介绍了使用R的cmprsk包的Fine-Gray检验和竞争风险模型。笔者认为读者在具体应用过程中应该注意两点：

第一，可选择性的使用Fine-Gray检验和竞争风险模型，如果终点事件存在竞争风险事件，并且很可能影响结论，那么使用这个模型是合适的，这个模型不一定比Cox模型更好，这两个模型应该是相辅相成的；

第二，竞争风险事件的考虑也受到竞争风险模型的限制。目前，仅将Cox模型的二元端点扩展为三重分类，即结果事件、删失事件和竞争风险事件。即便如此，也很难解读结果。因此，读者在选择统计方法时，应该进行更充分的评估和实验。

科研猫原创出品，任何形式转载，均需获授权

若您是读者，欢迎分享，无需授权

易法通律师事务所提供版权法律支持