使用SAS軟件分析競爭風(fēng)險模型

北醫(yī)仁智（北京）醫(yī)學(xué)科技發(fā)展有限公司醫(yī)學(xué)統(tǒng)計中心（100029） 陶 莊

·計算機應(yīng)用·

使用SAS軟件分析競爭風(fēng)險模型

北醫(yī)仁智（北京）醫(yī)學(xué)科技發(fā)展有限公司醫(yī)學(xué)統(tǒng)計中心（100029） 陶 莊

在生存分析中，無疑Kaplan-Meier（KM）估計是生存函數(shù)等指標(biāo)最流行的非參數(shù)估計方法，在實施這種方法的數(shù)據(jù)中，一部分出現(xiàn)了研究者感興趣的結(jié)局，而另一部分則沒有出現(xiàn)，這后一部分被處理成為刪失數(shù)據(jù)（censored data）。但在醫(yī)療實踐中，縱向數(shù)列并不總是僅僅出現(xiàn)研究者感興趣的事件，或稱真結(jié)局（true outcome），還會出現(xiàn)一些并不感興趣的結(jié)局。比如對于一個骨髓移植后的病人來說，復(fù)發(fā)是真結(jié)局，但是如果某個病人在復(fù)發(fā)前就死亡了，而死亡就不是研究者感興趣的結(jié)局，但是死亡將使得復(fù)發(fā)無法出現(xiàn)，因此死亡這個事件就成為了一個競爭事件（competing event），更為通俗的叫法是：死亡是復(fù)發(fā)的競爭風(fēng)險（competing risk）（其實更確切的說法是雙方互為競爭風(fēng)險）［1］。

對于這樣帶有競爭風(fēng)險的數(shù)據(jù)，早先的做法是將競爭事件也定義為刪失，然后直接使用KM估計。然而競爭事件與刪失數(shù)據(jù)還是有很大差別的：刪失是指數(shù)據(jù)的真結(jié)局沒有被觀察到，但該病例仍有可能發(fā)生真結(jié)局，只不過研究者不知道而已；而競爭事件是因為該事件的出現(xiàn)，使得真結(jié)局確實無法出現(xiàn)。如果此時將競爭事件簡單視為刪失數(shù)據(jù)，將會高估真結(jié)局的發(fā)生率［1］，也就是出現(xiàn)了所謂的競爭風(fēng)險偏倚（competing risk bias）。

競爭風(fēng)險偏倚頻繁地出現(xiàn)在醫(yī)學(xué)研究的文獻(xiàn)中，Koller等觀察了35篇使用KM估計的高影響因子的醫(yī)學(xué)論文，其中24篇（67%）被認(rèn)為可能存在競爭風(fēng)險偏倚［2］；而C van Walraven等人在一項專門對競爭風(fēng)險偏倚進行的研究中發(fā)現(xiàn)，有46%的文獻(xiàn)可能存在這種偏倚，同樣的，他們的樣本也來自于醫(yī)學(xué)領(lǐng)域的一些高分文獻(xiàn)［3］。

對于競爭風(fēng)險模型，通常使用的終點指標(biāo)是累計發(fā)生率函數(shù)（cumulative incidence function，CIF），最經(jīng)典的估計方法來自于Kalbfleisch和Prentice的著作［4］，它的表達(dá)式為：

其中j為某種結(jié)局，而k和i分別是分組與病人標(biāo)號。具體的估計方法通常采用下式：

由此Fine和Gray在1999年提出其危險率函數(shù)為：

1988年Gray在該定義的基礎(chǔ)上給出了對CIF在各組進行組間比較的檢驗方法，被稱為Gray檢驗［5］。Gray檢驗的理論過程相當(dāng)復(fù)雜，關(guān)鍵是構(gòu)建一個得分統(tǒng)計量：

然后根據(jù)這個得分的方差協(xié)方差矩陣進行計算，由于篇幅有限，此處不再展開，有興趣的讀者可以參看Gray的原文。

而對于帶有協(xié)變量的CIF檢驗，則是基于傳統(tǒng)的Cox比例風(fēng)險模型，它是由Fine和Gray在1999年提出的［6］。此時：

此處j是結(jié)局類型，Z是協(xié)變量，β是其系數(shù)。而對CIF的估計公式為：

