Cox回歸模型比例風險假定的檢驗方法研究

北京協(xié)和醫(yī)學院中國醫(yī)學科學院國家心血管病中心阜外醫(yī)院醫(yī)學統(tǒng)計部（100037） 嚴若華 李 衛(wèi)

?

Cox回歸模型比例風險假定的檢驗方法研究

北京協(xié)和醫(yī)學院中國醫(yī)學科學院國家心血管病中心阜外醫(yī)院醫(yī)學統(tǒng)計部（100037） 嚴若華 李 衛(wèi)△

生存分析（survival analysis）是一系列處理“事件發(fā)生時間”變量的統(tǒng)計方法總稱［1］，常用于研究疾病的發(fā)生、轉(zhuǎn)歸、痊愈和死亡。Cox比例風險回歸模型（Cox proportional hazards regression model）［2］作為生存分析中最重要的多因素分析方法之一，被廣泛應用于臨床隨訪資料的危險因素篩選及預測。

Cox回歸模型將風險函數(shù)表達為基準風險函數(shù)與相應協(xié)變量函數(shù)的乘積，通過描述不同人群在不同時刻的風險，來探索各危險因素對生存的影響。其基本形式為：

通過以上公式可以看出，該模型的參數(shù)估計不依賴于基準風險函數(shù)的分布類型，即對于不同個體X1和X2，其風險比

與基準風險函數(shù)無關，且不隨時間t發(fā)生變化。這就是Cox回歸模型最基本的比例風險（proportional hazards，PH）假定。此外，Cox回歸模型還要求滿足對數(shù)線性假定，即協(xié)變量與對數(shù)風險函數(shù)間呈線性［3］。

Cox回歸模型雖然使生存數(shù)據(jù)中的多因素分析成為可能，然而，由于它依賴于嚴格的假定條件，若資料無法滿足，則會較大程度地影響結果的解讀，甚至導致錯誤的結論。因此，在統(tǒng)計分析前，對PH假定的檢驗是重要且必要的。目前Cox回歸模型存在一定濫用現(xiàn)象，多數(shù)研究者在使用該方法時忽略了對PH假定的檢驗，影響了結果的真實性和可靠性。本文希望能夠?qū)ΜF(xiàn)有的PH假定的檢驗方法進行歸納總結，以提示讀者合理使用Cox回歸模型，選擇恰當?shù)姆椒z驗數(shù)據(jù)的適用性，建立穩(wěn)定有效的模型。

比例風險（PH）假定的檢驗方法

目前，檢驗Cox回歸模型PH假定的方法主要有圖示法［4-5］和假設檢驗法［6］兩種。圖示法包括：（1）Cox & K-M比較法，（2）累積風險函數(shù)法，（3）Schoenfeld殘差圖法；假設檢驗法包括：（1）時協(xié)變量法，（2）線性相關檢驗法，（3）加權殘差Score法；（4） Omnibus檢驗法。本文將對上述方法逐一做簡要介紹并附優(yōu)劣分析。

1.圖示法

（1）Cox & K-M比較法

即比較Cox回歸模型與其他非參數(shù)方法（如K-M曲線）估計的生存曲線的形態(tài)差異，若二者趨勢基本一致，且無交叉，則認為數(shù)據(jù)滿足PH假定。該方法最早由Cox［2］本人提出，后經(jīng)Harrel和Lee［7］等人拓展，可用于計數(shù)資料、等級資料和計量資料的分析，其中計量資料需先將變量離散化，再比較各組的Cox和K-M曲線結果。

（2）累積風險函數(shù)法

即觀察：①累積風險函數(shù)H（t|X）對于時間t的趨勢圖，若比例恒定；②累積風險函數(shù)H（t| X）對于基準累積風險函數(shù)H0（t）的趨勢圖，若斜率恒定；③對數(shù)累積風險函數(shù)lnH（t|X）對于時間t的趨勢圖，若相互平行；④對數(shù)累積風險函數(shù)差異lnH（t|X1）-lnH（t|X2）對于時間t的趨勢圖，若恒為常數(shù)。以上四者均可認為數(shù)據(jù)滿足PH假定。

對于Cox回歸模型，下列公式等價：

