宋德勝 ，李長平 ，2，劉媛媛 ，崔 壯 *，胡良平

（1.天津醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)教研室，天津 300070；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京 100029；3.軍事科學(xué)院研究生院，北京 100850*通信作者：崔 壯，E-mail：cuizhuang@tmu.edu.cn）

在臨床實踐中，調(diào)查者會在指定的時間段內(nèi)隨訪研究對象直到其發(fā)生預(yù)先指定的觀察事件。然而，部分研究對象會在隨訪期間出于某種原因而退出研究，例如，出現(xiàn)非觀察事件導(dǎo)致的死亡或者研究對象主動要求退出等。這些提前退出隨訪的情況被稱之為截尾。這種截尾會導(dǎo)致收集的數(shù)據(jù)不完整。因此，傳統(tǒng)的參數(shù)回歸模型并不適用于處理生存資料。目前，生存分析中最常用的回歸模型是英國統(tǒng)計學(xué)家D.R.Cox于1972年提出的Cox半?yún)?shù)回歸模型。該模型不要求生存時間滿足特定的概率分布，但要求生存資料滿足比例風(fēng)險假定。

1 比例風(fēng)險概念

經(jīng)典的Cox比例風(fēng)險回歸模型存在一個假設(shè)：不論基線風(fēng)險如何，在基線以后的任何時間點上，分別在影響因素的“暴露水平”與“非暴露水平”條件下的發(fā)生事件的風(fēng)險比是恒定的，換言之，所考察的影響因素對于所考察事件的效應(yīng)不會隨時間而改變。這就是比例風(fēng)險恒定假設(shè)［1］，簡稱為PH假設(shè)。例如，年齡為50歲時，男性發(fā)生心臟病的風(fēng)險是女性的2倍，那么60歲時，男性的風(fēng)險仍然是女性的2倍（說明：此例中的“影響因素”為“性別”，其“暴露水平”為“男性”，“非暴露水平”為“女性”）。但有很多臨床生存資料并不滿足此假設(shè)，此時，這種變量（即影響因素）效應(yīng)的風(fēng)險稱為非比例風(fēng)險。這種效應(yīng)隨時間變化的變量稱為時間依賴型變量，即時依協(xié)變量［2］。比例風(fēng)險恒定更普遍的情形是：設(shè)Xi為第i位受試對象的自變量向量、Xj為第j位受試對象的自變量向量，PH假設(shè)為第i位受試對象與第j位受試對象風(fēng)險之比僅與他們的自變量向量的取值之差呈比例關(guān)系，而與自變量向量在什么時間點取值無關(guān)。

2 使用圖示法檢驗比例風(fēng)險假定

判斷比例風(fēng)險假設(shè)是否成立的一個簡單的方法是圖示法，即通過觀察每一個定性自變量各水平條件下的Kaplan-Meier生存曲線圖是否存在交叉，如果存在交叉，則表示該定性自變量不滿足比例風(fēng)險假設(shè)。另外，對于特定的定性自變量的各水平組，繪制ln｛-ln［S（t）］｝與生存時間或生存時間的對數(shù)的關(guān)系圖，如果線段明顯不平行，說明該定性自變量不符合比例風(fēng)險假設(shè)。對于連續(xù)型自變量，可使用Schoenfeld殘差圖、Score殘差圖進行判斷，也可以將連續(xù)型自變量定性化，然后采取前述的圖示法。

對于已經(jīng)觀測到的事件時間，假如已知第i個對象的第k個協(xié)變量及其取值，則Schoenfeld殘差見式（1）［3］：

式（1）中，Xik是第i個對象第k個協(xié)變量的值，wik是給定事件時間的風(fēng)險集中協(xié)變量值的加權(quán)均數(shù)。若Schoenfeld殘差值為正，表示在對應(yīng)的死亡時間點，X的實際值高于預(yù)期值。繪制Schoenfeld殘差與生存時間的廣義線性回歸圖，若圖形呈現(xiàn)非零斜率，表示該變量不滿足比例風(fēng)險假設(shè)。

Martingale殘差定義見式（2）［4］：

