生存數(shù)據(jù)分析

1、第十八章生存分析和COX回歸,上海第二醫(yī)科大學 生物統(tǒng)計教研室,第一節(jié) 基本概念,在醫(yī)學,生物學研究中,常用到生存分析 (Survival Analysis)方法。例如對于腫瘤等疾病的療效及預后的考核,通常不用治愈率，有效率等表示，而用將來復發(fā)或死亡的時間長短表示,也即生存期來表示。 所謂生存期（survival time）是指從某個標準時刻(如發(fā)病,確診,開始治療或進行手術的時間)算起至死亡或復發(fā)為止的時間。,生存期不同于一般指標的二個特點： 1.有截尾數(shù)據(jù)（censored data) 隨訪中未能知道病人的確切生存時間，只知道病人的生存時間大于某時間。 （1）病人失訪或因其他原因而死亡-失

2、訪 （2）到了研究的終止期病人尚未死亡-終訪 截尾數(shù)據(jù)可記為t+,如： 4+ = 生存時間大于4年。 雖然截尾數(shù)據(jù)提供的信息是不完全的，但不能刪去，因為這不僅損失了資料，而且會造成偏性。,2. 生存期的資料一般不服從正態(tài)分布。 由于上述原因，常用的統(tǒng)計方法不適用，而要用特殊的統(tǒng)計方法。 生存分析是指對于生存期這一指標進行分析的一系列特殊的統(tǒng)計方法。,生存時間不一定專用于死與活的情況,生存時間(存活時間)可定義為從某種起始事件到達某終點事件所經(jīng)歷的時間跨度。例如急性白血病病人從治療開始到復發(fā)為止之間的緩解期;冠心病病人在兩次發(fā)作之間的時間間隔;已作輸卵管結扎的婦女從施行輸卵管吻合手術后至受孕的時

3、間間隔;在流行病學研究中,從開始接觸危險因素到發(fā)病所經(jīng)歷的時間等都可作為生存時間用作生存分析。 有時還收集一些有關因素(稱為自變量或協(xié)變量), 以分析這些協(xié)變量是否對生存時間有影響,影響的大小，是縮短或延長生存時間。這可以通過Cox回歸進行分析,因此,Cox 回歸可看成帶有協(xié)變量的生存分析。,包括： (1)開始觀察日期，終止觀察日期-生存時間 (2)結局（最終的觀察到的是死亡還是存活） 死于該病-完全數(shù)據(jù) 存活或死于其他原因-截尾數(shù)據(jù) 每個生存期數(shù)據(jù)要用2個變量表示：觀察到的生存時間和是否截尾（如：用1表示截尾，用0表示死亡；4+ 用4，1表示；4用4，0表示）。 (3)協(xié)變量-各種影響生存期

4、長短的因素。,隨訪資料的記錄：,第二節(jié) 描述生存時間分布規(guī)律的函數(shù),一. 生存率(Survival Rate) 又稱為生存概率或生存函數(shù),它表示一個病人的生存時間長于時間t的概率,用S(t) 表示： s(t)=P（Tt） 如5年生存率： s(5)=P（T5） 以時間t為橫坐標,S(t)為縱坐標所作的曲線稱為生存率曲線, 它是一條下降的曲線,下降的坡度越陡,表示生存率越低或生存時間越短,其斜率表示死亡速率。,1.2 概率密度函數(shù) (Probability Density Function) 簡稱為密度函數(shù),記為f(t),其定義為: f(t)=lim (一個病人在區(qū)間(t,t+t)內死亡概率/t)

5、 它表示死亡速率的大小。如以t為橫坐,f(t) 為縱坐標作出的曲線稱為密度曲線,由曲線上可看出不同時間的死亡速率及死亡高峰時間?？v坐標越大,其死亡速率越高,如曲線呈現(xiàn)單調下降,則死亡速率越來越小,如呈現(xiàn)峰值,則為死亡高峰。,1.3 風險函數(shù)(Hazard Function) 用h(t)表示,其定義為: h(t)=lim(在時間t生存的病人死于區(qū)間(t,t)的概率/t) 由于計算h(t)時,用到了生存到時間t,這一條件,故上式極限式中分子部分是一個條件概率?？蓪(t)稱為生存到時間t的病人在時間t的瞬時死亡率或條件死亡速率或年齡別死亡速率。當用t作橫坐標,h(t)為縱坐標所繪的曲線,如遞增,則