因此當PH假定成立時，累積風險函數(shù)圖中不同組間的比值H（t|X1）/ H（t| X2）=exp（（X1-X2）β）［11］、對數(shù)累積風險函數(shù)圖中不同組間的差異lnH（t|X1）-lnH（t|X2）=（X1-X2）β［12］應為常數(shù)，即顯示為互成比例或平行［13-15］。

此外，對于二值變量（X =0，1），Anderson［16］等提出可以采用H1（t）對H0（t）圖，即H-H圖來檢驗PH假定。由于H1（t）=H（t|X =1）=H0（t）eβ，因此H-H圖應為一條直線，截距為0，斜率為eβ。此方法同樣適用于連續(xù)變量和等級變量［17-18］。

上述方法雖然實現(xiàn)了以一條曲線代替兩條曲線來評價數(shù)據(jù)的適用性，但值得注意的是，在不滿足PH假定的情況下，＾β并不能靠lnH1（t）/ H0（t）進行簡單估計，因為此時lnH1（t）/ H0（t）≠lnh1（t）/ h0（t）。故當該值隨時間波動較大時，只能得出數(shù)據(jù)不服從PH假定的結論，不能以該值的形狀推斷β［5］。

（3）Schoenfeld殘差圖法

基于上述方法均與時間變量相關，Schoenfeld［19］定義了一個不依賴于t的偏殘差，以檢驗Cox回歸模型的PH假定。

令Ri為ti時刻（即第i個個體發(fā)生事件時）的風險集，則該個體的Schoenfeld殘差可表示為ri=XI-E （Xi|Ri），其中

為給定風險集下的條件期望［20］，可由最大偏似然估計（maximum partial likelihood estimate）得到的代入計算。由于為的解［2］，可以證明，當PH假定滿足時故ri對ti圖應在0附近波動。

此后，Grambsch和Therneau［21］對Scheonfeld殘差進行了標度調(diào)整，得到加權Scheonfeld殘差rid·Sβ，其中d為資料中的所有事件數(shù)。ri*的尺度與＾β一致，即ri*對ti圖在＾β附近波動，說明數(shù)據(jù)滿足PH假定。同時，上述兩位學者在另一篇文獻［22］中提到Score殘差，作為Scheonfeld殘差的擴展，還可用于更復雜的情況。

然而，Scheonfeld殘差圖中的散點趨勢難以評價，尤其對于二值變量，圖中僅有兩條近乎平行的趨勢線r1=1-E（1|R1）和r0=-E（0|R0），無法了解其波動情況［5］。LOWESS（locally-weighted scatterplot smoothing）平滑法［23-24］及其置信區(qū)間的估計［21］可能會對殘差圖的可讀性起到較大的幫助，并且將實際的PH檢驗結果從隨機趨勢中分離出來。

2.假設檢驗法

（1）時協(xié)變量（time-by-covariate interactions）法

Cox在建立回歸模型時，曾提出一個構造時協(xié)變量的方法來檢驗PH假定［2］，即在模型中加入一個交互作用項X·f（t），其中f（t）為時間函數(shù)，使模型擴展為：

此時，

對PH假定的檢驗即可轉(zhuǎn)化為對γ=0的檢驗。

此外，該方法對時間函數(shù)的選擇亦是十分關鍵的。線性函數(shù)（f（t）=a + bt）、對數(shù)函數(shù)（f（t）=a + blnt）、指數(shù)函數(shù)（f（t）=a + bexpt）等單調(diào)函數(shù)在以往的文獻中較為常見［2，10，26］，因為固定的參數(shù)a、b便于計算，且隨時間單調(diào)增/減的交互項也更易于解讀。對非單調(diào)時間函數(shù)的研究相對少見，Stablein等［27］曾考慮使用二次函數(shù)（f（t）=at + bt2）進行建模，Hess［5］提出也可采用分段線性函數(shù)擬合時間交互項。20世紀90年代初期，非參數(shù)模型逐漸被提出并使用，一些研究使用樣條函數(shù)估計時協(xié)變量［28-30］，無需確定時間函數(shù)的具體形式，可以避免模型的選擇錯誤導致的結果偏差，提高檢驗效率。但如何選擇合適的節(jié)點位置和數(shù)量將成為新的問題。

（2）線性相關檢驗法

