基于指數(shù)分布的生存資料盲態(tài)下樣本量再估計方法*

朱　典　王　陵　李嬋娟　蔣志偉　張　威　王治東　吳克堅　李　晨　夏結(jié)來△

【提　要】　目的　針對生存分析樣本量再估計的適應性設計臨床試驗,本文基于指數(shù)分布,提出盲態(tài)下生存資料的樣本量再估計方法。方法　采用蒙特卡洛模擬方法,預設4個期中分析時間點,以參數(shù)初始估計值M10為均數(shù)暫定搜索范圍為M10±0.5M10,在該搜索范圍內(nèi)產(chǎn)生指數(shù)分布隨機數(shù)對期中分析時的截尾生存數(shù)據(jù)進行填補,對填補后的數(shù)據(jù)進行極大似然估計,得到試驗參數(shù)的再估計值若落在搜索范圍外,則需更改搜索范圍以包含并重新進行搜索直至落在搜索范圍內(nèi),此時的即為試驗參數(shù)的再估計值,在此參數(shù)的基礎(chǔ)上重新估計樣本量,并比較4個期中分析點的樣本量再估計結(jié)果,確定1個最合適的期中分析時間點。 結(jié)果　經(jīng)過期中分析調(diào)整后的樣本量接近真實值,且期中分析時間點越向后移,樣本量估計結(jié)果越接近真實值,變異越小。 結(jié)論　建議在入組結(jié)束并完成1/4最短隨訪時間時進行一次期中分析重新估計樣本量,根據(jù)估計結(jié)果考慮是否增加樣本量。

傳統(tǒng)的臨床試驗在試驗開展前,研究者根據(jù)前期研究結(jié)果并結(jié)合相關(guān)文獻設定試驗關(guān)鍵參數(shù),對樣本量進行估計,以期望達到設定的檢驗效能。但是,試驗開始前對參數(shù)的設定可能存在較大偏差,以此為基礎(chǔ)估計的樣本量有可能不足以達到預設的檢驗效能。而在適應性設計中允許在試驗開始之后根據(jù)已得到的部分試驗數(shù)據(jù)進行期中分析,重新估計樣本量,以使試驗達到預設的檢驗效能[1]。

對于不同類型的數(shù)據(jù),需采用不同的樣本量再估計方法。對于正態(tài)資料,Gould[2]等人提出一種基于EM算法的盲態(tài)下樣本量再估計方法,在期中分析時對試驗數(shù)據(jù)的方差進行重新估計。對于二分類變量,Shih[3]等人提出一種基于分層理念的盲態(tài)下樣本量再估計方法,在期中分析時對兩組的事件發(fā)生率重新估計;Friede[4]等人提出一種基于負二項分布的盲態(tài)下樣本量再估計方法,在期中分析時以負二項分布為基礎(chǔ)建模,對總體率進行重新估計。對于生存資料,Togo[5]等人提出一種基于偏似然估計的樣本量再估計方法,在期中分析時對數(shù)據(jù)揭盲,并根據(jù)構(gòu)建的模型對兩組之間的風險比(HR)進行重新估計;Todd[6]等人提出一種基于外推法的盲態(tài)下生存數(shù)據(jù)樣本量再估計方法,假設兩組之間的風險比(HR)不變,試驗進展過程中在多個預設時間點上進行多次期中分析,并根據(jù)多次期中分析結(jié)果求效應的平均差別,推斷試驗結(jié)束時的生存率。

期中分析時對試驗揭盲會導致Ⅰ型錯誤的膨脹,從而影響試驗的完整性。而在盲態(tài)下進行期中分析,對Ⅰ型錯誤的影響可以忽略,無需對檢驗水準進行校正[1]。本研究提出的生存資料樣本量再估計方法是在盲態(tài)下只進行一次期中分析,該方法以極大似然估計為基礎(chǔ),結(jié)合Todd[6]等人的研究方法,假定生存數(shù)據(jù)服從指數(shù)分布,兩組之間的風險比(HR)不變,在期中分析時對其他試驗參數(shù)進行重新估計,提出生存資料樣本量再估計的新方法。

原理與方法

對于服從指數(shù)分布的生存資料,期中分析在盲態(tài)下進行,其分組未知,相當于求解混合指數(shù)分布的參數(shù)問題。EM算法是Dempster[7]等人提出一種對于不完全數(shù)據(jù)進行參數(shù)估計的方法,以極大似然估計為基礎(chǔ)進行迭代,最終得到一個局部最優(yōu)參數(shù)估計的迭代算法。朱利平[8]等人將EM算法應用到電子元件的壽命試驗中。研究表明,樣本的截尾比例越大,EM算法參數(shù)估計的效率越差,且估計結(jié)果往往依賴于初始值和收斂標準的設定[8-9]。在對生存資料進行期中分析時,往往產(chǎn)生大量的截尾數(shù)據(jù)。因此,本研究提出填補數(shù)據(jù)的方法,在期中分析時,用參數(shù)的初始估計值產(chǎn)生模擬數(shù)據(jù),并用模擬數(shù)據(jù)中大于期中分析時間點的數(shù)據(jù)填補截尾數(shù)據(jù),經(jīng)過填補后的數(shù)據(jù)包含一定比例的真實參數(shù)信息和一定比例的填補參數(shù)信息。本研究結(jié)合EM算法的思想,但考慮到其參數(shù)估計結(jié)果依賴于初始值的設定,因此,提出在期中分析時以初始估計值為均值確定搜索范圍(hunting zone,HZ)對實際數(shù)據(jù)進行填補,然后,對填補后的數(shù)據(jù)進行極大似然估計。若初始估計值非常接近參數(shù)真實值,則此時參數(shù)的初始估計值與極大似然估計值也非常接近,則在搜索范圍內(nèi)篩選出參數(shù)的極大似然估計值與初始估計值之差最小的參數(shù)極大似然估計值作為期中分析時的參數(shù)估計值。

