使用SAS软件分析竞争风险模型

北医仁智（北京）医学科技发展有限公司医学统计中心（100029） 陶 庄

·计算机应用·

使用SAS软件分析竞争风险模型

北医仁智（北京）医学科技发展有限公司医学统计中心（100029） 陶 庄

在生存分析中，无疑Kaplan-Meier（KM）估计是生存函数等指标最流行的非参数估计方法，在实施这种方法的数据中，一部分出现了研究者感兴趣的结局，而另一部分则没有出现，这后一部分被处理成为删失数据（censored data）。但在医疗实践中，纵向数列并不总是仅仅出现研究者感兴趣的事件，或称真结局（true outcome），还会出现一些并不感兴趣的结局。比如对于一个骨髓移植后的病人来说，复发是真结局，但是如果某个病人在复发前就死亡了，而死亡就不是研究者感兴趣的结局，但是死亡将使得复发无法出现，因此死亡这个事件就成为了一个竞争事件（competing event），更为通俗的叫法是：死亡是复发的竞争风险（competing risk）（其实更确切的说法是双方互为竞争风险）［1］。

对于这样带有竞争风险的数据，早先的做法是将竞争事件也定义为删失，然后直接使用KM估计。然而竞争事件与删失数据还是有很大差别的：删失是指数据的真结局没有被观察到，但该病例仍有可能发生真结局，只不过研究者不知道而已；而竞争事件是因为该事件的出现，使得真结局确实无法出现。如果此时将竞争事件简单视为删失数据，将会高估真结局的发生率［1］，也就是出现了所谓的竞争风险偏倚（competing risk bias）。

竞争风险偏倚频繁地出现在医学研究的文献中，Koller等观察了35篇使用KM估计的高影响因子的医学论文，其中24篇（67%）被认为可能存在竞争风险偏倚［2］；而C van Walraven等人在一项专门对竞争风险偏倚进行的研究中发现，有46%的文献可能存在这种偏倚，同样的，他们的样本也来自于医学领域的一些高分文献［3］。

对于竞争风险模型，通常使用的终点指标是累计发生率函数（cumulative incidence function，CIF），最经典的估计方法来自于Kalbfleisch和Prentice的著作［4］，它的表达式为：

其中j为某种结局，而k和i分别是分组与病人标号。具体的估计方法通常采用下式：

由此Fine和Gray在1999年提出其危险率函数为：

1988年Gray在该定义的基础上给出了对CIF在各组进行组间比较的检验方法，被称为Gray检验［5］。Gray检验的理论过程相当复杂，关键是构建一个得分统计量：

然后根据这个得分的方差协方差矩阵进行计算，由于篇幅有限，此处不再展开，有兴趣的读者可以参看Gray的原文。

而对于带有协变量的CIF检验，则是基于传统的Cox比例风险模型，它是由Fine和Gray在1999年提出的［6］。此时：

此处j是结局类型，Z是协变量，β是其系数。而对CIF的估计公式为：

而其中：

这里的G（X）指的是相应的Kaplan-Meier估计。

可以说，目前统计软件使用的，进行竞争风险模型的方法，主要基于以上这三篇文献。

在我国，明确使用竞争风险模型进行分析的中文医学文献寥寥无几，而介绍这种模型的文献也不多见，这些都可以通过CNKI等网站清晰地看到。有学者在2008年曾写过一篇介绍如何使用R软件进行分析的文章［7］，那时的SAS只能使用宏（macro）分析竞争风险模型，并不方便，而这种形势直到2013年SAS 9．4的发布才有所改变。本文并非是对竞争风险模型的理论介绍，而是将SAS如何进行该模型分析的发展历程与实战方法呈现给大家，以期对一线科研工作者有所帮助。

示例数据

为了更好地说明SAS进行竞争风险模型的方法和步骤，我们引用一个示例数据，它是一套真实研究数据的节选。这里的数据集仅包括3个变量55条观测，由于是示例数据，其变量含义并非十分重要。数据库的结构见表1。

表1 示例数据结构表

竞争风险模型SAS三部曲

