上海交通大學(xué)醫(yī)學(xué)院生物統(tǒng)計(jì)學(xué)教研室(200025) 張莉娜

有些事件可能在一個(gè)受試者身上多次發(fā)生。每?jī)纱芜B續(xù)事件之間的時(shí)間間隔稱為等待時(shí)間。這類資料稱為復(fù)發(fā)事件資料(recurrent event data)。其特點(diǎn)是上次事件的發(fā)生對(duì)下次事件發(fā)生的遲早有很大影響，即對(duì)一名受試者來(lái)說(shuō)，各次事件之間是相關(guān)的。

本文應(yīng)用臨床試驗(yàn)中一個(gè)實(shí)例分別用總時(shí)間模型(total time model)和邊際模型(marginal model)進(jìn)行復(fù)發(fā)事件資料的生存分析，并用穩(wěn)健的夾心方差估計(jì)方法對(duì)模型的回歸系數(shù)進(jìn)行估計(jì)，并給出PHREG的過(guò)程實(shí)現(xiàn)。

資 料

為了評(píng)價(jià)某新藥治療成人慢性乙型肝炎的有效性，采用多中心、隨機(jī)、雙盲雙模擬、陽(yáng)性藥物平行對(duì)照的Ⅱ期臨床試驗(yàn)，對(duì)238例受試者的HBV DNA進(jìn)行重復(fù)觀測(cè)。以HBV DNA陰轉(zhuǎn)作為終點(diǎn)事件(HBV DNA≤103)，記錄各次HBV DNA陰轉(zhuǎn)的時(shí)間。所有符合試驗(yàn)方案、依從性良好、試驗(yàn)期間未服用禁止用藥、完成CRF的病例納入PP(Per protocol)分析集，由于無(wú)效而提前退出的病例也納入PP分析，共218例受試者進(jìn)入PP分析集。

分別在治療前、治療4周、治療12周、、治療24周、治療36周、治療48周、治療52周各個(gè)時(shí)間點(diǎn)記錄每個(gè)受試者的HBV DNA，見表1。

表1 兩組慢性乙型肝炎患者各個(gè)時(shí)間點(diǎn)的HBV DNA情況

方 法

1．總時(shí)間模型

它把整個(gè)隨訪時(shí)間按終點(diǎn)事件的發(fā)生分為數(shù)個(gè)時(shí)間段，并考慮了終點(diǎn)事件發(fā)生的先后順序。一個(gè)受試者只有在發(fā)生前一次事件后才能進(jìn)入下一個(gè)危險(xiǎn)集，即受試者在k-1個(gè)事件發(fā)生后，才處于發(fā)生第k個(gè)事件的風(fēng)險(xiǎn)中，分層變量采用順序號(hào)賦值是為了保證這種事件發(fā)生的先后順序。所以此模型屬于條件模型(conditional model)。第k次事件的風(fēng)險(xiǎn)函數(shù)為:

其中 β1，β2，…，βp為第 k次終點(diǎn)事件的 p 個(gè)變量的回歸系數(shù);x1k，x2k，…，xpk為發(fā)生第k次終點(diǎn)事件的p個(gè)變量的取值;h0k(t)表示t時(shí)刻發(fā)生第k次終點(diǎn)事件的基線風(fēng)險(xiǎn);hk(t，xk)表示在第t時(shí)刻經(jīng)歷第k次終點(diǎn)事件的風(fēng)險(xiǎn)。此模型由Prentice，Williams和Peterson提出，故又稱為PWP模型。

總時(shí)間模型的數(shù)據(jù)格式整理如表2形式，放在數(shù)據(jù)集pwp中。

其中l(wèi)oghbv0為治療前HBV DNA的對(duì)數(shù)值，分層變量enum為陰轉(zhuǎn)事件序列號(hào)，etstart和etstop分別為每次陰轉(zhuǎn)事件開始觀察的時(shí)間和陰轉(zhuǎn)事件結(jié)束觀察的時(shí)間，status=1表示陰轉(zhuǎn)事件發(fā)生，status=0表示陰轉(zhuǎn)事件未發(fā)生。

