? Cox回歸模型(Cox regression model)是生存分析最常用的一種方法，由于其原理的簡單性、參數(shù)的易解釋性，且對生存時間無分布形式的要求，自從1972年由D. R. Cox提出以來，被廣泛應用于“時間-事件”數(shù)據(jù)(time-to-event data)，也是生存資料多因素分析的基本方法。其基礎模型要求生存資料滿足比例風險(Proportional Hazard)假定，故又稱為Cox比例風險回歸模型(Cox proportional hazards regression model)。

目錄：

一、什么是PH假定？

二、PH假定的檢驗方法

（一）KM曲線法

（二）?Schoenfeld殘差圖法

（三）累計鞅殘差圖-Supremum檢驗法

（四）構建協(xié)變量與時間的交互項

三、總結

一、什么是PH假定

? 首先，我們要了解Cox風險比例回歸模型的一般表達式：

其中，

t：生存時間；

X：協(xié)變量X1、X2、…、Xm向量；

h(t, X)：具有m個協(xié)變量的個體在t時刻的風險函數(shù)，即生存時間已經(jīng)達到t的個體在t時刻的瞬時風險率；

h0(t)：基準風險函數(shù)，指t時刻所有協(xié)變量取值全為0時的風險函數(shù)。由于在模型估計過程中其可以被消去，因此允許其未知，也可以不對其分布做任何假定；

h(t, X)/h0(t)為風險比HR(Hazard Ratio, HR)，其對數(shù)與協(xié)變量之間呈線性關系，

β1、β2、…、βm：模型中的待估參數(shù)，為各協(xié)變量的偏回歸系數(shù)，含義為其他協(xié)變量不變的情況下，協(xié)變量Xm每變化一個單位，其對數(shù)風險比的平均改變量。

比如，對于只含有一個二分類變量X1的情形，X1=1表示暴露于某因素，X1=0表示未暴露于某因素，此時模型可表達：

? 則暴露的相對風險比為：

? 顯然，不論基線風險如何，在基線以后的任何時間點上，分別在影響因素的“暴露水平”與“非暴露水平”條件下的發(fā)生事件的風險比是恒定的。換句話說：影響因素的“暴露水平”與“非暴露水平”條件下的發(fā)生事件的風險比在整個隨訪期間是恒定的，與基準風險函數(shù)h0(t)無關，也與時間t無關，即模型中自變量的效應值不隨時間而改變，稱為比例風險(Proportional Hazard)假定，簡稱PH假定。

二、 PH假定的檢驗

包括圖示法和假設檢驗法兩類，常用方法包括：

（一）KM曲線法：針對協(xié)變量為定性變量

（二）Schoenfeld殘差與log(t)的關系圖法：協(xié)變量為定量變量或定性變量

（三）累計鞅殘差圖-Supremum 檢驗法：協(xié)變量為定量變量或分類變量

（四）構建協(xié)變量與時間的交互項

?

（一）KM曲線法—適用定性協(xié)變量的PH假定的檢驗

通過繪制定性變量各水平組的Kaplan-Meier 生存率曲線圖(簡稱KM曲線)，如果不同水平組得KM曲線存在交叉，則表示該定性協(xié)變量不滿足比例風險假定。另外，也可以通過繪制ln{-ln[S(t)]}與生存時間t或對數(shù)生存時間log(t)關系圖，如果各組曲線明顯不平行，說明該定性自變量不符合比例風險假設。KM曲線法是一種比較粗糙、主觀的方法，而且即使不相交的曲線也可能不滿足PH假定，因此一般用于PH假定的初步分析。

SAS代碼：

proc lifetest data=Myeloma plots=(survival(cb=hw test atrisk(outside maxlen=13)) LLs);

??????? time Time * VStatus(0);

??????? strata platelet;

run;

?

關鍵代碼：

plots=(survival(cb=hw test atrisk(outside maxlen=13))LLs)

其中：

cb：顯示生存函數(shù)的置信帶（即不同時刻的置信區(qū)間）。默認值為CB=HW，即顯示Hall-Wellner 置信帶。

Test：顯示 STRATA 語句中指定的log-rank檢驗的P值。

LLs：LOGLOGS的簡寫。繪制估計幸存者函數(shù)的負對數(shù)對數(shù)與時間對數(shù)的關系圖。

?

結果：

兩組KM曲線和關系圖不存在交叉基本平行，可以認為滿足PH假定。

（二）Schoenfeld殘差圖法—定性、定量自變量通用

? 該方法通過繪制Y軸為每個協(xié)變量的縮放Schoenfeld殘差，X軸為對數(shù)生存時間log(t)的Schoenfeld殘差圖實現(xiàn)。若滿足等比例風險假設，則該協(xié)變量下的Schoenfeld殘差隨時間圍繞0上下波動，且擬合曲線與X軸基本平行，其斜率接近0，即Schoenfeld殘差與時間不存在相關性。同時，zph檢驗通過考察Schoenfeld殘差與時間的相關系數(shù)是否為0也進一步說明了等比例風險假設。