第一部：初探——%Cum Incid

2002年，SAS公司推出了基于LIFETEST过程的内部（SAS公司自己开发的）宏程序%Cum Incid，这基本上可以认为是SAS初次涉入竞争风险模型的分析领域。%Cum Incid的运行环境是SAS 9．1，其调用形式与SAS其他的宏没有不同，形式如下：

%Cum Incid（p1，p2，…，pn）；

其中各参数pi的形式与意义见表2。

表2 %Cum Incid中的参数形式及意义

对于本例，我们可以键入：

%Cum Incid的结果主要包括三个部分：

第一个部分是基于LIFETEST过程进行的一次KM分析，这里模型将竞争事件也作为删失数据处理，此时获得的寿命表中的Failure一项（即1-KM）有时也被称为真结局的原因别（cause-specific event，CSE）发生率。

第二个部分与第一部分类似，所不同的是将竞争事件也作为真结局处理。

第三个部分是真正的CIF。在这个部分里包括CIF的估计值列表以及CIF的统计图（图1），此时虽使用了PLOTCL选项，但曲线的置信带并非直接加到原图上，而是在不同的图中分层显示。需注意的是，“htm file”与“rtffile”必选其一，否则图形将被屏蔽。如果在选项中选择“CSE”，则在CIF表单中将同时包含第一部分的CSE内容，此时，使用者可以清楚地看到所谓的竞争风险偏倚。

图1 %Cum Incid所作出的CIF统计图（并不显示置信带）

目前%Cum Incid仍然可以在SAS Support上进行下载，下载的网址是：

http：//support．sas．com/kb/30/511．htm l

在SAS Support的网站上，对这个宏几乎没有什么像样的介绍，由此可以看出SAS在初涉该领域时的粗陋。而且该宏没有进行Gray检验的功能，这也极大的限制了%Cum Incid在实际中的应用。毕竟，谁又愿意用SAS估计完CIF后再用R跑一遍检验呢？

第二部：完善——%CIF

一直等待了10年，SAS才推出了%Cum Incid的升级产品%CIF。%CIF同样也是基于LIFETEST过程，其运行环境为SAS 9．2。SAS对%CIF的更新一直持续到2014年，考虑到此时SAS 9．4已经面世，显然该宏也可以顺利运行于9．4版。%CIF在样子上只是稍作改动，但是功能上终于加入了大众期盼已久的Gray检验，%CIF的下载地址是：

http：//support．sas．com/kb/45/addl/fusion＿45997＿15＿cif．txt

而且这次SAS显然非常重视%CIF，在其Support网站上的支持力度是%Cum Incid完全不能比拟的，甚至，SAS还专门为其撰写了应用文章［8］。%CIF的调用形式与%Cum Incid类似：

%CIF（p1，p2，…，pn）；

其中各参数pi的形式与意义见表3。

同样的，对于本例，我们可以在SAS程序编辑器里键入：

结果也与%Cum Incid类似，只不过在最后加入了Gray检验的结果。

表3 %CIF中的参数形式及意义

图2 %CIF所作出的CIF统计图（直接显示置信带）

表4 Gray检验的结果

尽管在10年的等待后，在功能中加入了Gray检验，但是SAS并未同时推出一款可以进行调整协变量的宏。

第三部：融合——SAS 9．4之PHREG

2013年7月10日，翘首以待的SAS 9．4终于面世了，在这个加入众多分析的新版本中就包含有竞争风险模型。这次，SAS基于的是PHREG过程，语句和选项都非常简洁，并且还增加了直接绘图功能，大大简化了分析步骤。SAS这次根本没有触及LIFETEST过程，因为在当model语句中只包含一个自变量的情况下，获得的估计就是前述的CIF［9］！

以下就直接使用示例数据介绍PHREG过程中相关基础选项，更多的选项和相关结果大家可以通过SAS帮助学习试用。

其中（3）是主程序，但是（1）（2）是必须的，因为如果没有这两部分程序，结果将仅出现回归系数估计及检验（也就是Gray检验）的结果，而没有CIF估计值及图形。