表2 總時(shí)間模型數(shù)據(jù)格式

2．邊際模型

所有受試者，不論個(gè)體實(shí)際發(fā)生了多少次事件，都被包括在所有事件的危險(xiǎn)集中，對(duì)每次事件都是從隨訪開始計(jì)時(shí)。對(duì)每一受試者，無(wú)論他實(shí)際發(fā)生了多少次終點(diǎn)事件，但在分層賦值中仍有6個(gè)值，此值即為隨訪次數(shù)，6是發(fā)生終點(diǎn)事件的最大可能次數(shù)。這一模型是在配合Cox分層模型的過(guò)程中加入了一個(gè)類似GEE算法的穩(wěn)健協(xié)方差矩陣。此矩陣的結(jié)構(gòu)是:

Cov(^β)=V(β)［L'L］V(β)

式中V(β)是關(guān)于β的信息矩陣的逆組成的P×P維矩陣，L是n×p維的記分殘差矩陣。此模型由Wei，Lin和Weissfeld提出，故又稱為WLW模型。

邊際模型的數(shù)據(jù)格式整理如表3形式，放在數(shù)據(jù)集wlw中。

表3 邊際模型數(shù)據(jù)格式

其中l(wèi)oghbv0為治療前HBV DNA的對(duì)數(shù)值，分層變量visit=i表示第i次隨訪，vtstop為每次隨訪結(jié)束觀察的時(shí)間，status=1表示陰轉(zhuǎn)事件發(fā)生，status=0表示陰轉(zhuǎn)事件未發(fā)生。

3．夾心方差估計(jì)(sandwich variance estimator)

Lin和Wei提出用最大似然法或偏最大似然法估計(jì)回歸系數(shù)的參數(shù)，并根據(jù)個(gè)體內(nèi)的相關(guān)性校正方差的估計(jì)值(夾心方差估計(jì)值)，使標(biāo)準(zhǔn)誤的大小更接近實(shí)際情況，是一個(gè)穩(wěn)健的估計(jì)方法。

用SAS9．1中的PHREG過(guò)程分別以總時(shí)間模型和邊際模型來(lái)實(shí)現(xiàn)復(fù)發(fā)事件資料的生存分析。

SAS程序如下:

title‘PWP total time model with common regression coefficients’;

proc phreg data=pwp covs(aggregate);model(etstart，etstop)*Status(0)=group loghbv0;id no;strata enum;

run;

title‘PWPtotal time model with stratum-specific regression coefficients’;

proc phreg data=pwp covs(aggregate);

model(etstart，etstop)*Status(0)=g1 g2 g3 g4 g5 g6 x1 x2 x3 x4 x5 x6;

strata Enum;

g1=group*(Enum=1);x1=loghbv0*(Enum=1);

g2=group*(Enum=2);x2=loghbv0*(Enum=2);

g3=group*(Enum=3);x3=loghbv0*(Enum=3);

g4=group*(Enum=4);x4=loghbv0*(Enum=4);

g5=group*(Enum=5);x5=loghbv0*(Enum=5);

g6=group*(Enum=6);x6=loghbv0*(Enum=6);

run;

title‘Wei-Lin-Weissfeld Model’;

proc phreg data=wlw covs(aggregate);

model vtstop*Status(0)=gg1 gg2 gg3 gg4 gg5 gg6 xx1 xx2 xx3 xx4 xx5 xx6;

gg1=group*(visit=1);xx1=loghbv0*(visit=1);

gg2=group*(visit=2);xx2=loghbv0*(visit=2);

gg3=group*(visit=3);xx3=loghbv0*(visit=3);

gg4=group*(visit=4);xx4=loghbv0*(visit=4);

gg5=group*(visit=5);xx5=loghbv0*(visit=5);

gg6=group*(visit=6);xx6=loghbv0*(visit=6);

strata visit;

id no;

run;

