生存資料回歸模型分析
——生存資料Cox比例風(fēng)險(xiǎn)回歸模型分析

姚婷婷 ，劉媛媛 ，李長(zhǎng)平 ，2*，胡良平

（1.天津醫(yī)科大學(xué)公共衛(wèi)生學(xué)院，天津 300070；2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會(huì)，北京 100029；3.軍事科學(xué)院研究生院，北京 100850

在對(duì)生存資料進(jìn)行分析時(shí)，若同時(shí)分析眾多因素對(duì)生存結(jié)局和生存時(shí)間的影響，需要采用多因素分析方法，而傳統(tǒng)的多因素分析方法并不適用，不能同時(shí)處理生存結(jié)局和生存時(shí)間，也不能充分利用刪失時(shí)間所提供的不完全信息。適用于生存數(shù)據(jù)的多因素生存分析方法主要有參數(shù)回歸模型和半?yún)?shù)回歸模型兩類，參數(shù)法需要以特定的分布為基礎(chǔ)建立回歸模型，應(yīng)用有其局限性，而半?yún)?shù)法的假定相對(duì)較少或較弱，特別是Cox比例風(fēng)險(xiǎn)回歸模型（Cox's proportional hazards regression model），不要求生存資料滿足特定的分布類型，是目前對(duì)生存資料進(jìn)行多因素分析最常用的方法之一。

1 概 述

1.1 Cox比例風(fēng)險(xiǎn)回歸模型簡(jiǎn)介

Cox比例風(fēng)險(xiǎn)回歸模型見式（1）：

在式（1）中，X1、X2、…、Xp為與生存時(shí)間可能有關(guān)的自變量（即影響因素），其中的自變量或影響因素可能是定量的或定性的，在整個(gè)觀察期內(nèi)不隨時(shí)間的變化而變化；h(t)為具有自變量X1、X2、…、Xp的個(gè)體在t時(shí)刻的風(fēng)險(xiǎn)函數(shù)；h0(t)為所有自變量為0時(shí)t時(shí)刻的風(fēng)險(xiǎn)函數(shù)，稱為基準(zhǔn)風(fēng)險(xiǎn)函數(shù)（baseline hazard function），是未知的；β1、β2、…、βp為各自變量的偏回歸系數(shù)，是一組未知的參數(shù)，需要根據(jù)實(shí)際的數(shù)據(jù)來估計(jì)。

Cox模型不直接考察生存函數(shù)S(t)與自變量的關(guān)系，而是利用生存函數(shù)S(t)與風(fēng)險(xiǎn)函數(shù)h(t)的關(guān)系，將風(fēng)險(xiǎn)函數(shù)h(t)作為因變量，間接反映自變量與生存函數(shù)S(t)的關(guān)系。該模型右側(cè)可分為兩個(gè)部分：一部分為h0(t)，它沒有明確的定義，分布無明確的假定，為非參數(shù)部分；另一部分是以p個(gè)自變量的線性組合為指數(shù)的指數(shù)函數(shù)，具有參數(shù)模型形式，其中回歸系數(shù)反映自變量的效應(yīng)，可通過樣本實(shí)際觀測(cè)值來估計(jì)。所以Cox比例風(fēng)險(xiǎn)回歸模型實(shí)為半?yún)?shù)模型（semi-parametric model）［1-6］。

1.2 模型假定及其檢驗(yàn)

由Cox比例風(fēng)險(xiǎn)回歸模型可知，任意兩個(gè)個(gè)體風(fēng)險(xiǎn)函數(shù)之比，即風(fēng)險(xiǎn)比（hazard ratio，HR）為：

該比值與h0(t)無關(guān)，也與時(shí)間t無關(guān)，即模型中自變量的效應(yīng)不隨時(shí)間的改變而改變，具有某種特定預(yù)后因素向量的患者的死亡風(fēng)險(xiǎn)與具有另一種特定預(yù)后因素向量的患者的死亡風(fēng)險(xiǎn)在所有時(shí)間點(diǎn)上都保持一個(gè)恒定的比例，這種情形被稱為比例風(fēng)險(xiǎn)（proportional hazard）假定，簡(jiǎn)稱PH假定。