在新的PHREG中，进行竞争风险模型的选项主要位于三个地方：

1．在proc phreg语句中增加了“plots＝”选项，用来描绘CIF图，其中的选项overlay＝stratum表示将各组的曲线置于一张图中。

2．在model选项中增加了“eventcode＝”选项，指明cg中哪个取值是真结局。

3．在baseline语句中加入了“cif＝”选项，此时将产生CIF的估计值并保存在“out＝”的数据集中，＿all＿相当于定制了CIF的估计值、CIF的标准误、以及95% CI的置信限4部分内容。

表5 参数估计及检验的结果

外部宏

SAS允许使用者开发自己的宏程序并搭载在SAS上使用。在这些针对竞争风险模型分析的宏中，比较有影响力的有：%Cum Inc（2003）［10］，%CUM INC POISSON（2008）［11］，%CIFCOX与%CIFSTRATA（2010）［12］，以及%PSHREG（2010）［13］。这些宏随着开发年代的不同，适应的SAS环境也不尽相同，方法也各有千秋，有兴趣的读者可以自行查找相关的文献，这里不做详论。

图3 PHREG过程中plots所作出的CIF统计图

总 结

中国有句古话叫“十年磨一剑”，这话用在SAS开发竞争风险模型的历程上无疑是非常恰当的，但不管怎么说，现在进行竞争风险模型已不再是R专美之事。本文除了介绍最新的SAS 9．4的相关内容外，一并介绍了%Cum Incid和%CIF，既是考虑到有一个SAS研发的完整性问题，也是因为这两个宏各自可适应的SAS版本不同，也许仍然可以解决我国相当一部分科研人员的实际问题。

［1］Klein JP，Moeschberger ML．Survival Analysis Techniques for Censored and Truncated Data．Second Edition．Springer-Verlag New York，Inc，2003：127-132．

［2］Koller MT，Raatz H，Steyerberg EW，et al．Competing risks and the clinical community：irrelevance or ignorance？Stat Med，2012，31：1089e97．

［3］van Walraven C，M cAlister FA．Competing risk bias is common in Kaplan-Meier estimates published in prominent medical journals．J Clin Epidem iol，2016，69：170-173．

［4］Kalbfleisch J，Prentice R．The statistical analysis of failure time data．John W iley＆Sons，New York，1980：168-169．

［5］Gray RJ．A class of K-sample tests for comparing the cumulative incidence of a competing risk．Annals of statistics，1988，16（3）：1141-1154．

［6］Fine JP，Gray RJ．A proportional hazardsmodel for the subdistribution of a competing risk．Journal of the American Statistical Association，1999，94：496-509．

［7］陶庄．使用R软件分析竞争风险模型简明攻略．中国卫生统计，2008，25（6）：80-81

［8］Lin G，So Y，Johnston J．Analyzing Survival Data w ith Competing Risks Using SAS®Software．Proceedings of the SAS®Global Forum 2012 Conference，Cary，NC：SAS Institute Inc，2012．

［9］Lin G，So Y，Johnston G．Using the PHREG Procedure to Analyze Competing-Risks Data．Proceedings of the SAS®Global Forum 2012 Conference，Cary，NC：SAS Institute Inc，2015．

［10］Rosthoj S，Andersen PK，Abildstrom S．SASmacros for estimation of the cumulative incidence functions based on acox regression model for competing risks survival data．Comput，Methods Progr．Biomed，2004，74：69-75．

［11］Waltoft BL．A SAS-macro for estimation of the cumulative incidence using Poisson regression．Comput Methods Progr Biomed，2009，93：140-147．

［12］Zhang X，Zhang MJ．SASmacros for estimation of directadjusted cumulative incidence curves under proportional subdistribution hazards models．Comput Methods Progr Biomed，2011，101：87-93．

［13］Kohl M，Plischke M，Leffondre K，et al．PSHREG：A SASmacro for proportional and nonproportional subdistribution hazards regression．Computer Methods and Programs in Biomedicine，2015，118：218-233．

（责任编辑：刘 壮）