復發事件資料的生存分析在臨床試驗中的應用及SAS實現

上海交通大學醫學院生物統計學教研室(200025) 張莉娜

有些事件可能在一個受試者身上多次發生。每兩次連續事件之間的時間間隔稱為等待時間。這類資料稱為復發事件資料(recurrent event data)。其特點是上次事件的發生對下次事件發生的遲早有很大影響，即對一名受試者來說，各次事件之間是相關的。

本文應用臨床試驗中一個實例分別用總時間模型(total time model)和邊際模型(marginal model)進行復發事件資料的生存分析，并用穩健的夾心方差估計方法對模型的回歸系數進行估計，并給出PHREG的過程實現。

資 料

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

分別在治療前、治療4周、治療12周、、治療24周、治療36周、治療48周、治療52周各個時間點記錄每個受試者的HBV DNA，見表1。

表1 兩組慢性乙型肝炎患者各個時間點的HBV DNA情況

方 法

1．總時間模型

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

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

總時間模型的數據格式整理如表2形式，放在數據集pwp中。

其中loghbv0為治療前HBV DNA的對數值，分層變量enum為陰轉事件序列號，etstart和etstop分別為每次陰轉事件開始觀察的時間和陰轉事件結束觀察的時間，status=1表示陰轉事件發生，status=0表示陰轉事件未發生。

表2 總時間模型數據格式

2．邊際模型

所有受試者，不論個體實際發生了多少次事件，都被包括在所有事件的危險集中，對每次事件都是從隨訪開始計時。對每一受試者，無論他實際發生了多少次終點事件，但在分層賦值中仍有6個值，此值即為隨訪次數，6是發生終點事件的最大可能次數。這一模型是在配合Cox分層模型的過程中加入了一個類似GEE算法的穩健協方差矩陣。此矩陣的結構是:

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

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

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

表3 邊際模型數據格式

其中loghbv0為治療前HBV DNA的對數值，分層變量visit=i表示第i次隨訪，vtstop為每次隨訪結束觀察的時間，status=1表示陰轉事件發生，status=0表示陰轉事件未發生。

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

Lin和Wei提出用最大似然法或偏最大似然法估計回歸系數的參數，并根據個體內的相關性校正方差的估計值(夾心方差估計值)，使標準誤的大小更接近實際情況，是一個穩健的估計方法。

用SAS9．1中的PHREG過程分別以總時間模型和邊際模型來實現復發事件資料的生存分析。

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;

程序解釋:先用總時間模型對復發事件生存資料進行分析，分別估計各變量的公共回歸系數和各變量在各分層的回歸系數，HBV DNA陰轉事件序列號enum為分層變量。

經歷 k(k=1，2，3，4，5，6)次陰轉事件的比例風險函數模型可以寫為:

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)

再用邊際模型對復發事件生存資料進行分析，估計各變量在各分層的回歸系數，隨訪次數visit為分層變量。

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

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)

對于每一個考慮納入模型的協變量需用k個指示變量表示，分別對應于k個時間段。

選項covs(aggregate)表示選用夾心方差估計方法。

結 果

試驗藥(group=2)與對照藥(group=1)出現HBV DNA陰轉的風險比為HR=1.184，P=0.056，即服用試驗藥比服用對照藥，HBV DNA的陰轉風險更大，陰轉時間更短，試驗藥優于對照藥。治療前HBV DNA的對數值每增加一個單位，出現HBV DNA陰轉的風險比為HR=0.848，P＜0.0001，即治療前的HBV DNA值越大，陰轉風險越小，陰轉時間越長。

表4 總時間模型的公共參數估計及假設檢驗

表5 總時間模型的分層參數估計及假設檢驗

g1，g2，g3，g4 的風險比 HR ＞1，g5，g6 的風險比HR＜1，說明經歷1～4次陰轉，服用試驗藥比服用對照藥，HBV DNA的陰轉時間更短;而對于經歷5和6次陰轉這兩個終點事件，服用試驗藥比服用對照藥，HBV DNA的陰轉時間更長。從檢驗水平分析，只有g1具有統計學意義。

由x1～x6的風險比可知，對于前5次HBV DNA陰轉事件，治療前的HBV DNA越大，陰轉風險越小，陰轉時間越長，且隨著陰轉次數的增加，風險比HR越來越接近1。而對于第6次陰轉事件來說治療前的HBV DNA越大，陰轉風險越大，即對于治療前HBV DNA值越大的受試者，在經歷了前5次陰轉事件后，最后一次陰轉的可能性也越大。

由gg1～gg6的風險比都大于1可知，對于每一次隨訪，每一個時間段，服用試驗藥后HBV DNA陰轉的風險都大于服用對照藥，且隨著隨訪時間的增加，風險比HR越來越小，直至接近于1，即隨著隨訪時間的增加，試驗藥和對照藥對于HBV DNA陰轉風險的影響越來越接近。而從檢驗水平分析，前三次隨訪的P值接近0.05，說明服藥后的前期(24周內)試驗組比對照組HBV DNA陰轉的風險大。

由xx1～xx6的風險比都小于1可知，對于每一次隨訪，治療前的HBV DNA越大，在該隨訪時間段的陰轉風險越小，且隨著隨訪時間的增加，風險比HR越來越大，直至接近1，即隨著時間的增加，治療前HBV DNA值的大小對陰轉風險的影響越來越小。

表6 邊際模型的分層參數估計及假設檢驗

討 論

本文中應用的總時間模型和邊際模型都可用來處理復發事件的生存時間資料，所不同的是生存時間的定義以及對數據結構方面的要求不同;并且通過定義協變量在各次復發事件對應不同的取值，實現同一協變量對應各次復發事件的風險比HR不同，以更符合實際情況。

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．蘇霞，錢國華，荀鵬程，等．臨床試驗中多維生存時間資料的分析．中國臨床藥理學與治療學，2006，11(9):1073-1077．

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

7．高惠璇，等編譯．SAS系統SAS/STAT軟件使用手冊．北京:中國統計出版社，1995，4．