而其中：

這里的G（X）指的是相應(yīng)的Kaplan-Meier估計。

可以說，目前統(tǒng)計軟件使用的，進行競爭風(fēng)險模型的方法，主要基于以上這三篇文獻(xiàn)。

在我國，明確使用競爭風(fēng)險模型進行分析的中文醫(yī)學(xué)文獻(xiàn)寥寥無幾，而介紹這種模型的文獻(xiàn)也不多見，這些都可以通過CNKI等網(wǎng)站清晰地看到。有學(xué)者在2008年曾寫過一篇介紹如何使用R軟件進行分析的文章［7］，那時的SAS只能使用宏（macro）分析競爭風(fēng)險模型，并不方便，而這種形勢直到2013年SAS 9．4的發(fā)布才有所改變。本文并非是對競爭風(fēng)險模型的理論介紹，而是將SAS如何進行該模型分析的發(fā)展歷程與實戰(zhàn)方法呈現(xiàn)給大家，以期對一線科研工作者有所幫助。

示例數(shù)據(jù)

為了更好地說明SAS進行競爭風(fēng)險模型的方法和步驟，我們引用一個示例數(shù)據(jù)，它是一套真實研究數(shù)據(jù)的節(jié)選。這里的數(shù)據(jù)集僅包括3個變量55條觀測，由于是示例數(shù)據(jù)，其變量含義并非十分重要。數(shù)據(jù)庫的結(jié)構(gòu)見表1。

表1 示例數(shù)據(jù)結(jié)構(gòu)表

競爭風(fēng)險模型SAS三部曲

第一部：初探——%Cum Incid

2002年，SAS公司推出了基于LIFETEST過程的內(nèi)部（SAS公司自己開發(fā)的）宏程序%Cum Incid，這基本上可以認(rèn)為是SAS初次涉入競爭風(fēng)險模型的分析領(lǐng)域。%Cum Incid的運行環(huán)境是SAS 9．1，其調(diào)用形式與SAS其他的宏沒有不同，形式如下：

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

其中各參數(shù)pi的形式與意義見表2。

表2 %Cum Incid中的參數(shù)形式及意義

對于本例，我們可以鍵入：

%Cum Incid的結(jié)果主要包括三個部分：

第一個部分是基于LIFETEST過程進行的一次KM分析，這里模型將競爭事件也作為刪失數(shù)據(jù)處理，此時獲得的壽命表中的Failure一項（即1-KM）有時也被稱為真結(jié)局的原因別（cause-specific event，CSE）發(fā)生率。

第二個部分與第一部分類似，所不同的是將競爭事件也作為真結(jié)局處理。

第三個部分是真正的CIF。在這個部分里包括CIF的估計值列表以及CIF的統(tǒng)計圖（圖1），此時雖使用了PLOTCL選項，但曲線的置信帶并非直接加到原圖上，而是在不同的圖中分層顯示。需注意的是，“htm file”與“rtffile”必選其一，否則圖形將被屏蔽。如果在選項中選擇“CSE”，則在CIF表單中將同時包含第一部分的CSE內(nèi)容，此時，使用者可以清楚地看到所謂的競爭風(fēng)險偏倚。

圖1 %Cum Incid所作出的CIF統(tǒng)計圖（并不顯示置信帶）

目前%Cum Incid仍然可以在SAS Support上進行下載，下載的網(wǎng)址是：

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

在SAS Support的網(wǎng)站上，對這個宏幾乎沒有什么像樣的介紹，由此可以看出SAS在初涉該領(lǐng)域時的粗陋。而且該宏沒有進行Gray檢驗的功能，這也極大的限制了%Cum Incid在實際中的應(yīng)用。畢竟，誰又愿意用SAS估計完CIF后再用R跑一遍檢驗?zāi)兀?/p>

第二部：完善——%CIF

一直等待了10年，SAS才推出了%Cum Incid的升級產(chǎn)品%CIF。%CIF同樣也是基于LIFETEST過程，其運行環(huán)境為SAS 9．2。SAS對%CIF的更新一直持續(xù)到2014年，考慮到此時SAS 9．4已經(jīng)面世，顯然該宏也可以順利運行于9．4版。%CIF在樣子上只是稍作改動，但是功能上終于加入了大眾期盼已久的Gray檢驗，%CIF的下載地址是：

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