6、表示條件死亡速率隨時間而增加,如平行于橫軸,則表示沒有隨時間而加速(或減少)死亡的情況。,風險函數(shù)的不同情況： 常數(shù)， 如：死于飛機失事。 下降， 如：急性損傷。 上升， 如：持續(xù)接觸危險因素。 澡盆樣，如：人的一生。,生存分析目的： （1）估計生存函數(shù)。 （2）比較各組的生存函數(shù)。 （3）研究影響生存期長短的因素。,第三節(jié) 生存率的估計方法,生存率S(t)的估計方法有參數(shù)法和非參數(shù)法。常用非參數(shù)法，非參數(shù)法主要有二個,即,乘積極限法與壽命表法,前者主要用于觀察例數(shù)較少而未分組的生存資料,后者適用于觀察例數(shù)較多而分組的資料,不同的分組壽命表法的計算結果亦會不同,當分組資料中每一個分組區(qū)間中最多

7、只有 1個觀察值時,壽命表法的計算結果與乘積極限法完全相同。,參數(shù)法可求出一個方程表示生存函數(shù)S(t)和時間t的關系，畫出的生存曲線是光滑的下降曲線。 非參數(shù)法只能得到某幾個時間點上的生存函數(shù)，再用直線聯(lián)起來，畫出的生存曲線是呈梯型的。,一. 乘積極限法(Product-Limit Method),簡稱為積限法或PL法,它是由統(tǒng)計學家Kaplan和Meier于1958年首先提出的, 因此又稱為Kaplan-Meier法, 是利用條件概率及概率的乘法原理計算生存率及其標準誤的。 設S(t)表示t年的生存率,s(ti/ti-1)表示活過ti-1年又活過 ti年的條件概率,例如s(1),s(2)分別

8、表示一年,二年的生存率,而s(2/1)表示活過一年者,再活一年的條件概率,據(jù)概率的乘法定律有: S(2)=S(1)S(2/1),一般地有 S(ti)=S(ti-1)S(ti/ti-1),例22.1 用某中藥加化療(中藥組)和化療(對照組)兩種療法治療白血病后, 隨訪記錄各患者的生存時間,不帶+號者表示已死亡,即完全數(shù)據(jù),帶+ 號者表示尚存活,即截尾數(shù)據(jù),試作生存分析。時間單位為月。 中藥組 10,2+,12+,13,18,6+,19+,26,9+,8+,6+,43+,9,4,31,24 對照組 2+,13,7+,11+,6,1,11,3,17,7,資料中藥組積限法計算生存率 時間 狀態(tài) 期初人

9、數(shù) 死亡人數(shù) 條件生存率 累積生 di di/ni(ni-di)累積生存 ti si ni di (ni-di)/ni 存率S(ti)ni(ni-di) 率標準誤 = 2 活 4 死 15 1 0.9333 0.9333 0.004762 0.004762 0.0644 6 活 6 活 8 活 9 死 11 1 0.9090 0.8485 0.009091 0.013853 0.0999 9 活 10 死 9 1 0.8889 0.7542 0.013889 0.027742 0.1256 12 活 13 死 7 1 0.8571 0.6465 0.023810 0.051551 0.1468

10、 18 死 6 1 0.8333 0.5387 0.033333 0.084885 0.1570 19 活 24 死 4 1 0.7500 0.4040 0.083333 0.168218 0.1657 26 死 3 1 0.6667 0.2694 0.166667 0.334885 0.1559 31 死 2 1 0.5000 0.1347 0.500000 0.834885 0.1231 43 活 ,二. 壽命表法(Life Table Method),適用于隨訪的病例數(shù)較多, 將資料按生存期進行分組,在分組的基礎上計算生存率 ,本法也能用于不分組的資料,此時計算結果與積限法相同。,某醫(yī)院

