姚婷婷 ,劉媛媛 ,李長平 ,2*,胡良平
(1.天津醫(yī)科大學公共衛(wèi)生學院,天津 300070;2.世界中醫(yī)藥學會聯(lián)合會臨床科研統(tǒng)計學專業(yè)委員會,北京 100029;3.軍事科學院研究生院,北京 100850
在對生存資料進行分析時,若同時分析眾多因素對生存結(jié)局和生存時間的影響,需要采用多因素分析方法,而傳統(tǒng)的多因素分析方法并不適用,不能同時處理生存結(jié)局和生存時間,也不能充分利用刪失時間所提供的不完全信息。適用于生存數(shù)據(jù)的多因素生存分析方法主要有參數(shù)回歸模型和半?yún)?shù)回歸模型兩類,參數(shù)法需要以特定的分布為基礎建立回歸模型,應用有其局限性,而半?yún)?shù)法的假定相對較少或較弱,特別是Cox比例風險回歸模型(Cox's proportional hazards regression model),不要求生存資料滿足特定的分布類型,是目前對生存資料進行多因素分析最常用的方法之一。
Cox比例風險回歸模型見式(1):
在式(1)中,X1、X2、…、Xp為與生存時間可能有關的自變量(即影響因素),其中的自變量或影響因素可能是定量的或定性的,在整個觀察期內(nèi)不隨時間的變化而變化;h(t)為具有自變量X1、X2、…、Xp的個體在t時刻的風險函數(shù);h0(t)為所有自變量為0時t時刻的風險函數(shù),稱為基準風險函數(shù)(baseline hazard function),是未知的;β1、β2、…、βp為各自變量的偏回歸系數(shù),是一組未知的參數(shù),需要根據(jù)實際的數(shù)據(jù)來估計。
Cox模型不直接考察生存函數(shù)S(t)與自變量的關系,而是利用生存函數(shù)S(t)與風險函數(shù)h(t)的關系,將風險函數(shù)h(t)作為因變量,間接反映自變量與生存函數(shù)S(t)的關系。該模型右側(cè)可分為兩個部分:一部分為h0(t),它沒有明確的定義,分布無明確的假定,為非參數(shù)部分;另一部分是以p個自變量的線性組合為指數(shù)的指數(shù)函數(shù),具有參數(shù)模型形式,其中回歸系數(shù)反映自變量的效應,可通過樣本實際觀測值來估計。所以Cox比例風險回歸模型實為半?yún)?shù)模型(semi-parametric model)[1-6]。
由Cox比例風險回歸模型可知,任意兩個個體風險函數(shù)之比,即風險比(hazard ratio,HR)為:
該比值與h0(t)無關,也與時間t無關,即模型中自變量的效應不隨時間的改變而改變,具有某種特定預后因素向量的患者的死亡風險與具有另一種特定預后因素向量的患者的死亡風險在所有時間點上都保持一個恒定的比例,這種情形被稱為比例風險(proportional hazard)假定,簡稱PH假定。
PH假定的檢驗方法包括圖示法和檢驗法。圖示法是通過觀察散點圖中散點的分布或趨勢來判斷生存時間是否滿足或近似滿足PH假定,主要有COX-KM生存曲線、基于累計風險函數(shù)的圖示法、Schoenfeld殘差圖、Score殘差圖;檢驗法是通過構造滿足PH假設下服從某一已知分布的統(tǒng)計量,基于檢驗統(tǒng)計量的數(shù)值大小和對應的P值來判斷是否滿足或近似滿足PH假定,主要有時間協(xié)變量法、線性相關檢驗、加權殘差Score檢驗、三次樣條函數(shù)法。因篇幅所限,在此僅介紹COX-KM生存曲線和基于累計風險函數(shù)的圖示法,對其他方法感興趣的讀者可參閱文獻[4-6]。
Cox比例風險回歸模型中偏回歸系數(shù)βi的實際意義是:設δi代表第i個自變量在兩個不同個體身上取值差量的絕對值,在其他自變量取值不變的條件下,變量δi每增加一個單位所引起的風險比的自然對數(shù),即lnHRi=βi。
當βi> 0時,HRi> 1,說明Xi增加時,風險函數(shù)增加,Xi為危險因素(其真正含義是:此類因素取高水平相對于取低水平風險增大);當βi<0時,HRi< 1,說明Xi增加時,風險函數(shù)下降,Xi為保護因素(其真正含義是:此類因素取高水平相對于取低水平風險減少);當βi=0時,HRi=1,說明Xi增加時,風險函數(shù)不變,Xi為對生存時間無影響的因素。
偏回歸系數(shù)β1、β2、…、βp的估計需借助偏似然函數(shù),用最大似然估計方法得到。偏似然函數(shù)的計算公式見式(3):
對偏似然函數(shù)取對數(shù),得到對數(shù)偏似然函數(shù)lnL,求lnL關于βj(j=1,2,…,p)的一階偏導數(shù)為0的解(通常需要采用非線性迭代算法,此處從略),便可獲得βj的最大似然估計值bj。
假設檢驗方法類似于logistic回歸分析,有似然比檢驗、Wald檢驗和Score檢驗,詳細理論此處從略;使用統(tǒng)計軟件計算時,參數(shù)估計與假設檢驗都可方便地實現(xiàn)。
【例1】30例膀胱腫瘤患者生存資料原始記錄見表1。研究者欲分析影響膀胱腫瘤患者生存時間(月)長短的因素,包括年齡、腫瘤分級、腫瘤大小和是否復發(fā),并根據(jù)影響因素的取值對不同患者的生存情況進行預測。
表1 30例膀胱腫瘤患者生存資料原始記錄
表1記錄了30例膀胱腫瘤患者的年齡、腫瘤分級、腫瘤大小和是否復發(fā)等情況。其中,年齡和生存時間是定量變量,腫瘤分級、腫瘤大小、是否復發(fā)和生存結(jié)局是定性變量。由于存在截尾數(shù)據(jù)、生存時間及生存結(jié)局,且涉及到多個影響因素,所以,此資料為多因素影響下的一元生存資料。具體地說,生存時間為定量結(jié)果變量、生存結(jié)局為定性結(jié)果變量,它們的信息將被整合在一起參與Cox比例風險模型的建模;而其他變量都屬于自變量或影響因素,其中年齡為定量自變量、腫瘤分級為多值有序自變量、腫瘤大小和是否復發(fā)為二值自變量。
創(chuàng)建一個名為“tjfx”的臨時數(shù)據(jù)集,數(shù)據(jù)步程序如下:
【SAS數(shù)據(jù)步程序說明】因篇幅所限,此處僅列出部分觀測。詳細數(shù)據(jù)見表1。
該資料中年齡為定量變量,將年齡轉(zhuǎn)化為二分類變量(<60歲和≥60歲),分別按年齡、腫瘤分級、腫瘤大小和是否復發(fā)這四個變量的水平分組,采用Kaplan-Meier法繪制生存曲線,程序如下:
【SAS程序說明】生存曲線可用lifetest過程繪制,method用于指定計算生存率的方法,PL表示生存率的乘積極限估計法,即Kaplan-Meier法,plots=(s)用于繪制生存曲線;time語句為lifetest過程的必需語句,設置生存時間變量和生存結(jié)局變量,括號內(nèi)為截尾變量的標示值;strata語句用于指定分層變量。
【SAS主要輸出結(jié)果】
圖1 年齡各水平下的生存曲線
圖2 腫瘤分級各水平下的生存曲線
圖3 腫瘤大小各水平下的生存曲線
圖4 是否復發(fā)的生存曲線
以上四圖顯示,除兩年齡組生存曲線在近30月處略有交叉外,其余各圖中曲線均基本無交叉,說明四個變量基本滿足PH假定,可以進行Cox比例風險回歸模型分析。
利用以下SAS程序,擬合Cox比例風險回歸模型,計算每位患者的預后指數(shù)pi及其所對應生存時間的生存率。
【SAS程序說明】phreg過程是實現(xiàn)Cox模型回歸分析的標準過程,其中class語句可以為分類變量設置啞變量,ref=選項用來指定參考水平,這里需注意的是:腫瘤分級為有序多分類變量,不同的腫瘤分級之間并非是嚴格的等距關系,因此也將其轉(zhuǎn)化為啞變量進行量化;model語句是必需語句,等號左邊為生存時間和生存結(jié)局變量(括號內(nèi)為截尾標識),右邊為協(xié)變量(即自變量),其中選項selection=forward|backward|stepwise|none|score用來指定變量篩選的方法,分別表示前進法、后退法、逐步法、全回歸模型和最優(yōu)子集法,sle=和sls=分別指定引入和剔除自變量的顯著性水平,rl用來指定輸出風險比hr的100(1-α)%置信限;程序中output語句創(chuàng)建一個新的SAS數(shù)據(jù)集report,含有為每一個觀測計算的一些統(tǒng)計量,SAS為每一個統(tǒng)計量定義一個關鍵字,如生存率和預后指數(shù)分別用survival和xbeta表示。選項order=data規(guī)定輸出的數(shù)據(jù)集中的觀測順序與輸入數(shù)據(jù)集中的順序一致。
【SAS主要輸出結(jié)果及解釋】
以上是模型的基本信息、分類變量的分類水平信息以及生存結(jié)局信息。
以上是經(jīng)逐步回歸后,用最終保留下來的協(xié)變量建立回歸模型計算出的模型擬合統(tǒng)計量及模型檢驗結(jié)果,結(jié)果表明模型成立,即因變量與自變量之間的關系可以用所建立的回歸方程來表示。
以上是模型的最大似然估計結(jié)果,包括參數(shù)估計值、估計值標準誤、Waldχ2值、P值、風險比HR及其95%置信區(qū)間。由似然估計結(jié)果得出風險函數(shù)的表達式為:
Cox回歸模型計算結(jié)果顯示:腫瘤分級、腫瘤大小和是否復發(fā)為膀胱腫瘤患者生存時間長短的重要影響因素。在腫瘤大小和是否復發(fā)保持不變的情形下,II級腫瘤患者死亡風險是I級腫瘤患者死亡風險的3.824倍(e1.34121),III級腫瘤患者死亡風險是I級腫瘤患者的31.848倍(e3.46098);在腫瘤分級和是否復發(fā)保持不變的情形下,腫瘤大于或等于3.0 cm者死亡風險是腫瘤小于3.0 cm者死亡風險的2.728倍(e1.00357);在腫瘤分級和腫瘤大小保持不變的情形下,腫瘤復發(fā)者死亡風險是未復發(fā)者死亡風險的2.786倍(e1.02450)。
風險函數(shù)表達式右邊變量的線性組合部分與風險函數(shù)成正比,其取值越大,風險越大,反映了一個個體的預后情況,稱為預后指數(shù)(prognostic index,PI),其值越大,則風險函數(shù)h(t)的取值就越大,預后越差。此案例的預后指數(shù)為:
以上是report數(shù)據(jù)集中的內(nèi)容,包括患者的基本信息、預后指數(shù)及其對應生存時間的生存率(由于篇幅限制,此處僅給出前5例患者的信息)。對于1號患者,腫瘤分級I級,腫瘤小于3.0 cm,未復發(fā),其預后指數(shù)為0,59個月的生存率為24.12%。
相對于非參數(shù)和參數(shù)回歸模型而言,半?yún)?shù)回歸模型兼有兩者的優(yōu)點,且不要求資料服從特定的概率分布,具有靈活性和穩(wěn)健性,而且現(xiàn)如今還沒有非常精準的方法判定待分析的生存資料中的生存時間服從何種分布,使得Cox比例風險回歸模型在醫(yī)學隨訪研究中得到廣泛的應用。雖然Cox比例風險回歸模型的適用范圍廣,但它依賴于嚴格的PH假定,若資料不滿足PH假定,則會較大程度上影響計算的結(jié)果,甚至導致錯誤的結(jié)論。因此,在統(tǒng)計分析前,對PH假定的檢驗是重要且必須的。
本文比較詳細地介紹了Cox比例風險回歸模型、構建模型需要滿足的PH假定及其判定方法,并通過一個實例,基于SAS軟件實現(xiàn)了Cox比例風險回歸模型的構建、校正混雜因素后的組間比較以及對每個個體進行預后指數(shù)和生存率的預測。在進行Cox回歸建模時,需注意只有滿足PH假定前提下,基于此模型的分析預測才是可靠有效的。其次,還應注意回歸所需樣本含量的問題,經(jīng)驗估算方法是至少需要相當于協(xié)變量個數(shù)10~15倍的陽性結(jié)局事件數(shù),或者根據(jù)Hsieh和Lavori給出的樣本量估算公式進行計算,詳細信息可參閱文獻[7]。