PH假定的檢驗(yàn)方法包括圖示法和檢驗(yàn)法。圖示法是通過觀察散點(diǎn)圖中散點(diǎn)的分布或趨勢(shì)來判斷生存時(shí)間是否滿足或近似滿足PH假定，主要有COX-KM生存曲線、基于累計(jì)風(fēng)險(xiǎn)函數(shù)的圖示法、Schoenfeld殘差圖、Score殘差圖；檢驗(yàn)法是通過構(gòu)造滿足PH假設(shè)下服從某一已知分布的統(tǒng)計(jì)量，基于檢驗(yàn)統(tǒng)計(jì)量的數(shù)值大小和對(duì)應(yīng)的P值來判斷是否滿足或近似滿足PH假定，主要有時(shí)間協(xié)變量法、線性相關(guān)檢驗(yàn)、加權(quán)殘差Score檢驗(yàn)、三次樣條函數(shù)法。因篇幅所限，在此僅介紹COX-KM生存曲線和基于累計(jì)風(fēng)險(xiǎn)函數(shù)的圖示法，對(duì)其他方法感興趣的讀者可參閱文獻(xiàn)［4-6］。

1.3 參數(shù)解釋

Cox比例風(fēng)險(xiǎn)回歸模型中偏回歸系數(shù)βi的實(shí)際意義是：設(shè)δi代表第i個(gè)自變量在兩個(gè)不同個(gè)體身上取值差量的絕對(duì)值，在其他自變量取值不變的條件下，變量δi每增加一個(gè)單位所引起的風(fēng)險(xiǎn)比的自然對(duì)數(shù)，即lnHRi=βi。

當(dāng)βi＞ 0時(shí)，HRi＞ 1，說明Xi增加時(shí)，風(fēng)險(xiǎn)函數(shù)增加，Xi為危險(xiǎn)因素（其真正含義是：此類因素取高水平相對(duì)于取低水平風(fēng)險(xiǎn)增大）；當(dāng)βi＜0時(shí)，HRi＜ 1，說明Xi增加時(shí)，風(fēng)險(xiǎn)函數(shù)下降，Xi為保護(hù)因素（其真正含義是：此類因素取高水平相對(duì)于取低水平風(fēng)險(xiǎn)減少）；當(dāng)βi=0時(shí)，HRi=1，說明Xi增加時(shí)，風(fēng)險(xiǎn)函數(shù)不變，Xi為對(duì)生存時(shí)間無影響的因素。

1.4 參數(shù)估計(jì)與假設(shè)檢驗(yàn)

偏回歸系數(shù)β1、β2、…、βp的估計(jì)需借助偏似然函數(shù)，用最大似然估計(jì)方法得到。偏似然函數(shù)的計(jì)算公式見式（3）：

對(duì)偏似然函數(shù)取對(duì)數(shù)，得到對(duì)數(shù)偏似然函數(shù)lnL，求lnL關(guān)于βj(j=1，2，…，p)的一階偏導(dǎo)數(shù)為0的解（通常需要采用非線性迭代算法，此處從略），便可獲得βj的最大似然估計(jì)值bj。

假設(shè)檢驗(yàn)方法類似于logistic回歸分析，有似然比檢驗(yàn)、Wald檢驗(yàn)和Score檢驗(yàn)，詳細(xì)理論此處從略；使用統(tǒng)計(jì)軟件計(jì)算時(shí)，參數(shù)估計(jì)與假設(shè)檢驗(yàn)都可方便地實(shí)現(xiàn)。

2 實(shí)例分析

2.1 問題與數(shù)據(jù)結(jié)構(gòu)

【例1】30例膀胱腫瘤患者生存資料原始記錄見表1。研究者欲分析影響膀胱腫瘤患者生存時(shí)間（月）長(zhǎng)短的因素，包括年齡、腫瘤分級(jí)、腫瘤大小和是否復(fù)發(fā)，并根據(jù)影響因素的取值對(duì)不同患者的生存情況進(jìn)行預(yù)測(cè)。

表1 30例膀胱腫瘤患者生存資料原始記錄

表1記錄了30例膀胱腫瘤患者的年齡、腫瘤分級(jí)、腫瘤大小和是否復(fù)發(fā)等情況。其中，年齡和生存時(shí)間是定量變量，腫瘤分級(jí)、腫瘤大小、是否復(fù)發(fā)和生存結(jié)局是定性變量。由于存在截尾數(shù)據(jù)、生存時(shí)間及生存結(jié)局，且涉及到多個(gè)影響因素，所以，此資料為多因素影響下的一元生存資料。具體地說，生存時(shí)間為定量結(jié)果變量、生存結(jié)局為定性結(jié)果變量，它們的信息將被整合在一起參與Cox比例風(fēng)險(xiǎn)模型的建模；而其他變量都屬于自變量或影響因素，其中年齡為定量自變量、腫瘤分級(jí)為多值有序自變量、腫瘤大小和是否復(fù)發(fā)為二值自變量。

