R笔记：生存分析之竞争风险模型[概念与实操]

转自个人微信公众号【Memo_Cleon】的统计学习笔记：R笔记：生存分析之竞争风险模型[概念与实操]。

在《生存分析之Kaplan-Meier曲线绘制与比较》文末，我们提到了竞争风险。由于出现了其他结局事件而阻碍或者改变了感兴趣的目标结局事件的发生（概率），相对目标结局事件而言，这里的其他结局事件就是竞争风险事件。实际上除了以全因死亡为研究结局，以死因、复发等为结局的研究都广泛存在着竞争风险。比如以心血管死亡为研究结局，其他原因的死亡（如车祸意外、非心血管疾病死亡）发生时，心血管死亡就不会再出现了，但如果将其他原因的死亡按删失（censoring）处理，就意味着只要随访时间足够长，这些其他原因导致的死亡对象仍会再发生心血管疾病死亡，这显然与事实不符。同样以疾病复发为结局事件的研究中，研究对象的死亡同样导致该对象不会再发生疾病的复发。也就是说当失效时间与删失时间不独立时。将竞争风险作为删失处理是不恰当的。进一步你可能会问，删失是不是目标结局事件的竞争事件呢？当受试者失去随访时（发生删失），这些受试者仍然有发生目标结局事件的风险，只是研究者无法记录发生的准确时间而已，而竞争风险的发生是阻碍了目标结局事件的发生（如死亡发生后，疾病不会再出现复发）。换句话说，删失只是阻止了你观察到目标结局事件，而竞争事件则是阻止了目标结局事件的发生，所以删失并不是竞争事件，竞争事件的分析需要用一些新的方法。

竞争风险模型就是用于这种多终点的生存数据的分析。在竞争风险模型中，一个关键点是，一个受试者只能发生一个结局事件类型，这样该事件的发生就排除了后续再发生其他事件类型。

本次笔记内容：（1）竞争风险模型理论基础；（2）R操作示例：累积发生率单因素比较；（3）R操作示例：特定原因风险回归；（4）R操作示例：子分布风险回归/累积发生率回归。

【1】竞争风险模型理论基础

此部分有助于理解，不感兴趣的可只关注下面这几个主要总结性的结论：

• 不存在竞争风险时，累积发生率函数与生存函数、风险函数、累积风险函数可以通过简单的函数直接转换，对累积发生率的研究可以转变成结局事件生存函数、风险函数/累积风险函数的分析，expβ可以作为考察协变量对累积发生率的相对影响的度量。当存在竞争风险时，累积发生率不仅取决于目标结局事件的风险函数，而且与竞争事件的风险函数也有关系，此时累积发生率不能简单地转变为对风险函数和生存函数的研究；
• 特定原因风险集排除了发生竞争事件的观察人群，特定原因的HR可以反应暴露对目标结局事件的影响的关系，而竞争事件可以从风险集中移除观察人群来从侧面影响这种关系。子分布风险集的构建包含了没发生任何事件的观察人群和那些已经发生竞争事件的人群；
• 特定原因的风险比可以表示暴露的出现使得结局发生的风险是原来的expβ倍（或者比原来增加100(expβ-1)%)，但其值>1并不意味着暴露组的累积发生率大于非暴露组；子分布风险比表示暴露的出现使得结局的累积发生率增加（HRsd>1）或者减少（HRsd<1）。

竞争风险的中文学习教材很少，生存分析中的概念本就晦涩难懂，加上有些翻译的也不准确（个人以为是这样），很容易给人造成困扰。所以还是要认真的理解基本的概念。

基本概念：

在经典的生存分析中有几个非常重要的概念：生存时间T、风险集、失效函数F(t)、生存函数S(t)、风险函数h(t)。不论有没有竞争风险，这都是非常重要的概念，只是存在竞争风险时有些指标需要进行修正。

死亡概率密度函数f(t)是一个概率密度函数（Probability Density Function,PDF），表示所有观察对象在t时刻的瞬时死亡率（instantaneous rate）【注：此处的死亡代表的是失效事件，并不一定是真正意义上的死亡】。

