使用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．

（責任編輯：劉 壯）