而且這次SAS顯然非常重視%CIF，在其Support網(wǎng)站上的支持力度是%Cum Incid完全不能比擬的，甚至，SAS還專門為其撰寫了應(yīng)用文章［8］。%CIF的調(diào)用形式與%Cum Incid類似：

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

其中各參數(shù)pi的形式與意義見表3。

同樣的，對于本例，我們可以在SAS程序編輯器里鍵入：

結(jié)果也與%Cum Incid類似，只不過在最后加入了Gray檢驗的結(jié)果。

表3 %CIF中的參數(shù)形式及意義

圖2 %CIF所作出的CIF統(tǒng)計圖（直接顯示置信帶）

表4 Gray檢驗的結(jié)果

盡管在10年的等待后，在功能中加入了Gray檢驗，但是SAS并未同時推出一款可以進行調(diào)整協(xié)變量的宏。

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

2013年7月10日，翹首以待的SAS 9．4終于面世了，在這個加入眾多分析的新版本中就包含有競爭風(fēng)險模型。這次，SAS基于的是PHREG過程，語句和選項都非常簡潔，并且還增加了直接繪圖功能，大大簡化了分析步驟。SAS這次根本沒有觸及LIFETEST過程，因為在當(dāng)model語句中只包含一個自變量的情況下，獲得的估計就是前述的CIF［9］！

以下就直接使用示例數(shù)據(jù)介紹PHREG過程中相關(guān)基礎(chǔ)選項，更多的選項和相關(guān)結(jié)果大家可以通過SAS幫助學(xué)習(xí)試用。

其中（3）是主程序，但是（1）（2）是必須的，因為如果沒有這兩部分程序，結(jié)果將僅出現(xiàn)回歸系數(shù)估計及檢驗（也就是Gray檢驗）的結(jié)果，而沒有CIF估計值及圖形。

在新的PHREG中，進行競爭風(fēng)險模型的選項主要位于三個地方：

1．在proc phreg語句中增加了“plots＝”選項，用來描繪CIF圖，其中的選項overlay＝stratum表示將各組的曲線置于一張圖中。

2．在model選項中增加了“eventcode＝”選項，指明cg中哪個取值是真結(jié)局。

3．在baseline語句中加入了“cif＝”選項，此時將產(chǎn)生CIF的估計值并保存在“out＝”的數(shù)據(jù)集中，＿all＿相當(dāng)于定制了CIF的估計值、CIF的標(biāo)準(zhǔn)誤、以及95% CI的置信限4部分內(nèi)容。

表5 參數(shù)估計及檢驗的結(jié)果

外部宏

SAS允許使用者開發(fā)自己的宏程序并搭載在SAS上使用。在這些針對競爭風(fēng)險模型分析的宏中，比較有影響力的有：%Cum Inc（2003）［10］，%CUM INC POISSON（2008）［11］，%CIFCOX與%CIFSTRATA（2010）［12］，以及%PSHREG（2010）［13］。這些宏隨著開發(fā)年代的不同，適應(yīng)的SAS環(huán)境也不盡相同，方法也各有千秋，有興趣的讀者可以自行查找相關(guān)的文獻(xiàn)，這里不做詳論。

圖3 PHREG過程中plots所作出的CIF統(tǒng)計圖

總 結(jié)

中國有句古話叫“十年磨一劍”，這話用在SAS開發(fā)競爭風(fēng)險模型的歷程上無疑是非常恰當(dāng)?shù)?，但不管怎么說，現(xiàn)在進行競爭風(fēng)險模型已不再是R專美之事。本文除了介紹最新的SAS 9．4的相關(guān)內(nèi)容外，一并介紹了%Cum Incid和%CIF，既是考慮到有一個SAS研發(fā)的完整性問題，也是因為這兩個宏各自可適應(yīng)的SAS版本不同，也許仍然可以解決我國相當(dāng)一部分科研人員的實際問題。

［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軟件分析競爭風(fēng)險模型簡明攻略．中國衛(wèi)生統(tǒng)計，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．

（責(zé)任編輯：劉 壯）