極大似然法是混合分布參數(shù)估計的常用方法,期中分析時對服從混合指數(shù)分布的生存數(shù)據(jù)應用極大似然法,重新估計混合指數(shù)分布的參數(shù)進行樣本量再估計。期中分析時的數(shù)據(jù)是含有缺失數(shù)據(jù)的觀測數(shù)據(jù),若用K表示分組(未知),X表示觀測數(shù)據(jù),則包含分組的完整數(shù)據(jù)用Y表示:Y=(X,K)。λ是指數(shù)分布的參數(shù),f(K|X,λ)是給定觀測數(shù)據(jù)X=Xi,λ=λi的條件下缺失數(shù)據(jù)K的條件密度。通過對觀測數(shù)據(jù)X的對數(shù)似然函數(shù)lnL(λ|X)求極大值得到λ的極大似然估計值。期中分析時的生存時間為觀測數(shù)據(jù),用X表示,xi表示每個受試者的生存時間。受試者的分組情況用K表示,Ki=1表示受試者接受對照組治療,Ki=0表示受試者接受試驗組治療,因?qū)φ战M與試驗組的樣本量比例為1:1,所以Ki服從二項分布,π(Ki=1)=π(Ki=0)=0.5,但是,期中分析是在盲態(tài)下進行,因此Ki是不能被觀測到的隨機變量。用λc表示對照組的風險率,λt表示試驗組的風險率,則對照組受試者的生存時間服從指數(shù)分布f1i,其概率密度函數(shù)為:

f1i=λce-λcxi(λc>0),

試驗組受試者的生存時間服從指數(shù)分布f2i,其概率密度為:

f2i=λte-λtxi(λt>0),

則xi服從混合指數(shù)分布fi,其概率密度函數(shù)為:

fi=0.5f1i+0.5f2i.

在混合指數(shù)分布fi中,記θ=(λc,λt),則xi和Ki的聯(lián)合分布為:

g(xi,Ki,θ)=(0.5f1i)Ki(0.5f2i)1-Ki,

因此,Ki在xi給定時的條件分布為:

xi給定時θ的似然函數(shù)為:

對數(shù)似然函數(shù)為:

給定參數(shù)初始估計值θ0(λc0,λt0)時,該對數(shù)似然函數(shù)的期望為:

通過極大化Q(θ,θ0)函數(shù)可得到參數(shù)θ(λc,λt)的極大似然估計值[8]。

在以事件為終點的隨機對照臨床試驗中,研究者根據(jù)以往數(shù)據(jù)估計受試者的生存時間,進而轉(zhuǎn)換為風險率λ進行樣本量計算。指數(shù)分布下風險率λ與中位生存時間M的轉(zhuǎn)換公式為:

本研究中定義隨訪時間為最后一個受試者入組開始至研究結(jié)束所經(jīng)歷的時間。由于入組時間不同,因此最后一個入組受試者的隨訪時間決定了總研究時長,且所有受試者均隨訪至研究結(jié)束時間點,此時樣本量計算公式為(1)～(3)所示[10]:

(1)

(2)

(3)

公式(1)～(3)中,λc和λt分別是對照組和試驗組的風險率,T0是入組時間,T是總研究時長,k是試驗組和對照組的樣本含量比例,nc是對照組樣本量。

采用Monte Carlo模擬構(gòu)建基于指數(shù)分布的生存資料模型,對模擬產(chǎn)生的生存資料進行l(wèi)og-rank檢驗,探討本研究方法對生存資料樣本量再估計的準確性和穩(wěn)定性。本研究所有數(shù)據(jù)均由模擬產(chǎn)生,采用SAS 9.1.3統(tǒng)計分析軟件編寫程序,并進行數(shù)據(jù)分析。

1.參數(shù)設置

對于一個基于log-rank檢驗的雙臂臨床試驗,對照組與試驗組樣本比例為1:1,生存時間服從指數(shù)分布。試驗的真實參數(shù)如下,對照組的中位生存時間為M1=1年,對照組與試驗組的風險比HR為0.65和0.75,試驗入組時間為T0=0.5年,受試者在入組期內(nèi)均勻入組,最短隨訪時間為t=2年,總研究時長為T=T0+t=2.5年,檢驗水準α為雙側(cè)0.05,檢驗效能1-β為0.8,則對照組的風險率λc為ln2/M1,試驗組的中位生存時間為M2=M1/HR,風險率λt為ln2/M2,試驗所需樣本量為N。