程序解釋:先用總時(shí)間模型對(duì)復(fù)發(fā)事件生存資料進(jìn)行分析，分別估計(jì)各變量的公共回歸系數(shù)和各變量在各分層的回歸系數(shù)，HBV DNA陰轉(zhuǎn)事件序列號(hào)enum為分層變量。

經(jīng)歷 k(k=1，2，3，4，5，6)次陰轉(zhuǎn)事件的比例風(fēng)險(xiǎn)函數(shù)模型可以寫為:

H(t，g1，g2，g3，g4，g5，g6，x1，x2，x3，x4，x5，x6)=h0k(t)exp(β11g1+β12g2+ β13g3+ β14g4+ β15g5+ β16g6+ β21x1+β22x2+ β23x3+ β24x4+ β25x5+ β26x6) (k=1，2，3，4，5，6)

再用邊際模型對(duì)復(fù)發(fā)事件生存資料進(jìn)行分析，估計(jì)各變量在各分層的回歸系數(shù)，隨訪次數(shù)visit為分層變量。

第 k(k=1，2，3，4，5，6)次隨訪發(fā)生陰轉(zhuǎn)事件的比例風(fēng)險(xiǎn)模型可以寫為:

H(t，gg1，gg2，gg3，gg4，gg5，gg6，xx1，xx2，xx3，xx4，xx5，xx6)=h0k(t)exxp(β11gg1+ β12gg2+ β13gg3+β14gg4+ β15gg5+ β16gg6+ β21xx1+ β22xx2+ β23xx3+ β24xx4+ β25xx5+ β26xx6)(k=1，2，3，4，5，6)

對(duì)于每一個(gè)考慮納入模型的協(xié)變量需用k個(gè)指示變量表示，分別對(duì)應(yīng)于k個(gè)時(shí)間段。

選項(xiàng)covs(aggregate)表示選用夾心方差估計(jì)方法。

結(jié) 果

試驗(yàn)藥(group=2)與對(duì)照藥(group=1)出現(xiàn)HBV DNA陰轉(zhuǎn)的風(fēng)險(xiǎn)比為HR=1.184，P=0.056，即服用試驗(yàn)藥比服用對(duì)照藥，HBV DNA的陰轉(zhuǎn)風(fēng)險(xiǎn)更大，陰轉(zhuǎn)時(shí)間更短，試驗(yàn)藥優(yōu)于對(duì)照藥。治療前HBV DNA的對(duì)數(shù)值每增加一個(gè)單位，出現(xiàn)HBV DNA陰轉(zhuǎn)的風(fēng)險(xiǎn)比為HR=0.848，P＜0.0001，即治療前的HBV DNA值越大，陰轉(zhuǎn)風(fēng)險(xiǎn)越小，陰轉(zhuǎn)時(shí)間越長(zhǎng)。

表4 總時(shí)間模型的公共參數(shù)估計(jì)及假設(shè)檢驗(yàn)

表5 總時(shí)間模型的分層參數(shù)估計(jì)及假設(shè)檢驗(yàn)

g1，g2，g3，g4 的風(fēng)險(xiǎn)比 HR ＞1，g5，g6 的風(fēng)險(xiǎn)比HR＜1，說(shuō)明經(jīng)歷1～4次陰轉(zhuǎn)，服用試驗(yàn)藥比服用對(duì)照藥，HBV DNA的陰轉(zhuǎn)時(shí)間更短;而對(duì)于經(jīng)歷5和6次陰轉(zhuǎn)這兩個(gè)終點(diǎn)事件，服用試驗(yàn)藥比服用對(duì)照藥，HBV DNA的陰轉(zhuǎn)時(shí)間更長(zhǎng)。從檢驗(yàn)水平分析，只有g(shù)1具有統(tǒng)計(jì)學(xué)意義。