累积失效概率：也称为累积失效分布函数、累积发生率函数（cumulative incidence function，CIF）、失效函数(failure function)，记为F(t)，F(t)=Pr(T≤t)，表示观察对象生存时间T不大于某时刻t的概率，或者说直到时刻t，结局事件发生的累积概率。这是一个概率累积分布函数（Cumulative Distribution Function,CDF）。

生存函数（survival function）：常用S(t)表示，是用来描述生存时间T分布的一个指标，定义为随机变量T越过时间点t的概率：S(t)=Pr(T>t)，表示时刻t仍未发生终点事件（如死亡）的概率，即t时刻仍生存的概率。t=0时，S(t)=1，随着时间的延长，S(t)取值逐渐减小。生存函数是一个累积函数，也称为生存累积分布函数或累积生存概率（cumulative survival probability），也简称为生存率（survival rate）。从定义上来说概率则更符合这种累积函数的定义；函数表达的是一种关系，生存函数可以理解为随时间变化的生存时间概率；而率一般说明某种现象发生的强度，常与时间有关系，将累积生存概率翻译为生存率个人觉得不利于理解。

风险函数（Hazard function）h(t)：定义为f(t)/S(t)，也称为风险率（hazard rate），到t时刻仍存活的观察对象在接下来瞬间（t时刻）发生结局事件（如死亡、复发等）的概率。h(t)是相对于t时刻仍存活的人，而并非所有人。当这个“瞬时”刚好等于1单位时，此时的风险函数就是t时刻仍然存活的个体在此后一个单位时间段内的死亡概率。风险函数也是风险的一种指标，虽然它表达的是瞬时死亡率，从定义上看但它既不是概率密度，也不是累积概率分布。风险函数的累积分布函数是累积风险函数H(t)。

当不存在竞争风险时，

由于此时CIF可以直接和生存函数或者累计风险函数直接关联，对CIF的研究可以转变为对目标结局事件生存函数或者累积风险函数的分析，只是我们往往关注的是生存函数S(t)=Pr(T>t and event type 1) 。