11、1946年1月1日到1951年12月31日收治的126例胃癌病例,生存情況如表22.2,試用壽命表法估計生存率。 表22.2 126例胃癌患者壽命表法估計生存率 時間(年) 期初例數(shù) 死亡例數(shù) 失訪例數(shù) 截尾例數(shù) 有效例數(shù) 條件生存率 累積生存率 di di/ni(ni-di)累積生存 ti ni di ui wi ni S(ti/ti-1) S(ti) ni(ni-di) 率標準誤 = 0- 126 47 4 15 116.5 0.5966 0.5966 5.80510-3 5.80510-3 0.0455 1- 60 5 6 11 51.5 0.9029 0.5386 2.08810-3

12、7.89310-3 0.0479 2- 38 2 0 15 30.5 0.9344 0.5033 2.30110-3 0.0102 0.0508 3- 21 2 2 7 16.5 0.8788 0.4423 8.35910-3 0.0186 0.0602 4- 10 0 0 6 7.0 1.0000 0.4423 0 0.0186 0.0602 5- 4 0 0 4 2.0 1.0000 0.4423 0 0.0186 0.0602 ,壽命表法估計生存率步驟如下: 1.將觀察例數(shù)按時間段(年)0-,1-,2-,劃分,分別計數(shù)期初例數(shù),死亡,失訪, 截尾例數(shù)列入表22.2的1-5列。事實上,從第