該方法起源于Schoenfeld［19］提出的偏殘差概念，后經(jīng)Harrel和Lee［7］改進，逐漸發(fā)展為一種簡便的PH假定檢驗方法。其原理主要是：Schoenfeld殘差不依賴于時間變量，因此Schoenfeld殘差與生存時間的秩次線性無關。此時只需檢驗相關系數(shù)ρ=0（檢驗統(tǒng)計量，即可證明數(shù)據(jù)滿足PH假定。

線性相關檢驗法也同樣適用于其他形式的殘差，如鞅殘差（martingale residual）［22，31］等。定義Mi=δi-H0（ti）exp（Xβ）為第i個個體的鞅殘差，則Mi可認為是該個體觀察到的事件數(shù)與Cox回歸模型下期望的事件數(shù)之差，即資料中存在但未被預測到的部分。在PH假定成立的條件下，鞅殘差同樣具有與時間t無關的特性，因此可用于相關性檢驗；此外，鞅殘差還可用于對數(shù)線性假定的識別［32］。

（3）加權殘差Score法

此方法類似于上述兩種方法的綜合。令β（t）=β+ f（t）γ為時間相關的參數(shù)，即

則可以證明［21］，上述方程的Schoenfeld殘差，其期望約為Eri≈Sβ（ti）f（ti）γ，對β（t）=β的檢驗等價于對Schoenfeld殘差的廣義最小二乘檢驗。

其中

由此可以構造一個自由度為p的漸進χ2統(tǒng)計量

若能夠得出β（t）與時間無關，則表明數(shù)據(jù)滿足PH假定。

（4）Omnibus檢驗法

一些研究者［7，14］認為，檢驗Cox回歸模型的PH假定，還可以通過分割生存時間，即預先將時間按順序劃分為幾個互不相交的區(qū)間，在每個區(qū)間內(nèi)分別擬合Cox模型，以比較不同區(qū)間內(nèi)的回歸系數(shù)是否一致。

假設將時間劃分為q個區(qū)間（τ0，τ1），…，（τq-1，τq），其中τ0=0，τq=∞，令

則分段的Cox回歸模型可表達為

該方法最初由Moreau等［33-34］提出，后經(jīng)O′Quigley和Pessione［35］拓展，成為時協(xié)變量法的一個特例［15，23］。此法要求每一區(qū)間內(nèi)的事件數(shù)均衡可比，且不能太少。針對如何確定區(qū)間分割點的位置和數(shù)量，曾被廣泛探討，一些作者［25，35］認為選擇事件發(fā)生時間的分位點是合適的，但該考慮受到刪失機制的影響，若生存數(shù)據(jù)不滿足隨機刪失，則難以給出令人信服的結論；同時區(qū)間的數(shù)量也依賴于樣本量的大小和生存時間的分布。有作者提出，可以采用事后分析來確定分割點［5］，但這樣的做法必然會導致統(tǒng)計效率的下降，影響結果的可信度。

討 論

比例風險假定作為Cox回歸模型非常重要的應用條件之一，應成為模型建立的基本前提。目前，國際上已有一些研究在進行生存分析時提到了PH假定的檢驗［36-38］，然而，絕大多數(shù)使用Cox回歸模型的文章中，仍然缺乏對該假定的考慮。截至2015年7月，在Pubmed中搜索“Cox regression”，共析出文獻24167篇，但其中同時提到“proportional hazards assumption”的文獻僅有29篇，占0.12%；即便放寬搜索條件，在提到“Cox”的121342篇文獻中，也僅有413篇（0.34%）提到了“assumption”，可見在使用Cox回歸模型時，大多數(shù)作者對其假定的檢驗并不關心。另一方面，在29篇提及比例風險假定的文章中，探討方法學的文章達16篇之多，僅13篇臨床研究在建立Cox回歸模型時考慮了PH假定的適用性，其中8篇得到不符合PH假定的結論，并依此對模型進行了調(diào)整，4篇認為變量符合PH假定，還有1篇對PH假定進行了提及，但未進行檢驗。研究者對于PH假定的重視程度還有待提高，目前尚無公認的檢驗方法，也給PH假定的推廣造成了困難。

本文對目前常見的PH假定檢驗方法進行了總結，其中不乏一些簡單直觀的方法，便于使用。本文希望通過描述，讓更多人對PH假定有更深入的理解，使之日后納入Cox回歸模型的評價成為可能。