f(t）、F(t)、S(t)、h(t)、H(t)之间的数学关系如下：

单因素的生存分析常用到乘积极限法（Kaplan-Meier法），用1-S(t)来估计粗发生率；当考察协变量的影响时（多因素分析），常用到Cox风险比例回归模型。在Cox模型中，

由此衍生出风险比（HR）的概念：当固定其他因素，自变量每增加一个单位，风险比(HR)是expβ，或者：当固定其他因素，自变量每增加一个单位，事件发生风险增加100(expβ-1)%。Cox模型通过HR（即expβ）作为考察协变量对风险函数的相对影响，由于不存在竞争风险时累积风险函数与累积发生率直接相关，HR也可以考察协变量对CIF的相对影响。

---

当存在竞争事件时，在结局事件发生之前我们并不知道明天和意外哪个先到，是先发生目标结局还是先竞争事件。你可能会问“1年内发生心血管死亡的概率是多少”，而不是“1年前没有发生任何结局事件，接下来发生心血管死亡而不是非心血管死亡的概率”，因此我们更关注的是失效函数F(t)=P(T≤t and event type 1), 也就是大名鼎鼎的累积发生率函数（cumulative incidence function，CIF）。第k类结局事件的累积发生函数CIF_k(t)=Pr(T≤t,D=k),定义为“the probability of experiencing the kth event before time t and before the occurrence of a different type of event”，直译为“时间t以前且在发生其他类结局事件之前发生k类结局事件的概率”，或者这样更好理解：没有发生其他结局事件的受试者在时间t之前发生k类结局事件的概率。此时因为受竞争风险的影响，发生竞争事件的受试者不再发生目标结局事件，CIF_k与第k个结局事件和竞争事件都有关系，即CIF_k也不再等于1-S_k(t)。

通过特定原因的风险函数h_k(t)来计算CIF_k比较麻烦。对每个结局事件分别构建经典的Cox比例风险模型，将其中发生竞争结局事件的观察人群定义为删失，就可以获得不同原因的结局事件的风险函数h_k(t)：

表示无任何结局事件的观察人群中，k结局事件的瞬时发生率。f_k(t)是目标结局事件k发生的概率密度函数，d_tk和n_t分别表示t时刻结局事件k的发生人数和t时刻前的风险人数。这个特定原因的风险函数跟不考虑竞争风险的风险函数非常类似，说白了就是分别计算不同结局的风险函数而已。

这里的不同原因的结局事件的风险函数h_k(t)，cause-specifific hazard function，在众多学习材料上不知为何翻译为原因别风险函数，刚看见这个概念的时候一脸的懵逼，我一度还以为是笔误。其实个人觉得翻译为特定原因的风险函数或者特定结局的风险函数更贴合实际意义一些。

当存在竞争风险时，第k类结局的累积风险函数H_k(t)和生存函数S_k(t)定义为：

虽然S_k(t)可以估计，但是只有当竞争事件的时间分布和删失分布是独立的时候，才可以解释为边际生存分布（KM法生存率），而且在这种情况下，边际分布描述的是在竞争事件不发生的情况下的事件时间分布。

总体生存函数S(t)是在时间 t 没有发生任何结局事件的概率，用特定原因风险函数表达为：

第k类结局事件的累积发生函数CIF_k(t)=Pr(T≤t,D=k),即时刻t前因原因k发生失效的概率，可用特定原因的风险函数表达为：

当t趋于无穷大时，CIF_k=P(D=k)。

此处各学习材料上都会花笔墨解释为什么CIF_k≠1-S_k(t)。为什么呢？因为

两者的区别在于计算CIF用到的总体的生存函数S(t)，而1-S_k(t)计算用到的是S_k(t)，也就是说k结局事件的CIF不光与k结局事件的风险函数相关， 而且与每一个结局的风险函数相关，而1-S_k(t)则只与第k个结局事件的风险函数相关。另外，因为S(t)≤S_k(t),所以CIF_k≤1-S_k(t)，只有当不存在竞争风险时CIF_k=1-S_k(t)。

由于CIF_k与每个原因的风险函数都有关系，通过特定原因的风险函数h_k(t)来计算CIF_k比较麻烦。

注：I_k(t)即CIF_k，λ_k(t_j)即时刻j的h_k(t)。

1999年Fine和Gray提出累积发生率回归模型（CIF回归模型，也被称为Fine-Gray模型）用于直接估计协变量对目标结局事件累积发病率的影响，该模型定义了子分布风险函数（subdistribution hazard function）h_sdk(t)：

表示t时刻未发生第k类结局事件的观察人群中，k结局事件的瞬时发生率（瞬时概率强度）。与特定原因的风险函数不同，子分布风险函数将发生了其他结局事件的观察人群仍然在风险集中。此公式中的F_k(t)=CIF_k=P(T≤t,D=k)、S_k(t)=P(T>t,D=k)、fk(t)=dF_k(t)/dt分别是子分布累积发生率、子生存函数和子密度函数，注意P(D=k)=S_k(t)+F_k(t)。

此时t时间内结局事件k的累积风险函数CIF_k=F_k(t)与常规生存分析相同：

很明显，暴露因素对k结局事件发生率的影响方向与其对子分布风险函数/子分布累积分布函数的影响方向一致。

---

当考虑协变量的影响时（多因素分析），竞争风险模型常常考虑特定原因的风险回归模型和子分部的风险回归模型。

特定原因的风险回归模型：

其中h_k(t)即我们前面讲过的特定原因的风险函数（cause-specific hazard function）：

h_k0(t)表示k结局的基线风险，而β表示协变量对k结局的影响。

在经典的Cox比例风险模型中，我们通过expβ获得HR，HR可作为考察协变量对风险函数、累积风险函数、累积发生率的相对影响的指标。同样采用特定原因的风险回归模型可以获得特定原因的风险比HR_k=expβ，HR_k表示协变量每增加一个单位，k结局事件的特定原因风险比的变化。当要分析协变量如何影响特定原因的结局事件这种病因学问题，将竞争风险视为删失是可行的。由于存在竞争风险时，目标结局事件的CIF与每一个结局的风险函数相关，因此HR_k不能作为协变量对目标结局事件的CIF的影响指标。某个解释变量的出现增加了目标结局事件的发生风险，但存在竞争风险时，风险函数的增加并不能等同于累积事件发生率的增加。也就是说当存在竞争风险时，某个暴露因素导致了结局事件的风险增加，但不一定会导致其累积发生率CIF的增加。

Fine和Gray提出的累积发生率回归模型（CIF回归模型）解决了这个问题，该模型也被称为Fine-Gray模型，用于直接估计协变量对目标结局事件累积发病率的影响。

CIF回归模型：

其中h_sdk(t)即前面讲过的子分布风险函数（subdistribution hazard function）：

采用CIF回归模型可以获得子分布风险比（HR_sdk），即HR_sdk=expβ，协变量每增加一个单位，k结局事件的子分布风险比变化。由于CIF_k与子分布风险h_sdk(t)直接相关，分布风险率（HR_sdk）可以作为协变量影响CIF_k的直接度量指标。

特定原因的风险回归模型一般用于解决流行病病因学问题，而CIF回归模型则多用于回答于临床绝对发生率问题。

参考文献：

H. Putter,et al.Tutorial in biostatistics: Competing risks and multi-state models.Statist. Med. 2007; 26:2389–2430.

Bryan L,et al.Competing Risk Regression Models for Epidemiologic Data.Am J Epidemiol. 2009 Jul 15;170(2):244-56.

Peter CA,et al.Introduction to the Analysis of Survival Data in the Presence of Competing Risks.Circulation. 2016;133:601–609. 

示例来自《生存分析之Kaplan-Meier曲线绘制与比较》，数据略作调整，删除变量中有缺失值的个案。本次笔记先进行单因素分析，比较不同病理学分期的多发性骨髓瘤死于癌症和非癌症的累积发生率。然后再进行多因素分析，考察考虑年龄、性别、种族、家族遗传后，病理学分期对累积发生率的影响，多因素分析可以采用特定原因的风险回归和CIF回归，当然这两种方法也可用于单因素分析。

软件采用R4.03，程序包主要用到cmprsk和survival。

cmprsk可用于Subdistribution Analysis of Competing Risks，用到函数cuminc、crr、plot等；survival可用于生存分析，主要用到函数coxph。

 【2】R操作示例：累积发生率单因素比较

cuminc{cmprsk}:Cumulative Incidence Analysis,Estimate cumulative incidence functions from competing risks data and test equality across groups.程序包cmprsk用于Subdistribution Analysis of Competing Risks，其通过函数cuminc获得的CIF也是子分布风险函数来获得的。

cuminc(ftime, fstatus, group, strata, rho=0, cencode=0,subset, na.action=na.omit)

命令清单如下：

library(haven)

surv <- read_dta("D:/Temp/survivalm.dta") #从stata文件中导入数据

library(cmprsk)

attach(surv)  #将sury附加到搜索路径，这样后面直接使用sury中的变量名称即可，而不用采用surv$变量名的形式

CIF<-cuminc(time,cause,stage,0) #CIF<-cuminc(ftime=time,fstatus=cause,group=stage,cencode=0)

CIF #print(CIF)

legend<-c("stage I_CancerRelated", "stage II_CancerRelated", "stage III_CancerRelated", "stage I_notCancerRelated", "stage II_notCancerRelated", "stage III_notCancerRelated")  #图例命名

plot(CIF, curvlab=legend,ylim=c(0,0.5),xlab="time_days", ylab="cumulative incidence function", lty=1:length(CIF), color=1:6)

detach(surv)

Gray检验结果显示不同病理学分期的多发性骨髓瘤死于癌症的累积发生率不同（Z=27.701，P<0.001）,死于非癌症的累积发生率不同（Z=24.552，P<0.001）。

【3】R操作示例：特定原因风险回归

能进行经典的Cox风险比例回归分析的软件都可以进行特定原因的风险回归，只需要对不同原因的结局事件分别拟合即可，拟合的时候把对应的发生竞争结局的观察个体在发生竞争事件时视作删失来处理。R中可采用比较经典的coxph {survival}来处理。

Surv{survival}：Create a Survival Object.

Surv(time, time2, event,type=c('right', 'left', 'interval', 'counting', 'interval2', 'mstate'), origin=0)

该函数的难点在于定义失效事件，event: The status indicator, normally 0=alive, 1=dead. Other choices are TRUE/FALSE (TRUE = death) or 1/2 (2=death). For interval censored data, the status indicator is 0=right censored, 1=event at time, 2=left censored, 3=interval censored. For multiple endpoint data the event variable will be a factor, whose fifirst level is treated as censoring. Although unusual, the event indicator can be omitted, in which case all subjects are assumed to have an event.

建议通过event参数直接采用逻辑码“==”来指定失效事件。The function tries to distinguish between the use of 0/1 and 1/2 coding for censored data via the condition if (max(status)==2). If 1/2 coding is used and all the subjects are censored, it will guess wrong. In any questionable case it is safer to use logical coding, e.g., Surv(time,status==3) would indicate that ’3’ is the code for an event. For multi-state survival the status variable will be a factor, whose fifirst level is assumed to correspond to censoring.

coxph {survival}：Fit Proportional Hazards Regression Model.

coxph(formula, data=, weights, subset,na.action, init, control,ties=c("efron","breslow","exact"),

singular.ok=TRUE, robust,model=FALSE, x=FALSE, y=TRUE, tt, method=ties,id, cluster, istate, statedata, nocenter=c(-1, 0, 1), ...)

命令清单如下：

library(haven)

surv <- read_dta("D:/Temp/survivalm.dta") #从stata文件中导入数据

surv$stage<-factor(surv$stage)  #gender、age、relative_with_cancer_history都是采用0/1赋值的二分类变量，不设为因子不影响结果

library(survival)

coxph(Surv(time,cause==1)~gender+age+relative_with_cancer_history+stage,data=surv) #癌症死亡结局

coxph(Surv(time,cause==2)~gender+age+relative_with_cancer_history+stage,data=surv) #非癌症死亡结局

结果显示，在校正其他变量的作用后，相对于stage1，stage2的死于癌症估计系数是1.09，是一个风险因子（HR=2.976,P=0.003）;stage3也是死于癌症的一个风险因子（HR=5.367,P<0.001）。而对于非癌症死亡，stage3相对stage1是一个风险因子（HR=3.503,P<0.001），stage2却不是（HR=1.165,P=0.702）。值得注意的是，年龄是非癌症死亡的危险因子（HR=1.053,P<0.001），但不是癌症死亡的危险因子（HR=1.00,P=0.922）。

【4】R操作示例：子分布风险回归/累积发生率回归

累积发生率回归可以使用crr{cmprsk}：

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

该函数比较令人讨厌的地方是没有提供公式模型。如果要进行多因素分析，需要构建协变量矩阵，因子变量还需要设置哑变量。可借助函数model.matrix() 来完成，model.matrix() 可将对象中的因子变量转换成0/1哑变量，而连续变量保持不变，该函数可参见《哑变量》一文。

命令清单如下：

library(haven)

surv <- read_dta("D:/Temp/survivalm.dta") #从stata文件中导入数据

surv$stage<-factor(surv$stage)  #gender、age、relative_with_cancer_history都是采用0/1赋值的二分类变量，不设为因子不影响结果

library(cmprsk)

multvar<-surv[2:6]

#multvar<-surv[c("age","gender","race","relative_with_cancer_history","stage")]

attach(multvar) #将multvar附加到搜索路径

cov<-model.matrix(~age+gender+race+relative_with_cancer_history+stage)[,-1] #构建多因素矩阵，model.matrix(~A+B+…)[,-1] 将为因子编码生成变量0/1哑变量，并保持连续变量不变。最后的[,-1]从 model.matrix 的输出中删除常数项

detach(multvar)

sdr1<-crr(surv$time,surv$cause,cov, failcode=1, cencode=0)

summary(sdr1, conf.int = 0.95) #结局1

sdr2<-crr(surv$time,surv$cause,cov, failcode=2, cencode=0)

summary(sdr2, conf.int = 0.95) #结局2

结果显示，在校正其他变量的作用后，相对于stage1，stage2和stage3的死于癌症的发生率更高（子分布风险比HR2=2.927，P=0.0033；HR3=4.470，P<0.001）；相对于stage1，stage2的死于非癌症的发生率变化不大（子分布风险比HR2=1.124，P=0.77），stage3的死于非癌症的发生率更高（子分布风险比HR3=3.157，P<0.001）。

转自个人微信公众号【Memo_Cleon】的统计学习笔记：R笔记：生存分析之竞争风险模型[概念与实操]。

END