該方法應用廣泛，對定性和定量資料都適用。

SAS代碼：

PROC PHREG data=Myeloma zph(global transform=log) ;?????? CLASS? platelet frac;

??????? MODEL time*vstatus(0)=LogBUN HGB Platelet Age LogWBC Frac LogPBM Protein Scalc/ selection=s;

RUN;

?

關鍵代碼：

zph(global transform=log)

其中：

zph：請求基于加權殘差的診斷，以檢查比例風險假設。

global：計算全局相關性測試。

結果：

首先看zph檢驗結果，最下面一行全局假設(_Global_）的P=0.0351，P<0.05說明模型中有變量不滿足PH假定，改變量為LogBun。

變量LogBun的Schoenfeld殘差圖顯示，擬合曲線和X軸不平行，可以認為不滿足PH假定。

變量HGB的Schoenfeld殘差圖也顯示，擬合曲線和X軸基本平行，滿足PH假定。

（三）累計鞅殘差圖-Supremum檢驗法——定性、定量協(xié)變量通用

該方法通過繪制各變量的對數(shù)生存時間-累積鞅殘差圖，同時模擬若干條模擬數(shù)據(jù)的路徑圖，若實際路徑圖在模擬路徑圖的范圍內(nèi)，則可認為滿足等比例風險假定。此外，還可以通過“Kolmogorov-Type Supremum”檢驗（簡稱Supremum檢驗）進行等比例風險假定進行假設檢驗。Supremum檢驗通過對若干條（通常為1000條）路徑與實際路徑進行比較得出P值，若P≤α則認為不滿足等比例風險假定。

SAS代碼：

PROC PHREG data=Myeloma ;

??????? CLASS? platelet frac;

??????? MODEL time*vstatus(0)=LogBUN HGB Platelet Age LogWBC Frac LogPBM Protein Scalc/ selection=s;

??????? Assess var=(LogBUN) ph /seed=222 resample=1000;

???????

RUN;

關鍵代碼：

Assess ph /seed=222 resample=1000;

ASSESS：要求進行Cox回歸模型的充分性檢驗。通過這個語句，可檢驗一個或多個協(xié)變量的函數(shù)形式。PH選項要求進行比例風險假設。

Resample 選項即要求 SAS進行1000次模擬計算，并根據(jù)模擬計算結果進行Supremum 檢驗。

Seed：指定抽樣隨機種子數(shù)。

?

結果：

圖為對數(shù)生存時間-累積鞅殘差路徑圖及Supremum檢驗示例，實線為實際累積鞅殘差路徑圖，虛線為模擬數(shù)據(jù)的路徑圖，左下方Pr為Supremum檢驗所得P值。

變量HGB仍然滿足PH假定(P=0.5510)，然而logBUN不滿足PH假定(P=0.0490)。

（四）構建協(xié)變量與時間的交互項

將欲考察協(xié)變量和對數(shù)時間log(t)的以乘積交互項的形式納入模型，若交互項的偏回歸系數(shù)有統(tǒng)計學意義，則不滿足PH假定。

PROC PHREG data=Myeloma ;

??????? CLASS? platelet frac;

??????? MODEL time*vstatus(0)=LogBUN LogBUNt HGB Platelet Age LogWBC Frac LogPBM Protein Scalc/ selection=s;

??????? LogBUNt=LogBUN*log(time);

???????

RUN;

三、總結

推薦Schoenfeld殘差圖法和累計鞅殘差圖-Supremum檢驗法，其通過圖示法與假設檢驗法相結合的形式使結果更為可靠，而且適用范圍廣泛，適用于定量或定性協(xié)變量，。

?

下期：PH假定不成立，我們該怎么解決？

R語言如何實現(xiàn)PH風險假設？

?

?

參考文獻：

[1]Cox D R. Regression models and life‐tables[J]. Journal of the Royal Statistical Society: Series B (Methodological), 1972, 34(2): 187-202.

[2]Grambsch P M, Therneau T M. Proportional hazards tests and diagnostics based on weighted residuals[J]. Biometrika, 1994, 81(3): 515-526.

[3]Lin D Y, Wei L J, Ying Z. Checking the Cox model with cumulative sums of martingale-based residuals[J]. Biometrika, 1993, 80(3): 557-572.

[4]Borgan ?, Zhang Y. Using cumulative sums of martingale residuals for model checking in nested case‐control studies[J]. Biometrics, 2015, 71(3): 696-703.

?

文字編輯、程序編輯：天涯二毛君

審稿：老陳

關注微信公眾號，獲取更多相關內(nèi)容！