式（2）中，Nj（t）代表t時刻個體j是否經(jīng)歷了某事件的指示變量，Yj(t)是t時刻前個體j是否在觀察中的指示變量，是回歸系數(shù)向量，zj(t)是t時刻，第j個個體的協(xié)變量向量，是累積基準風(fēng)險函數(shù)的Breslow估計。因此，Martingale殘差可能具有超額事件數(shù)，并且這些殘差的總和等于0。在滿足比例風(fēng)險假設(shè)的前提下，如果利用該殘差與時間作圖，則可以觀察到該殘差隨時間圍繞一條水平線波動。

通過SAS的PHREG過程中ASSESS語句繪制累積得分殘差與時間的圖形以檢驗比例風(fēng)險假設(shè)。圖形中每條曲線的值開始于0且終止于0，這是一種Brownian過程。在比例風(fēng)險的假設(shè)下，該選項產(chǎn)生若干隨機路徑。將隨機產(chǎn)生的路徑與實際數(shù)據(jù)相對比，若變量的實際路徑在隨機路徑范圍內(nèi)，則表示該協(xié)變量服從比例風(fēng)險假設(shè)。反之，則不服從比例風(fēng)險假設(shè)。但是，非比例風(fēng)險假設(shè)的形式則不清楚。

3 SAS軟件實現(xiàn)

3.1 創(chuàng)建數(shù)據(jù)集

Krall等［5］對一項多發(fā)性骨髓瘤研究的數(shù)據(jù)進行了分析。65例患者接受了烷化劑治療，其中48例在研究期間去世，17例存活。創(chuàng)建的數(shù)據(jù)集Myeloma包含變量如下：Time（預(yù)后生存時間）；Vstatus（患者狀態(tài)，0表示存活，1表示死亡）；在診斷時被認為與生存時間有關(guān)的變量，如LogBUN（血尿素氮水平）；HGB（血紅蛋白水平）；Platelet（血小板水平，0表示非正常，1表示正常）；Age（年齡）；LogWBC（白細胞水平的對數(shù)）；Frac（骨折，0表示未發(fā)生，1表示發(fā)生）；LogPBM（骨髓漿細胞對數(shù)百分比）；Protein（蛋白質(zhì)水平）；Scalc（血鈣水平）。

數(shù)據(jù)集創(chuàng)建程序如下：

3.2 比例風(fēng)險假設(shè)檢驗

3.2.1 Kaplan-Meier生存曲線以及生存函數(shù)負對數(shù)的對數(shù)與時間的對數(shù)關(guān)系圖

以下SAS語句使用圖示法中的Kaplan-Meier生存曲線以及生存函數(shù)負對數(shù)的對數(shù)與時間的對數(shù)關(guān)系圖判斷定性變量是否符合比例風(fēng)險假設(shè)。

【程序說明】ODS語句指定后續(xù)PROC LIFETEST產(chǎn)生的圖形使用HTML BLUECML的樣式顯示，圖形的DPI設(shè)置為300。PROC LIFETEST語句中，PLOTS選項的survival以及l(fā)ls要求繪制Kaplan-Meier生存曲線與生存函數(shù)負對數(shù)的對數(shù)與時間的對數(shù)關(guān)系圖。TIME語句指定生存時間與截尾指示變量（0表示截尾），STRATA語句指定需要考察是否滿足PH假定的定性變量。

圖1 不同PLATELET水平下的生存曲線

圖2 不同PLATELET水平下的LOG［-LOG（生存函數(shù)）］與LOG（time）的變化趨勢

圖3 不同F(xiàn)RAC水平下的生存曲線

圖4 不同F(xiàn)RAC水平下的LOG［-LOG（生存函數(shù)）］與LOG（time）的變化趨勢

圖1顯示的是不同platelet水平下的生存曲線；圖2顯示的是不同platelet水平下的LOG［-LOG（生存函數(shù)）］與LOG（time）之間的折線圖；圖3顯示的是不同F(xiàn)RAC水平下的生存曲線；圖4顯示的是不同F(xiàn)RAC水平下的LOG［-LOG（生存函數(shù)）］與LOG（time）之間的折線圖。圖1中，兩條生存曲線無交叉；圖2中，兩條線沒有明顯交叉的趨勢，因此可認為變量PLATELET滿足比例風(fēng)險假設(shè)；圖3中，兩條生存曲線存在交叉情況；圖4中，兩條曲線亦存在相交，因此可以認為FRAC不滿足比例風(fēng)險假設(shè)。