本文將檢驗PH假定的方法分為圖示法和假設檢驗法兩種。其中圖示法具有清晰、實用等特點，無需復雜的計算，當發(fā)現(xiàn)數(shù)據(jù)可能不滿足PH假定時，還可以指導恰當?shù)臋z驗統(tǒng)計量的選擇。圖示法除了上文中提到的三種外，尚有UCP（updated covariate percentile）圖［25］、Aalen圖［39-41］、Arjas圖［42］等，由于計算過程相對繁瑣，且無相應的軟件模塊可以直接繪制，故在此不做贅述。

相比于圖示法，很大程度上依賴于研究者的主觀判斷，假設檢驗法可以更加客觀準確地給出令人信服的結論。本文探討了四種較為常見的假設檢驗方法，Nicholas［6］的模擬研究表明，時協(xié)變量法、線性相關檢驗法和加權殘差Score法在檢驗PH假定時均有較高的準確率；其他可行的檢驗方法還包括：廣義矩檢驗［43］、線性秩檢驗［44］、Score檢驗［22，45］。假設檢驗法雖然更為正式可靠，但檢驗結果受到樣本量的影響，當樣本量較小時，檢驗的靈敏度降低，樣本量過大時，又可能因為概率原因拒絕原本成立的假設。因此圖示法在已有的研究中更受歡迎，一般只需判斷既定模型近似滿足PH假定即可。各方法的優(yōu)劣分析見表1。

由表1可知，各方法特點不同，適應于不同條件下的應用。在探索性分析階段，建議采用Cox & K-M比較法或累積風險函數(shù)法進行簡單的比例風險描述，若符合PH假定，即可直接建模；當發(fā)現(xiàn)數(shù)據(jù)偏離PH假定時，再考慮計算Schoenfeld殘差進行繪圖或檢驗，以便為改進Cox回歸模型提供幫助。

表1　比例風險檢驗方法的優(yōu)劣分析

參考文獻

［1］Langova K.Survival analysis for clinical studies.Biomed Pap Med Fac Univ Palacky Olomouc Czech Repub，2008，152：303-307.

［2］Cox DR.Regression models and life-tables（with discussion）.Journal of the Royal Statistical Society，Series B：Methodology，1972，34：187-220.

［3］Anderson PK.Statistical models based on counting processes.Berlin：Springer-Verlag，1993.

［4］Hess KR.Graphical methods for assessing violations of the proportional hazards assumption in Cox regression.Statistics in medicine，1995，14：1707-1723.

［5］余紅梅，何大衛(wèi).檢查Cox模型比例風險假定的幾種圖示法.中國衛(wèi)生統(tǒng)計，2000，17：215-218.

［6］Nicholas HN.An empirical comparison of statistical tests for assessing the proportional hazards assumption of Cox′s model.Statistics in Medicine，1997，16：611-626.

［7］Harrel FE，Lee KL.Verifying assumptions of the Cox proportional hazards mode1.Proceedings of the Eleventh Annual SAS Users Group International Conference，1986：823-828.

［8］Breslow NE.Covariance analysis of censored survival data.Biometrics，1974，30：88-89.

［9］Link CL.Confidence intervals for the survival function using Cox′s proportional-hazard model with covariates.Biometrics，1984，40：601-609.

［10］Kalbfleisch JD，Prentice RL.The Statistical analysis of failure time data.New York：Wiley，1980.

［11］Muenz LR.Comparing survival distributions：a review for nonstatisticians II.Cancer Investigation，1983，1：537-545.

［12］Schumacher M.Two-sample tests of Cramer-von Mises-and Kolmogorov-Smirnov-type for randomly censored data.International Statistical Review，1984，52：263-281.

［13］Lagakos SW.The graphical evaluation of explanatory variables in proportional hazard regression models.Biometrika，1981，68：93-98.

［14］Kay R.Goodness of fit methods for the proportional hazards regression model：a review.Revue Epidemiol Sante Publique，1984，32：185-198.

［15］Gray RJ.Some diagnostic methods for Cox regression models through hazard smoothing.Biometrics，1990，46：93-102.

［16］Andersen PK.Testing goodness of fit of Cox′s regression and life model.Biornetrics，1982，38：67-77.

［17］Wei LJ.Testing goodness of fit for proportional hazards model with censored observations.Journal of the American Statistical Association，1984，79：649-652.

［18］Gill RD，Schumacher M.A simple test of the proportional hazards assumption.Biometrika，1987，74：289-300.

