【原】生存分析之不满足风险比例假定的竞争风险模型：STATA还是R，总有一款适合你！

Memo_Cleon 2022-07-10 发布于上海

展开全文

示例：从TCGA（The Cancer Genome Atlas, 癌症基因组图谱，https://portal.gdc./）下载的多发性骨髓瘤CoMMpass研究中的临床资料（2021年9月），数据经过了清理。变量包括id、age（年龄）、gender（性别，0=female；1=male）；race（种族：0=white；1=others）、FamHist（家族癌症史：0=No；1=Yes）、stage（病理学分期：1=I期；2=II期；3=III期）、Trt（治疗方案，0=药物联合治疗，1=制剂+干细胞移植），status（状态：0=alive；1=dead）、cause（死亡原因：0=Alive；1=Not Cancer Related；2=Cancer Related）；time（确诊后生存时间，days）。治疗多发性骨髓瘤时在药物治疗中加入干细胞移植有助于降低患者的癌症死亡率？

疾病介绍：中国多发性骨髓瘤诊治指南（2022年修订），https://share./zLN0Uqi9，密码：Memocl

严格来说除了以全因死亡为研究结局，以死因、复发等为结局的研究都广泛存在着竞争风险。非癌症相关死亡的发生的时候，癌症死亡就不可能出现了，也就是说癌症相关死亡和非癌症相关死亡互为竞争事件。<<非比例风险的Cox回归模型_分层分析法>>中没有使用全因死亡作为死亡事件只是为了示例非比例风险的Cox回归模型的分层分析法，因为如果以全因死亡作为研究结局，stage等变量都满足风险比例假定，就不需要进行分层分析了。

生存分析的竞争风险模型笔者前面也做过一次长篇笔记：《R笔记：生存分析之竞争风险模型[概念与实操]》。竞争风险模型有两种分析方法：（1）特定原因风险模型，（2）子分布风险模型。能进行经典的Cox风险比例回归分析的软件都可以进行特定原因的风险回归，只需要对不同原因的结局事件分别拟合、将竞争结局视作删失来处理即可。特定原因的风险模型一般用于解决流行病病因学问题，某因素的特定原因的HR>1意味着该因素是一个风险因子，但该因素不一定会导致结局事件累积发生率升高。子分布CIF回归模型则多用于回答于临床绝对发生率问题，某因素的子分布HR>1意味着该因素导致结局事件累计发生率升高。这些我们在《R笔记：生存分析之竞争风险模型[概念与实操]》一文的理论部分做过详细的介绍。

不论是特定原因风险模型还是子分布风险模型都是需要满足风险比例假定的，两者的翻译可见一斑：Cause-specific Proportional Hazards Model、Subdistribution Proportional Hazards Model。遇到不满足风险比例假定的情况，可以考虑时依系数法、分段模型法和分层分析法，可分别参见《非比例风险的Cox回归模型_时依系数法》、《非比例风险的Cox回归模型_分段模型》和《非比例风险的Cox回归模型_分层分析法》。同样的，这些处理方法也适用于竞争风险模型。笔者在《非比例风险的Cox回归模型_分层分析法》文中的操作方法实际上可以看做是特定原因的风险模型分析的一部分，只是对结果的解读上会有一些变化：风险高不代表结局发生率高。

个人觉得STATA的竞争风险回归做得非常棒，这个过程实际上是子分布风险模型。当遇到不满足风险比例假定的情况，只需要启用其中的时变协变量（time-varying covariate）选项与相应的时间函数即可。

分析前可以先查看一下数据的基本情况：

分析>>生存分析>>摘要统计，检验和表格>>生存时间数据汇总：

【生存设置】声明生存时间数据，时间变量选入time，失效变量选入cause，失效值输入2（Cancer Related Dead），确定。

确定；

分组单独摘要选入Trt，确定；

分组单独摘要选入stage，确定。

结果中首先显示对应菜单操作的stata命令。数据集共纳入631例，存在76个失效事件，共有575747个分析人次，随访时间最长个体是1984天。

本例删除了缺失值，实际分析631例。由于较早出现了截尾数据，总体中位生存时间无法估计，Ⅲ期中位生存时间1753天。

我们仍以stage不满足风险比例假定为例，进行竞争风险回归：

统计>>生存分析>>回归模型>>竞争风险回归

【生存设置】声明生存时间数据，时间变量选入time，失效变量选入cause，失效值输入2（Cancer Related Dead）。在查看数据信息声明过了此处就不用再额外声明了。

自变量选入：age、gender、race、FamHist、stage和Trt。注意stage以因子变量纳入，可直接输入i.stage，或者通过11~14步来菜单操作；

竞争风险事件中变量选择cause，数值输入1；

选中时变协变量复选框，选入stage，对应命令tvc中指定的变量；时变协变量的乘数本例保持不变，这个参数代表了时间变量进入时依协变量的形式，对应命令中的texp参数，默认是不做变换，如果做对数变换可以输入ln(_t)，还可通过该参数依据分界点进行分段模型的建立，比如可以输入: _t>1000，可将模型分为两段。