在進行期中分析的臨床試驗中,一般進行1～2次期中分析,本模擬研究考慮最大可能性,設定4個期中分析時間點(interim,I),并比較4個期中分析時間點的樣本量再估計結(jié)果,最終選定一個合理的期中分析時間點。4個期中分析時間點分別為I1-I4,I1是入組結(jié)束時(T0),I2是完成1/4隨訪時間時(T0+0.25t),I3是完成1/2隨訪時間時(T0+0.5t),I4是完成3/4隨訪時間時(T0+0.75t)。

2.模擬步驟

(1)根據(jù)預先設定的參數(shù),運用SAS隨機函數(shù)ranexp(seed)先產(chǎn)生參數(shù)為M1、M2的指數(shù)分布數(shù)據(jù)x,受試者均勻入組,入組間隔為tn=T0/N0,入組順序即為隨機數(shù)的序號ob。令int=I-T0,令y=x-(N0-ob)×tn,將其中y>int者設定為y=int。

(2)重復步驟(1)產(chǎn)生參數(shù)為HM10、HM20的指數(shù)分布數(shù)據(jù)x,用其中y>int者填補參數(shù)為M1,M2數(shù)據(jù)中y=int的截尾數(shù)據(jù),即為期中分析的數(shù)據(jù)。每個搜索范圍內(nèi)有多個HM10,以HZ1(0.25～0.75)為例,搜索間隔為0.1年,則有0.25、0.35、0.45、0.55、0.65、0.75共6種填補參數(shù),根據(jù)HR求得HM20,并以HM10、HM20為參數(shù)產(chǎn)生6組混合指數(shù)分布隨機數(shù),用該6組隨機數(shù)中y>int的數(shù)據(jù)分別填補一次截尾數(shù)據(jù),產(chǎn)生6組期中分析數(shù)據(jù)。

模擬結(jié)果

表1　不同HR時4個期中分析點不同HZ的樣本量再估計結(jié)果

表2　不同HR時4個期中分析點HZ11的樣本量再估計結(jié)果

圖1　不同HR時再估計樣本量與期中分析時間點、HZ之間的關(guān)系

圖2　不同HR時再估計樣本量與期中分析時間點、HZ(包含參數(shù)真實值)之間的關(guān)系

圖3　各期中分析時間點上HZ11-HZ5再估計樣本量的α(a)和power(b)

模擬示例

討　　論

藥物臨床試驗的投入巨大,而樣本量是否充足直接關(guān)系到試驗的成敗。若在試驗設計階段對試驗參數(shù)的估計存在較大的不確定性,建議采用適應性設計,預設期中分析點對樣本量進行再估計,根據(jù)估計結(jié)果考慮是否增加樣本量。一般期中分析時只考慮是否增加樣本量以維持把握度,而不考慮降低樣本量[1]。Todd[6]等人提出的盲態(tài)下樣本量再估計方法,需要設定多個期中分析時間點進行多次期中分析,具體實施有一定的難度。本研究提出的盲態(tài)下生存資料樣本量再估計方法,只需要進行一次期中分析。本研究假設對照組與試驗組的樣本量之比為1:1,風險比為HR,α為0.05,power為0.8,并假設5個不同的M10,并以M10為均值暫定搜索范圍為M10±0.5M10,比較4個期中分析點上不同搜索范圍下樣本量再估計結(jié)果。模擬研究表明,在期中分析點I2～I4,當搜索范圍包含參數(shù)真實值時,經(jīng)過期中分析調(diào)整后的樣本量接近真實值,且樣本量估計結(jié)果變異較小。在I1時間點時經(jīng)填補的期中分析數(shù)據(jù)包含真實參數(shù)信息較少,搜索結(jié)果與填補參數(shù)相關(guān)性較大,因此樣本量再估計結(jié)果呈現(xiàn)出與搜索范圍值大小相關(guān)的現(xiàn)象,但在I2～ I4時間點時未出現(xiàn)該現(xiàn)象。

綜上所述,并考慮到實際應用中,若樣本量再估計結(jié)果大于初始樣本量,需要繼續(xù)追加樣本量,并完成相同的隨訪時間t,此時期中分析時間點越晚,則研究時間就越長,研究成本也隨之增加。因此,相似背景下建議在I2(完成1/4隨訪時間)時進行期中分析,根據(jù)參數(shù)估計結(jié)果重新確定搜索范圍,以獲得更加準確的樣本量。若是其他差異較大的背景時,1/4隨訪時間作為期中分析時間點劃分的參考依據(jù),可進行適當?shù)恼{(diào)整。若樣本量再估計結(jié)果比試驗設計之初的樣本量少,則維持原定樣本量不變;反之,則根據(jù)樣本量再估計結(jié)果考慮增加樣本量以維持把握度。本研究的不足之處是在HR恒定不變的假設下進行模擬研究,后續(xù)研究將會考慮HR改變時的情況。