2.2 創(chuàng)建SAS數(shù)據(jù)集

創(chuàng)建一個(gè)名為“tjfx”的臨時(shí)數(shù)據(jù)集，數(shù)據(jù)步程序如下：

【SAS數(shù)據(jù)步程序說明】因篇幅所限，此處僅列出部分觀測(cè)。詳細(xì)數(shù)據(jù)見表1。

2.3 利用SAS/STAT中LIFETEST過程繪制生存曲線

該資料中年齡為定量變量，將年齡轉(zhuǎn)化為二分類變量（＜60歲和≥60歲），分別按年齡、腫瘤分級(jí)、腫瘤大小和是否復(fù)發(fā)這四個(gè)變量的水平分組，采用Kaplan-Meier法繪制生存曲線，程序如下：

【SAS程序說明】生存曲線可用lifetest過程繪制，method用于指定計(jì)算生存率的方法，PL表示生存率的乘積極限估計(jì)法，即Kaplan-Meier法，plots=（s）用于繪制生存曲線；time語句為lifetest過程的必需語句，設(shè)置生存時(shí)間變量和生存結(jié)局變量，括號(hào)內(nèi)為截尾變量的標(biāo)示值；strata語句用于指定分層變量。

【SAS主要輸出結(jié)果】

圖1 年齡各水平下的生存曲線

圖2 腫瘤分級(jí)各水平下的生存曲線

圖3 腫瘤大小各水平下的生存曲線

圖4 是否復(fù)發(fā)的生存曲線

以上四圖顯示，除兩年齡組生存曲線在近30月處略有交叉外，其余各圖中曲線均基本無交叉，說明四個(gè)變量基本滿足PH假定，可以進(jìn)行Cox比例風(fēng)險(xiǎn)回歸模型分析。

2.4 利用SAS/STAT中PHREG過程擬合Cox比例風(fēng)險(xiǎn)回歸模型

利用以下SAS程序，擬合Cox比例風(fēng)險(xiǎn)回歸模型，計(jì)算每位患者的預(yù)后指數(shù)pi及其所對(duì)應(yīng)生存時(shí)間的生存率。

【SAS程序說明】phreg過程是實(shí)現(xiàn)Cox模型回歸分析的標(biāo)準(zhǔn)過程，其中class語句可以為分類變量設(shè)置啞變量，ref=選項(xiàng)用來指定參考水平，這里需注意的是：腫瘤分級(jí)為有序多分類變量，不同的腫瘤分級(jí)之間并非是嚴(yán)格的等距關(guān)系，因此也將其轉(zhuǎn)化為啞變量進(jìn)行量化；model語句是必需語句，等號(hào)左邊為生存時(shí)間和生存結(jié)局變量（括號(hào)內(nèi)為截尾標(biāo)識(shí)），右邊為協(xié)變量（即自變量），其中選項(xiàng)selection=forward|backward|stepwise|none|score用來指定變量篩選的方法，分別表示前進(jìn)法、后退法、逐步法、全回歸模型和最優(yōu)子集法，sle=和sls=分別指定引入和剔除自變量的顯著性水平，rl用來指定輸出風(fēng)險(xiǎn)比hr的100(1-α)%置信限；程序中output語句創(chuàng)建一個(gè)新的SAS數(shù)據(jù)集report，含有為每一個(gè)觀測(cè)計(jì)算的一些統(tǒng)計(jì)量，SAS為每一個(gè)統(tǒng)計(jì)量定義一個(gè)關(guān)鍵字，如生存率和預(yù)后指數(shù)分別用survival和xbeta表示。選項(xiàng)order=data規(guī)定輸出的數(shù)據(jù)集中的觀測(cè)順序與輸入數(shù)據(jù)集中的順序一致。

【SAS主要輸出結(jié)果及解釋】

以上是模型的基本信息、分類變量的分類水平信息以及生存結(jié)局信息。

以上是經(jīng)逐步回歸后，用最終保留下來的協(xié)變量建立回歸模型計(jì)算出的模型擬合統(tǒng)計(jì)量及模型檢驗(yàn)結(jié)果，結(jié)果表明模型成立，即因變量與自變量之間的關(guān)系可以用所建立的回歸方程來表示。