默认结果显示子风险比率SHR，如果想显示系数，可在主对话框中的[报告]选项卡中选择[报告系数，而不是子风险比率]的复选框。

实际分析631例，发生失效事件76个，竞争事件62个，删失293个。伪似然值=-432.41，Chi2=44.39，P<0.001，表明模型是有统计学意义的，至少有一个变量的系数不为0。在校正其他因素的作用后，治疗多发性骨髓瘤时在药物治疗中加入干细胞移植确实有助于降低患者的癌症死亡率（Z=-3.69，P<0.001）。分段模型的逻辑、分界点、PH假定检验等可参考笔者在《非比例风险的Cox回归模型_分段模型》中的介绍，结果详细解读可参考《R笔记：生存分析之竞争风险模型[概念与实操]》，就不再赘述啦。

有读者留言，想用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)

几个关键参数说明：cov1指定的是需要分析和校正的因素，cov2指定的是其实就是时变协变量（类似于stata中tvc参数），而tf则是时间变量的函数，该函数与时变协变量的乘积作为时依协变量纳入模型（类似于stata中texp参数）；cengroup指定分层变量；failcode指定感兴趣的结局编码, cencode指定删失编码；na.action默认对有缺失值的个案进行删除处理；gtol迭代停止收敛容差；maxiterinit是允许的最大迭代次数；init指定回归参数初始值。

library(readxl)

MM<- read_excel("D:/Temp/PCTC20210924.xlsx")

MM<-na.omit(MM)

MM$stage<-factor(MM$stage)

#age、gender、race、FamHist、Trt都是采用0/1赋值的二分类变量，不设为因子不影响结果

library(cmprsk)

library(survival)

mm<-MM[c("age","gender","race","FamHist","stage","Trt")]

cov<-model.matrix(~age+gender+race+FamHist+stage+Trt,data=mm)[,-1]

#crr函数要求变量的输入形式是矩阵。model.matrix将为因子编码生成变量0/1哑变量，并保持连续变量不变。最后的[,-1]从model.matrix的输出中删除常数项

head(cov)

sdcr<-crr(ftime=MM$time,fstatu=MM$cause,cov1=cov, cov2=MM["stage"], tf=function(uft) uft, failcode=2, cencode=0)

#tf构造了一个时间变量函数uft，这个函数是关于失效时间的，该命令中不做任何变换，直接与cov2中指定的变量stage进行相乘来构建进入模型的时依协变量。模型中将包含age、gender、race、FamHist、stage哑变量、Trt以及一个关于stage的时依协变量

#特别注意cov2=MM["stage"] ，而不能使用cov2=cov[,5]。因为cov矩阵中stage被编码成了哑变量，cov[,5]相当于stage2。变量Trt也不再是cov[,6]而是cov[,7]

summary(sdcr,conf.int=0.95)

结果中的stage*tf1就是时依协变量，同前面stata的运行结果完全一致。在校正其他因素的作用后，治疗中加入干细胞移植有助于降低多发性骨髓瘤患者的癌症死亡率

（Z=-3.689，P=0.00023）。

同样，死于非癌症相关死亡也可以做出竞争风险分析：

sdcr2<-crr(ftime=MM$time,fstatu=MM$cause,cov1=cov, cov2=MM["stage"], tf=function(uft) uft, failcode=1, cencode=0)

summary(sdcr2,conf.int=0.95)

下面这些时依协变量的构造方法也很有用：

##构造age与时间平方的时依协变量age*time+age*time²

sdcr_age2<-crr(ftime=MM$time,fstatu=MM$cause,cov1=cov, cov2=cbind(cov[,1],cov[,1]), tf=function(uft) cbind(uft,uft*uft), failcode=2, cencode=0)

##构造stage的时依协变量，失效时间经过自然对数变换

sdcr_stage<-crr(ftime=MM$time,fstatu=MM$cause,cov1=cov, cov2=cov2=MM["stage"], tf=function(uft) log(uft), failcode=2, cencode=0)

##同时构建多个时依协变量，失效时间不做变换

sdcr_gender.Trt<-crr(ftime=MM$time,fstatu=MM$cause,cov1=cov, cov2=cov[,c(2,7)], tf=function(uft) cbind(uft,uft), failcode=2, cencode=0)

#sdcr_age.stage.Trt<-crr(ftime=MM$time,fstatu=MM$cause,cov1=cov, cov2=cbind(cov[,1],cov2=MM["stage"],cov[,7]), tf=function(uft) cbind(uft,uft,uft), failcode=2, cencode=0)

时依协变量也用于风险比例假定的检验：

sdcr.phtest<-crr(ftime=MM$time,fstatu=MM$cause,cov1=cov, cov2=cbind(cov[,-c(5:6)],MM["stage"]), tf=function(uft) cbind(uft,uft,uft,uft,uft,uft), failcode=2, cencode=0)