宋德勝 ，李長(zhǎng)平 ，2，劉媛媛 ，崔 壯 *，胡良平

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

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

1 原理闡述

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

對(duì)于單個(gè)預(yù)測(cè)變量（即自變量）的檢驗(yàn)統(tǒng)計(jì)量見（1）：

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

對(duì)于p個(gè)預(yù)測(cè)變量的總體檢驗(yàn)的統(tǒng)計(jì)量見式（2）：

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

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

2 SAS軟件實(shí)現(xiàn)

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

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

2.2 比例風(fēng)險(xiǎn)假設(shè)檢驗(yàn)

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

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

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

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

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

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

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

表2是ZPH選項(xiàng)給出的比例風(fēng)險(xiǎn)假設(shè)檢驗(yàn)結(jié)果。可以看出模型的全局檢驗(yàn)結(jié)果：卡方統(tǒng)計(jì)量為6.6966（P＜0.05），因此，模型中至少存在一個(gè)變量不滿足比例風(fēng)險(xiǎn)假設(shè)。而LogBUN的卡方統(tǒng)計(jì)量為6.6837（P＜0.05），即LogBUN不滿足比例風(fēng)險(xiǎn)假設(shè)。

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

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

3 討論與小結(jié)

3.1 討論

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

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

3.2 小結(jié)

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