宋德勝 ，李長平 ，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）

在生存分析中，生存函數(shù)是指發(fā)生事件的時間超過指定時間的概率［1］，這個概率取決于兩個方面：基線風(fēng)險函數(shù)和效應(yīng)參數(shù)。前者體現(xiàn)了在基線協(xié)變量下，事件發(fā)生的風(fēng)險如何隨時間變化；后者表示發(fā)生事件的風(fēng)險如何隨協(xié)變量變化。因此，若效應(yīng)參數(shù)隨時間發(fā)生變化，使用經(jīng)典的Cox比例風(fēng)險回歸模型估計的生存函數(shù)會存在較大偏差。鑒于上述情況，對比例風(fēng)險假設(shè)的檢驗也就成為了應(yīng)用Cox比例風(fēng)險模型的重要前提。目前，對于比例風(fēng)險的檢驗主要包括圖示法（前文已闡述，本文不再贅述）和擬合優(yōu)度檢驗法。前者雖然簡便直觀，但是存在主觀性。因此，本文從擬合優(yōu)度的角度來介紹比例風(fēng)險假設(shè)的檢驗方法。

1 原理闡述

Grambsch和Therneau證明了縮放Schoenfeld殘差在診斷比例風(fēng)險假設(shè)方面的作用［2］。在比例風(fēng)險假設(shè)成立的條件下，每個模型中的協(xié)變量的縮放Schoenfeld殘差與生存時間的函數(shù)斜率為0?；诖?，可以構(gòu)建檢驗統(tǒng)計量。

對于單個預(yù)測變量（即自變量）的檢驗統(tǒng)計量見（1）：

式（1）中，rs是縮放Schoenfeld殘差，g(t)是檢驗之前預(yù)先指定的時間函數(shù)形式，δi是事件的指示變量，Δ是事件總數(shù)，是目標預(yù)測變量的參數(shù)估計值的方差。該統(tǒng)計量近似服從自由度為1的χ2分布。

對于p個預(yù)測變量的總體檢驗的統(tǒng)計量見式（2）：

式（2）中，ri是目標預(yù)測變量未縮放的Schoenfeld殘差向量。該統(tǒng)計量服從自由度為p的漸近χ2分布［3］。

除了上述檢驗方法外，也可以在SAS軟件RHREG過程的MODEL語句中加入?yún)f(xié)變量與時間函數(shù)的交互項，使用Wald檢驗來判斷協(xié)變量是否符合比例風(fēng)險假設(shè)。

2 SAS軟件實現(xiàn)

2.1 示例數(shù)據(jù)說明

本文使用的數(shù)據(jù)集Myeloma來源于多發(fā)性骨髓瘤的研究數(shù)據(jù)［4］。原始數(shù)據(jù)、變量說明和SAS數(shù)據(jù)集等內(nèi)容參見前文，此處從略。

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

以下程序演示了如何通過SAS軟件的PHREG過程進行Schoenfeld殘差檢驗以及如何向模型中加入自變量與時間函數(shù)交互項檢驗比例風(fēng)險假設(shè)。

【程序說明】①是基于加權(quán)Schoenfeld殘差進行檢驗的SAS程序。Proc phreg表示調(diào)用PHREG過程，Data選項指定要分析的數(shù)據(jù)集為myeloma，后面的ZPH選項（注意：此選項在SAS 9.3中無效）要求進行比例風(fēng)險假設(shè)檢驗，使用的是加權(quán)Schoenfeld殘差；ZPH選項中的transform指定時間的變換形式，本例使用的是對數(shù)變換；Class語句指定platelet和frac是分類變量；Model語句指定要擬合的模型變量；等號左邊time*vstatus（0）中的time是生存時間，vstatus指定截尾指示變量，括號中的0表示vstatus=0為截尾；等號右邊指定要分析的自變量，Selection=s表示使用逐步選擇進行變量篩選。②是在Cox回歸模型中加入自變量與時間函數(shù)的交互項來進行檢驗，語句和①基本一致。不同之處是在MODEL語句中加入了自變量與時間函數(shù)的交互項，即Log-BUNt、HGBt、Aget、LogWBCt、LogPBMt、Proteint 和Scalct。MODEL語句后面指定交互項的具體計算式。

輸出結(jié)果分別見表1、表2、表3。

表1 Cox比例風(fēng)險回歸模型中參數(shù)的最大似然估計結(jié)果

表2 ZPH選項給出的比例風(fēng)險假設(shè)檢驗結(jié)果

表3 Cox非比例風(fēng)險回歸模型中參數(shù)的最大似然估計結(jié)果

表1是使用逐步回歸進行變量篩選后的結(jié)果?？梢钥吹阶罱K模型中保留的變量為LogBUN和HGB（P＜0.05）；危險比分別為5.336、0.888，它們是基于Cox比例風(fēng)險回歸模型計算得到的。

表2是ZPH選項給出的比例風(fēng)險假設(shè)檢驗結(jié)果?？梢钥闯瞿Ｐ偷娜謾z驗結(jié)果：卡方統(tǒng)計量為6.6966（P＜0.05），因此，模型中至少存在一個變量不滿足比例風(fēng)險假設(shè)。而LogBUN的卡方統(tǒng)計量為6.6837（P＜0.05），即LogBUN不滿足比例風(fēng)險假設(shè)。

表3為最終的最大似然估計結(jié)果，可以看出，LogBUN和LogBUN與時間對數(shù)的交互項均保留在模型中（P＜0.05），因此，可以得出LogBUN不滿足比例風(fēng)險假設(shè)。

綜合以上兩種檢驗結(jié)果，可以看出，LogBUN不滿足比例風(fēng)險假設(shè)。

3 討論與小結(jié)

3.1 討論

殘差可以體現(xiàn)構(gòu)建模型的擬合程度。由于生存分析資料的特殊性，其殘差的計算相比一般的參數(shù)回歸分析方法更為復(fù)雜。在生存資料中，常見的殘差包括Martingale殘差、Deviance殘差、Score殘差以及Schoenfeld殘差。其中，Schoenfeld殘差在檢驗比例風(fēng)險假設(shè)時更為常用。因此，本文介紹了一種基于Schoenfeld殘差進行比例風(fēng)險假設(shè)檢驗的方法。該方法克服了圖示法中只能直觀判斷某變量是否滿足比例風(fēng)險假設(shè)的缺點，以構(gòu)建統(tǒng)計量的形式給出了一種客觀的數(shù)值度量方式。這種方法在本文中是通過使用SAS軟件PHREG過程ZPH選項實現(xiàn)的，但這個選項在SAS/STAT 13.1（即SAS 9.4）以后版本才可用。使用SAS/STAT 13.1以前版本的讀者可考慮使用Chen［3］編寫的PH_score_test宏程序來實現(xiàn)，或者考慮構(gòu)建協(xié)變量與時間函數(shù)的交互項，將之放入MODEL語句中進行擬合的方式來檢驗比例風(fēng)險假設(shè)是否成立。但這種方法存在一個局限，即無法確定合理的時間函數(shù)形式。

在明確存在時間依賴型協(xié)變量以后，可以考慮擴展的Cox回歸模型來進行后續(xù)的模型構(gòu)建，也可以選用RMST回歸［5］來進行擬合。

3.2 小結(jié)

本文主要介紹了使用殘差檢驗法和構(gòu)建協(xié)變量與時間函數(shù)的交互項、從客觀角度進行檢驗的方法。這兩種方法在本文示例中檢驗的結(jié)果一致。但考慮到時間函數(shù)形式無固定標準，讀者需要根據(jù)實際資料靈活選擇檢驗方法。