［19］Schoenfeld D.Partial residuals for the proportional hazards regression model.Biornetrika，1982，69：239-241.

［20］Cox DR.Partial likehood.Biometrika，1975，62：269-276.

［21］Grambsch PM，Therneau TM.Proportional hazard tests and diagnostics based on weighted residuals.Biometrika，1994，81：515-526.

［22］Therneau TM，Grambsch PM.Martingale-based residuals for survival models.Biomatrika，1990，77：147-160.

［23］Pettitt AN，Daud IB.Investigating time dependence in Cox′s proportional hazards model.Applied Statistics，1990，39：313-329.

［24］Cleveland WS.Robust locally weighted regression and smoothing scatterplots.Journal of the American Statistical Association，1979，74：829-836.

［25］Thaler HT.Nonparametric estimation of the hazard ratio.Journal of the American Statistical Association，1984，79：290-293.

［26］Anderson JA，Senthilselvan A.A two-step regression model for hazard functions.Applied Statistics，1982，31：44-51.

［27］Stablein DM，Carter WH，Novak JW.Analysis of survival data with nonproportional hazard functions.Controlled Clinical Trials，1981，2：149-159.

［28］Gray RJ.Flexible methods for analyzing survival data using splines，with applications to breast cancer prognosis.Journal of the American Statistical Association，1992，87：942-951.

［29］Hess KR.Assessing time-by-covariate interactions in proportional hazards regression model using cubic spline functions.Statistics in Medicine，1994，13：1045-1062.

［30］Heinzl H，Kaider A，Zlabriger G.Assessing interactions binary timedependent covariates with time in Cox proportional hazards regression models using cubic spline functions.Statistics in medicine.1996，15：2589-2601.

［31］Kay R.Proportional hazards regression models and the analysis of censored survival data.Applied Statistics，1977，26：227-237.

［32］Marzec L，Marzec P.Generalized martingale-residual processes for goodness-of-fit inference in Cox′s type regression models.The Annals of Statistics，1997，25：683-714.

［33］Moreau T，O′Quigley J，Mesbah M.A global goodness-of-fit statistic for the proportional hazards model.Applied Statistics，1985，34：212-218.

［34］Moreau T，O′Quigley J，Lellouch J.On D Schoenfeld′s approach for testing the proportional hazards assumption.Biometrics，1986，73：513-515.

［35］O′Quigley J，Pessione F.Score tests for homogeneity of regression effect in the proportional hazards model.Biometrics，1989，45：135-144.

［36］Santori G，F(xiàn)ontana I，Bertocchi M，et al.Application and validation of Cox regression models in a single-center series of double kidney transplantation.Transplantation Proceedings.2010，42：1098-1103.

［37］Abadi A，Saadat S，Yavari P，et al.Comparison of Aalen′s additive and Cox proportional hazards models for breast cancer survival：analysis of population-based data from British Columbia，Canada.Asian Pac J Cancer Prev，2011，12：3113-3116.

［38］Abadi A，Yavari P，Dehghani-Arani M，et al.Cox models survival analysis based on breast cancer treatments.Iran J Cancer Prev，2014，7：124-129.

［39］Aalen OO.Interaction between life history events.Nonparametric analysis for prospective and retrospective data in the presence of censoring.Scandinavian Journal of Statistics，1980，7：161-171.

［40］Mau J.On a graphical method for the detection of time-dependent effects of covariates in survival data.Applied Statistics，1986，35：245-255.

［41］Henderson R，Milner A.Aalen plots under proportional hazards.Applied Statistics，1991，40：401-409.

［42］Arjas E.A graphical method for assessing goodness of fit in Cox′s proportional hazards mdoel.Journal of the American Statistical Association，1988，83：204-212.

［43］Horowitz JL，Neulnann GR.A generalized moments specification test of the proportional hazards model.Journal of the American Statistical Association，1992，87：234-240.

［44］Chappell R.A note on linear rank test and Gill and Schumacher′s tests of proportionality.Biometrika，1992，79：l99-201.

［45］Lin DY，Wei LJ，Ying Z.Checking the Cox model with cumulative sums of martingale-based residuals.Biometrika，1993，80：557-572.

（責任編輯：郭海強）

通信作者：△李衛(wèi)，E-mail：liwei@ mrbc-nccd.com