以上是模型的最大似然估計(jì)結(jié)果，包括參數(shù)估計(jì)值、估計(jì)值標(biāo)準(zhǔn)誤、Waldχ2值、P值、風(fēng)險(xiǎn)比HR及其95%置信區(qū)間。由似然估計(jì)結(jié)果得出風(fēng)險(xiǎn)函數(shù)的表達(dá)式為：

Cox回歸模型計(jì)算結(jié)果顯示：腫瘤分級(jí)、腫瘤大小和是否復(fù)發(fā)為膀胱腫瘤患者生存時(shí)間長(zhǎng)短的重要影響因素。在腫瘤大小和是否復(fù)發(fā)保持不變的情形下，II級(jí)腫瘤患者死亡風(fēng)險(xiǎn)是I級(jí)腫瘤患者死亡風(fēng)險(xiǎn)的3.824倍（e1.34121），III級(jí)腫瘤患者死亡風(fēng)險(xiǎn)是I級(jí)腫瘤患者的31.848倍（e3.46098）；在腫瘤分級(jí)和是否復(fù)發(fā)保持不變的情形下，腫瘤大于或等于3.0 cm者死亡風(fēng)險(xiǎn)是腫瘤小于3.0 cm者死亡風(fēng)險(xiǎn)的2.728倍（e1.00357）；在腫瘤分級(jí)和腫瘤大小保持不變的情形下，腫瘤復(fù)發(fā)者死亡風(fēng)險(xiǎn)是未復(fù)發(fā)者死亡風(fēng)險(xiǎn)的2.786倍（e1.02450）。

風(fēng)險(xiǎn)函數(shù)表達(dá)式右邊變量的線性組合部分與風(fēng)險(xiǎn)函數(shù)成正比，其取值越大，風(fēng)險(xiǎn)越大，反映了一個(gè)個(gè)體的預(yù)后情況，稱為預(yù)后指數(shù)（prognostic index，PI），其值越大，則風(fēng)險(xiǎn)函數(shù)h(t)的取值就越大，預(yù)后越差。此案例的預(yù)后指數(shù)為：

以上是report數(shù)據(jù)集中的內(nèi)容，包括患者的基本信息、預(yù)后指數(shù)及其對(duì)應(yīng)生存時(shí)間的生存率（由于篇幅限制，此處僅給出前5例患者的信息）。對(duì)于1號(hào)患者，腫瘤分級(jí)I級(jí)，腫瘤小于3.0 cm，未復(fù)發(fā)，其預(yù)后指數(shù)為0，59個(gè)月的生存率為24.12%。

3 討論與小結(jié)

3.1 討論

相對(duì)于非參數(shù)和參數(shù)回歸模型而言，半?yún)?shù)回歸模型兼有兩者的優(yōu)點(diǎn)，且不要求資料服從特定的概率分布，具有靈活性和穩(wěn)健性，而且現(xiàn)如今還沒有非常精準(zhǔn)的方法判定待分析的生存資料中的生存時(shí)間服從何種分布，使得Cox比例風(fēng)險(xiǎn)回歸模型在醫(yī)學(xué)隨訪研究中得到廣泛的應(yīng)用。雖然Cox比例風(fēng)險(xiǎn)回歸模型的適用范圍廣，但它依賴于嚴(yán)格的PH假定，若資料不滿足PH假定，則會(huì)較大程度上影響計(jì)算的結(jié)果，甚至導(dǎo)致錯(cuò)誤的結(jié)論。因此，在統(tǒng)計(jì)分析前，對(duì)PH假定的檢驗(yàn)是重要且必須的。

3.2 小結(jié)

本文比較詳細(xì)地介紹了Cox比例風(fēng)險(xiǎn)回歸模型、構(gòu)建模型需要滿足的PH假定及其判定方法，并通過一個(gè)實(shí)例，基于SAS軟件實(shí)現(xiàn)了Cox比例風(fēng)險(xiǎn)回歸模型的構(gòu)建、校正混雜因素后的組間比較以及對(duì)每個(gè)個(gè)體進(jìn)行預(yù)后指數(shù)和生存率的預(yù)測(cè)。在進(jìn)行Cox回歸建模時(shí)，需注意只有滿足PH假定前提下，基于此模型的分析預(yù)測(cè)才是可靠有效的。其次，還應(yīng)注意回歸所需樣本含量的問題，經(jīng)驗(yàn)估算方法是至少需要相當(dāng)于協(xié)變量個(gè)數(shù)10～15倍的陽性結(jié)局事件數(shù)，或者根據(jù)Hsieh和Lavori給出的樣本量估算公式進(jìn)行計(jì)算，詳細(xì)信息可參閱文獻(xiàn)［7］。