3.2.2 Schoenfeld殘差與log（time）的關(guān)系圖

SAS程序如下：

【程序說明】Proc phreg調(diào)用PHREG過程進行分析，data選項指定要分析的數(shù)據(jù)集是myeloma，zph選項（注意：此選項在SAS 9.3中無效）要求進行比例風(fēng)險檢驗，括號中的notest表示不進行相關(guān)檢驗，fit指定是否呈現(xiàn)光滑曲線擬合結(jié)果，本例中使用了Spline，要求進行懲罰B樣條曲線擬合。Class語句指定分類變量為platelet和frac。Model語句進行模型構(gòu)建。Time*vstatus（0）中的time表示生存時間，當(dāng)vstatus=0時代表截尾。等號右邊為一系列自變量。Selection=s表示進行逐步選擇篩選變量。

輸出結(jié)果見表1、圖5、圖6。

表1 基于最大似然法估計回歸參數(shù)

圖5 變量logBun的縮放Schoenfeld殘差與時間對數(shù)的關(guān)系圖

圖6 變量HGB的縮放Schoenfeld殘差與時間對數(shù)的關(guān)系圖

【結(jié)果解釋】最大似然估計結(jié)果顯示，逐步回歸篩選后，模型中剩余的變量為logBun和HGB。從這兩個自變量的縮放Schoenfeld殘差與時間對數(shù)的關(guān)系圖可以看出，logBun擬合的曲線斜率明顯不為0，而HGB擬合的曲線斜率基本為0。因此，LogBun不符合比例風(fēng)險假設(shè)，HGB符合比例風(fēng)險假設(shè)。

3.2.3 使用PHREG過程ASSESS語句判斷比例風(fēng)險假設(shè)是否成立

SAS程序如下：

【程序說明】ASSESS語句要求進行Cox回歸模型的充分性檢驗。通過這個語句，可檢驗一個或多個協(xié)變量的函數(shù)形式。PH選項要求進行比例風(fēng)險假設(shè)。Seed選項設(shè)置了種子數(shù)，以保證結(jié)果的重現(xiàn)性。對于用戶指定的每個協(xié)變量［通過VAR=（變量列表）指定］，ASSESS語句會繪制已觀測的累積Martingale殘差與解釋變量值的關(guān)系圖，并模擬若干殘差圖形（通過NAPTHS=n指定）。

輸出結(jié)果見圖7、圖8。

圖7 變量LogBun的模擬路徑圖

圖8 變量HGB的模擬路徑圖

【結(jié)果解釋】LogBun變量的實際路徑在模擬路徑外，而HGB變量的實際路徑均在模擬路徑范圍內(nèi)；模擬路徑圖右下角Kolmogorov-type supremum檢驗結(jié)果顯示LogBun的P＜0.05。因此，LogBun不滿足比例風(fēng)險假設(shè)；而HGB變量滿足比例風(fēng)險假設(shè)。

4 小結(jié)

Cox比例風(fēng)險回歸模型常用于分析生存數(shù)據(jù)，需滿足比例風(fēng)險假設(shè)，但這在實際的生存資料中往往不能滿足。因此，本文介紹了在臨床試驗中簡便的比例風(fēng)險假設(shè)的檢驗方法，即圖示法。Kaplan-Meier生存曲線圖和LLS生存函數(shù)負對數(shù)的對數(shù)與時間對數(shù)的關(guān)系圖是最常用的用于直觀判斷分類變量是否滿足比例風(fēng)險假設(shè)的工具；對于定量變量，通過觀察Schoenfeld殘差與時間函數(shù)的關(guān)系圖，可大致判斷定量變量是否滿足比例風(fēng)險假設(shè)。另外，本文也使用了SAS程序PHREG過程ASSESS語句中的PH選項和RESAMPLE選項檢驗比例風(fēng)險假設(shè)是否成立。