由x1～x6的風(fēng)險(xiǎn)比可知，對(duì)于前5次HBV DNA陰轉(zhuǎn)事件，治療前的HBV DNA越大，陰轉(zhuǎn)風(fēng)險(xiǎn)越小，陰轉(zhuǎn)時(shí)間越長(zhǎng)，且隨著陰轉(zhuǎn)次數(shù)的增加，風(fēng)險(xiǎn)比HR越來(lái)越接近1。而對(duì)于第6次陰轉(zhuǎn)事件來(lái)說(shuō)治療前的HBV DNA越大，陰轉(zhuǎn)風(fēng)險(xiǎn)越大，即對(duì)于治療前HBV DNA值越大的受試者，在經(jīng)歷了前5次陰轉(zhuǎn)事件后，最后一次陰轉(zhuǎn)的可能性也越大。

由gg1～gg6的風(fēng)險(xiǎn)比都大于1可知，對(duì)于每一次隨訪，每一個(gè)時(shí)間段，服用試驗(yàn)藥后HBV DNA陰轉(zhuǎn)的風(fēng)險(xiǎn)都大于服用對(duì)照藥，且隨著隨訪時(shí)間的增加，風(fēng)險(xiǎn)比HR越來(lái)越小，直至接近于1，即隨著隨訪時(shí)間的增加，試驗(yàn)藥和對(duì)照藥對(duì)于HBV DNA陰轉(zhuǎn)風(fēng)險(xiǎn)的影響越來(lái)越接近。而從檢驗(yàn)水平分析，前三次隨訪的P值接近0.05，說(shuō)明服藥后的前期(24周內(nèi))試驗(yàn)組比對(duì)照組HBV DNA陰轉(zhuǎn)的風(fēng)險(xiǎn)大。

由xx1～xx6的風(fēng)險(xiǎn)比都小于1可知，對(duì)于每一次隨訪，治療前的HBV DNA越大，在該隨訪時(shí)間段的陰轉(zhuǎn)風(fēng)險(xiǎn)越小，且隨著隨訪時(shí)間的增加，風(fēng)險(xiǎn)比HR越來(lái)越大，直至接近1，即隨著時(shí)間的增加，治療前HBV DNA值的大小對(duì)陰轉(zhuǎn)風(fēng)險(xiǎn)的影響越來(lái)越小。

表6 邊際模型的分層參數(shù)估計(jì)及假設(shè)檢驗(yàn)

討 論

本文中應(yīng)用的總時(shí)間模型和邊際模型都可用來(lái)處理復(fù)發(fā)事件的生存時(shí)間資料，所不同的是生存時(shí)間的定義以及對(duì)數(shù)據(jù)結(jié)構(gòu)方面的要求不同;并且通過(guò)定義協(xié)變量在各次復(fù)發(fā)事件對(duì)應(yīng)不同的取值，實(shí)現(xiàn)同一協(xié)變量對(duì)應(yīng)各次復(fù)發(fā)事件的風(fēng)險(xiǎn)比HR不同，以更符合實(shí)際情況。

1．Hosmer DW，Lemeshow S．Applied survival analysis:regression modeling of time to event data．John Wiley ＆ Sons，lnc．New Yark，1999:308-316．

2．Shoukri MM，Pause CA．Statistical methods for health sciences．Second Edition．CRC Press LLC，1999:369-375．

3．Kleinbaum DG，Klein M．Survival analysis:a self learning text．Second Edition．Springer Science+Business Media，lnc，2005:331-390．

4．Allison PD．Survival analysis using SAS:a practical guide．Second Edition．Cary，NC，USA:SAS Institute Inc，2010:260-281．

5．蘇霞，錢國(guó)華，荀鵬程，等．臨床試驗(yàn)中多維生存時(shí)間資料的分析．中國(guó)臨床藥理學(xué)與治療學(xué)，2006，11(9):1073-1077．

6．余松林，向惠云，編著．重復(fù)測(cè)量資料分析方法與SAS程序．北京:科學(xué)出版社，2004，3．

7．高惠璇，等編譯．SAS系統(tǒng)SAS/STAT軟件使用手冊(cè)．北京:中國(guó)統(tǒng)計(jì)出版社，1995，4．