13、二個時間段開始,期初人數(shù)ni 系由下式算得: ni=ni-1-di-ui-wi 例如第二行,即時間段1-,有 n2=126-47-4-15=60 2.計算各時間段期初實際觀察例數(shù),(亦稱有效例數(shù))ni ni=ni-ui/2-wi/2 上式表明該時間段期初例數(shù)中的失訪,及截尾例數(shù)只計其半時,即得有效例數(shù)。 如第一行,n1=126-4/2-15/2=116.5 3.分別用(22.5)(22.6)(22.7)式計算條件生存率S(ti/ti-1),累積生存率s(ti)及其標準誤。 計算結果已列于表22.2中,第7,8,11列,表中9,10二列系用于第11列的計算。 例如時間段0-中 S(ti/ti-1

14、)=(116.5-47)/116.5=0.5966 S(ti)=10.5966=0.5966 SE(S(ti)=0.59665.80510-3=0.0455 故一年生存率的估計為0.59660.0455 同樣二年生存率的估計為0.53860.0479 由于壽命表法與積限法的累積生存率及其標準誤的計算公式完全相同,所以,當分組資料中每一個分組區(qū)間中最多只有1個觀察值時,壽命表法就是積限法。,第四節(jié) 生存率的比較,當有兩個或兩個以上的生存分布時,我們常需比較它們是否來自同一生存分布,此時的假設檢驗為: H0:樣本所來自的總體生存分布相同。 H1:樣本所來自的總體生存分布不相同。 可選用的檢驗方法有

15、:Logrank法, 廣義Wilcoxon法,和Cox-Mantel法等。當拒絕H0時,認為幾個生存分布不相同。,當不需要整體比較，而只要比較個別時間點上幾組生存率時可用下面方法： （1）兩個生存率比較 生存率S1和S2，其方差為V1和V2 用卡方檢驗： 2 =（S1-S2）2 / （V1+V2） df=1,（2）兩個以上兩個生存率比較 生存率S1，S2和S3，方差為V1，V2和V3 用卡方檢驗： 權重W1=1/V1, W2=1/V2，W3=1/V3 加權平均生存率： S=（W1*S1+W2*S2+W3*S3） /（W1+W2+W3） 2 =W1*(S1-S)2 +W2*(S2-S)2 +W3

16、*(S3-S)2 df=3-1,3.1 Logrank檢驗(Log Rank Test) 當比較的幾個樣本生存分布,全部為完全數(shù)據(jù)時,本檢驗又稱為Savage檢驗。 Logrank檢驗的計算步驟如下: 1.將兩樣本的生存數(shù)據(jù)混合,由小到大排列,并給以秩次i1, 當截尾數(shù)據(jù)與完全數(shù)據(jù)數(shù)值相同時,截尾數(shù)據(jù)排列在后。并設兩樣本含量分別為m1,m2,總例數(shù)n=m1+m2。 例22.1中藥組與對照組生存數(shù)據(jù)排列結果見表22.3中第1,2列。 2.列出所比較的兩組中任一個組的序號i2(本處選用中藥組),記入表22.3中第3列。 3.列出死亡例的序號i3(見表22.3中第4列)。 4.計算非截尾數(shù)據(jù)(完全數(shù)

17、據(jù))各時間點處于危險狀態(tài)的例數(shù)r,它表示該時刻時還剩下多少例數(shù)。r系由與i3相應的i1值計算而得。 r=n-i1+1 (22.10) 例如與生存期7(月)相應的r值系由 r=26-9+1=18 算得,見表中第5列。,5.對秩次i1作logrank變換,即計算logrank變換值W,其算法為 秩次為i1 的序號為i3非截尾數(shù)據(jù)的W值為 i3 W= 1/rj-1 (22.11) j=1 秩次為i1的截尾數(shù)據(jù),首先判斷它在那二個非截尾數(shù)據(jù)之間,如果它在序號i3與i3+1之間則W為 i3 W= 1/rj (22.12) j=1 特別地,當截尾數(shù)據(jù)在第一個非截尾數(shù)據(jù)之前時,取W=0,幾個截尾數(shù)據(jù)落在同樣

18、序號的非截尾數(shù)據(jù)之間時,它們具有相同的W值。 例如表22.3中第1個數(shù)據(jù)為非截尾(已死亡)則由(22.11)式得 W=1/26-1=-0.96 第2,3數(shù)據(jù)都是截尾數(shù)據(jù)(存活),它處于序號i3=1與2之間據(jù)(22.12)式,有 W=1/26=0.04 第4個數(shù)據(jù)為非截尾,則由(22.11)式 W=1/26+1/23-1=-0.92 余類推(見表22.3第6列),6.計算所指定的組別(本例為中藥組,序號為i2)的logrank變換值之和T T=W i2 上式的連加系在指定的i2范圍內相加。其均數(shù)與方差分別為 E(T)=m1/n W V(T)=m1m2/n(n-1)(W-E(T)/m1)2 (22

19、.14)式中連加系在全部觀察值上完成,m1系指所指定的組別的例數(shù),(本例為中藥組m1=16),n為總例數(shù)。 Z=T-E(T)/V(T) Z服從標準正態(tài)分布,故可由Z0.05=1.96,Z0.01=2.58作出統(tǒng)計推斷。 本例資料有T=3.822,E(T)=0.440210-6, V(T)=3.1755,Z=2.145,故P0.05, 拒絕H0,認為兩種療法生存分布不相同。,3.2 Cox-Mantel檢驗(Cox-Mantel Test),又稱廣義Savage檢驗(Generalized Savage Test),可用于兩個或多個生存分布的比較。仍用例22.1的資料說明本檢驗的計算過程,為敘述

20、方便現(xiàn)將中藥組稱為A組,對照組稱為B組。本檢驗的H0,H1同前。,計算步驟為 1.將兩組生存數(shù)據(jù)混合由小到大排列,當截尾數(shù)據(jù)與非截尾數(shù)據(jù)數(shù)值相同時,截尾數(shù)據(jù)排列在后。并指明各生存數(shù)據(jù)的狀態(tài)(死或活)及所屬組別(見表22.4中1-3列) 2.列出A,B兩組各生存時間上的期初人數(shù)及死亡人數(shù)分別以n1i,d1i,n2i,d 2i表示(見表22.4中第4-7列)。 3.在完全數(shù)據(jù)的相應行中計算合并死亡率Pi Pi=(d1i+d2i)/(n1i+n2i) (22.16) 4.在兩組中任選一組(本處用B組)計算 各生存時間點上的期望死亡人數(shù),它由該組期初人數(shù)乘以合并死亡率而得 E(d2i)=n2iPi (

21、22.17) 參見表22.4中8.9兩列。 5.所指定的組別(本處可B組)死亡人數(shù)的期望值與方差為 E(d2i)=n2iPi (22.18) V(d2i)=n1in2i/(n1i+n2i-1)Pi(1-Pi) (22.19) 計算服從標準正態(tài)分布的統(tǒng)計量Z Z=d2i-E(d2i)/V(d2i) (22.20) 可據(jù)Z0.05=1.96,Z0.01=2.58,作出統(tǒng)計推斷。 代入本例資料有 d2i=7, E(d2i)=3.212284, V(d2i)=1.916190, Z=2.7363 故 P0.01,拒絕H0,認為兩種療法的生存期不相同。,3.3 廣義Wilcoxon檢驗(Generali

22、zed Wilcoxon Test),又稱為Breslow檢驗法(Breslow Test),可用于兩個或兩個以上生存分布的比較,其H0,H1同前。本處仍用例22.1資料說明其計算過程。,1.將兩組生存數(shù)據(jù)混合由小到大排列,當截尾數(shù)據(jù)與完全數(shù)據(jù)數(shù)值相同時,截尾數(shù)據(jù)排列在后。并寫出每個生存數(shù)據(jù)的狀態(tài)(死或活)及所屬組別(見表22.5中第1-3列)。 2.用積限估計法對兩組合并資料估計生存率(即累積生存率)S(ti),列于表中第4列 3.計算各生存時間點的計分值Ui。 Ui= S(ti-1)+S(ti)-1 觀察值為完全數(shù)據(jù) S(ti)-1 觀察值為截尾數(shù)據(jù) (22.21) 其中S(0)=1 例如

23、第一個時間點Ti=1是完全數(shù)據(jù)(死亡),故U1=1+0.9615-1=0.9615, 第2,3個時間點為截尾數(shù)據(jù) U2=U3=0.9615-1=-0.0385 第4個時間點為完全數(shù)據(jù) U4=0.9615+0.9197-1=0.8812 余類推。 4.計算任一組的計分值之和的絕對值,本處選擇A組,其計分值已記第6列中, T=Ui(A) (22.22) T的期望值為0,方差為 V(T)=m1m2U2/(m1+m2)(m1+m2-1) (22.23) (22.23)式中U2系指全部生存時間點的U值平方和, Z=T/V(T) (22.24) Z服從標準正態(tài)分布,故可據(jù)Z0.05=1.96,Z0.01=

24、2.58作出統(tǒng)計推斷。 本例資料有T=2.8712,U2=6.6559,V(T)=1.6384,Z=2.243,P0.05,拒絕H0, 認為兩種治療方法的生存期不相同。,第五節(jié) 估計和比較生存函數(shù)的SAS程序 用LIFETEST過程,第六節(jié) COX回歸,COX回歸用于研究各種因素（稱為協(xié)變量，或伴隨變量等）對于生存期長短的關系，進行多因素分析。 h(t,x)=h0(t)exp(1x1 + 2x2 + mxm ） X1，X2，Xm是協(xié)變量 1 ，2，m是回歸系數(shù)，由樣本估計而得。 I 0表示該協(xié)變量是危險因素，越大使生存時間越短 I 0表示該協(xié)變量是保護因素，越大使生存時間越長,h(t,x)=h0(t)exp(1x1 + 2x2 + mxm ） h0(t)為基礎風險函數(shù)，它是全部協(xié)變量X1，X2，Xm都為0或標準狀態(tài)下的風險函數(shù)，一般是未知的。 h(t,x)表示當各協(xié)變量值X固定時的風險函數(shù)，它和h0(t)成比例，所以該模型又稱為比例風險模型（proportional hazard model) COX回歸模型不用于估計生存率，主要用于因素分析。,COX回歸的應用： 和LOGISTIC回歸相似 （1）因素分析 分析哪些因素（協(xié)變量）對生存期的長短有顯著作用。 對各偏回歸系數(shù)作顯著