黃興輝 呂晶晶 楊紫荊 侯雅文 陳 征△

【提 要】 目的 針對區(qū)間刪失生存數(shù)據(jù)的分析研究，提出限制平均生存時間(restricted mean survival time，RMST)的估計(jì)和兩組比較檢驗(yàn)。方法 利用修正EM算法進(jìn)行迭代并得到生存率估計(jì)值，并基于此估計(jì)值構(gòu)建RMST檢驗(yàn)統(tǒng)計(jì)量求得P值，通過Monte-Carlo模擬驗(yàn)證其統(tǒng)計(jì)性能。結(jié)果 本文提出的區(qū)間刪失生存數(shù)據(jù)中RMST檢驗(yàn)法的I類錯誤在0.05附近波動，且其檢驗(yàn)效能與現(xiàn)有常用方法Sun模型相當(dāng)。結(jié)論 針對區(qū)間刪失生存數(shù)據(jù)的組間比較，本文提出RMST檢驗(yàn)法不僅準(zhǔn)確估計(jì)各組生存率，且通過計(jì)算RMST能夠直觀解釋組間差異大小并由此作出統(tǒng)計(jì)推斷，具有較好的統(tǒng)計(jì)性能，為臨床研究者和病人提供決策依據(jù)。

在醫(yī)學(xué)臨床研究中，通常需要對事件發(fā)生時間進(jìn)行觀測，并進(jìn)一步分析評價(jià)。當(dāng)只知道事件發(fā)生在某一特定區(qū)間內(nèi)，而不知道其確切的時間點(diǎn)時，將這類數(shù)據(jù)稱為區(qū)間刪失數(shù)據(jù)[1-3]，表示為T∈(L,R)，其中T表示個體的生存時間，L表示刪失區(qū)間的下界，R表示上界。目前對于區(qū)間刪失數(shù)據(jù)，現(xiàn)有一類廣義log-rank檢驗(yàn)較為常用，例如Sun模型[4-5]，該方法基于單個區(qū)間內(nèi)實(shí)際死亡數(shù)與理論死亡數(shù)的差值構(gòu)造統(tǒng)計(jì)量，但是該類方法難以直觀解釋差異大小且可能會出現(xiàn)不收斂的問題；有研究[6-8]表明，限制平均生存時間(restricted mean survival time，RMST)表示從研究起始時間到一個特定時間點(diǎn)τ上的平均生存時間，即對應(yīng)于時間點(diǎn)τ前的生存曲線下面積大小，并且針對該指標(biāo)差值或者比值進(jìn)行檢驗(yàn)稱為RMST檢驗(yàn)，能夠簡單描述生存狀態(tài)下平均生存時間的差異大小，在進(jìn)行組間比較時應(yīng)用廣泛[9-10]，特別是當(dāng)風(fēng)險(xiǎn)率成比例假定不成立或者感興趣的事件數(shù)較少時，該指標(biāo)能夠更加有效反映預(yù)后效果[11-12]。本文利用修正EM算法進(jìn)行區(qū)間刪失數(shù)據(jù)的生存率估計(jì)，由此構(gòu)造RMST檢驗(yàn)統(tǒng)計(jì)量，對應(yīng)提出區(qū)間刪失生存數(shù)據(jù)下的RMST檢驗(yàn)法，通過模擬探索其穩(wěn)健性和適用性并進(jìn)行實(shí)例分析。

方法介紹

1.區(qū)間刪失生存數(shù)據(jù)中的RMST估計(jì)

(1)

對似然函數(shù)取對數(shù)，并對pj求偏導(dǎo)數(shù)得

i=1,…,n;l=j=1,…,m

(2)

如圖1所示，對某戒毒中心的940名靜脈注射毒品患者經(jīng)過戒毒治療后的HIV感染數(shù)據(jù)進(jìn)行分析[14]，假設(shè)截止時間點(diǎn)為τ(最后一個觀察者的生存時間或者是最后一個事件者的死亡時間或者是某個給定的時間點(diǎn)，本研究取兩組中最大區(qū)間端點(diǎn)的較小值)，則限制性生存時間為X=min(T,τ)，由此可得X的期望，即限制性平均生存時間可以表示如下：

(3)

S(t)為生存函數(shù)，μ為區(qū)間[0,τ]上生存函數(shù)的積分，即對應(yīng)于生存函數(shù)S(t)曲線下面積。

進(jìn)一步計(jì)算E(X2)，即

E(X2)=E(T2|T≤τ)Pr(T≤τ)+τ2Pr(T>τ)

(4)

由于Pr(T≤τ)=1-S(τ)，則

(5)

故X的方差可以估計(jì)如下：

Var(X)=RSDST2=E(X2)-[E(X)]2

(6)

其中，RSDST為限制性標(biāo)準(zhǔn)差。

圖1 經(jīng)戒毒治療后的HIV感染

圖1A為通過修正EM算法分別估計(jì)兩組的生存曲線圖，圖1B、圖1C陰影部分面積分別為兩組限制性平均生存時間的估計(jì)值，其中限制性平均生存時間及其方差的估計(jì)值如下：

(7)

(8)

2.區(qū)間刪失生存數(shù)據(jù)中的RMST檢驗(yàn)法

比較兩組間生存率的差異，原假設(shè)是在任意時刻點(diǎn)t上，兩組對應(yīng)的生存率相等，即H0:S1(t)=S2(t)，備擇假設(shè)H1:S1(t)≠S2(t)。根據(jù)限制性平均生存時間構(gòu)造統(tǒng)計(jì)量，即在時間點(diǎn)τ上計(jì)算兩組限制性平均生存時間的差值Δ：

(9)

S0(t)，S1(t)分別為第一組、第二組的生存函數(shù)；μ0，μ1分別為第一組、第二組的限制性平均生存時間。

根據(jù)修正EM算法可得：

根據(jù)delta法進(jìn)一步估計(jì)其方差：

(10)

在原假設(shè)、大樣本下該檢驗(yàn)統(tǒng)計(jì)量服從t分布或正態(tài)分布[10]。

模擬研究

本文采用Monte-Carlo模擬來探索區(qū)間刪失生存數(shù)據(jù)中RMST檢驗(yàn)法的檢驗(yàn)效能和I類錯誤[5]，評價(jià)其穩(wěn)健性和適用性。生存時間T：由參數(shù)為(α+βzi)的指數(shù)分布產(chǎn)生，zi是組別指示變量(zi=0表示對照組，1為試驗(yàn)組)，即對照組、試驗(yàn)組的生存時間T分別服從風(fēng)險(xiǎn)率為α、α+β的指數(shù)分布，且β表示試驗(yàn)組與對照組風(fēng)險(xiǎn)率的差值。刪失區(qū)間產(chǎn)生：首先產(chǎn)生分別服從均勻分布(1,θ1)和(1,θ2)的相互獨(dú)立U1,U2，其中θ1,θ2是大于1的常數(shù)；然后定義U為U1四舍五入后的值，定義V為max(U1+U2,U+1)并取四舍五入后的值，由此產(chǎn)生刪失區(qū)間(U,V]；模擬通過θ1,θ2的不同取值控制左刪失(生存時間T小于刪失區(qū)間下限)、區(qū)間刪失(生存時間T位于刪失區(qū)間之間)及右刪失(生存時間T大于刪失區(qū)間上限)的比例。本研究設(shè)置α=0.2，n1=n2=30,100,200，循環(huán)次數(shù)設(shè)為1000次。

表1展示了不同樣本量、刪失百分比和β取值時兩種檢驗(yàn)法的I類錯誤和檢驗(yàn)效能，當(dāng)α=0.2，β=-0.1，-0.05，0，0.05，0.1時，即固定第一組生存時間T服從風(fēng)險(xiǎn)率為0.2的指數(shù)分布，第二組生存時間T分別服從風(fēng)險(xiǎn)率分別為0.1，0.15，0.2，0.25，0.3的指數(shù)分布。表1第一列分別對應(yīng)模擬數(shù)據(jù)中左刪失、區(qū)間刪失、右刪失所占比例，本文考慮了(1/3,1/3,1/3)和(1/4,1/2,1/4)兩種刪失率組合。從表1可見，當(dāng)β=0時為I類錯誤：在不同刪失比例、不同樣本量組合下，兩種方法的I類錯誤都在0.05附近波動，RMST檢驗(yàn)法在刪失比例為(1/3,1/3,1/3)且樣本量為(30,30)時I類錯誤偏小。當(dāng)β≠0時為檢驗(yàn)效能：樣本量為(30,30)時，兩種方法的檢驗(yàn)效能交替較高；樣本量為(100,100)和(200,200)時，當(dāng)β=-0.1和0.1，Sun模型較高，當(dāng)β=-0.05和0.05，RMST檢驗(yàn)法較高。

表1 I類錯誤和檢驗(yàn)效能

*：α=0.2。

總體來說，本文提出的區(qū)間刪失生存數(shù)據(jù)中RMST檢驗(yàn)法與常用的Sun模型均能較好控制I類錯誤，并且除了小樣本外具有較高的檢驗(yàn)效能，具有較好的統(tǒng)計(jì)性能。

實(shí)例分析

對某戒毒中心的940名靜脈注射毒品患者經(jīng)戒毒治療后的HIV感染數(shù)據(jù)進(jìn)行分析[14]，該研究起始時間是戒毒治療的開始，終點(diǎn)事件定義為發(fā)生HIV感染，研究中心以月為單位定期檢查血清情況，由此確定在每個觀察區(qū)間內(nèi)是否發(fā)生終點(diǎn)事件。因發(fā)生HIV感染的確切時間并不能被直接觀察到，即可將這些患者感染HIV看作區(qū)間刪失型數(shù)據(jù)；若患者在觀察期內(nèi)未發(fā)生終點(diǎn)事件，則為右刪失數(shù)據(jù)。根據(jù)性別進(jìn)行分組，男性組有759例患者，女性組有181例患者。

表2 RMST檢驗(yàn)法的實(shí)例分析結(jié)果

*：Sunχ2=3.448,P=0.063。

由圖1和表2可見：當(dāng)τ取兩組最大區(qū)間端點(diǎn)的較小值(即τ=166)時，在男性中限制性平均生存時間μ0=68.106個月(圖1B中生存曲線下面積值)，在女性中限制性平均生存時間是μ1=55.534個月(圖1C中生存曲線下面積值)，故男性與女性的限制性平均生存時間差值為12.572個月，且其對應(yīng)的95%可信區(qū)間為(2.822,22.322)，RMST檢驗(yàn)統(tǒng)計(jì)量z=2.527，P=0.011，提示在戒毒治療開始至第166個月，男性組、女性組的平均HIV感染時間有統(tǒng)計(jì)學(xué)差異，且男性較女性長12.572個月；假如研究者關(guān)心的時間點(diǎn)是截止到第50個月(即τ=50)時，男性組、女性組的平均感染時間分別為36.928、 29.859個月，其差值及95%可信區(qū)間為7.069(3.598,10.540)，z=3.992，P<0.001；同理，假如研究者關(guān)心的時間點(diǎn)是截止到第100個月(即τ=100)時，男性組和女性組的平均感染時間分別為54.721、45.596個月，其差值及95%可信區(qū)間為9.125(2.365,15.885)，z=2.646，P=0.008。對于Sun模型，χ2=3.448，P=0.063，該檢驗(yàn)法只能給出兩組差異的整體性檢驗(yàn)結(jié)果，即提示兩組HIV感染時間沒有統(tǒng)計(jì)學(xué)差異。

討 論

本研究結(jié)合修正EM算法，在區(qū)間刪失生存數(shù)據(jù)中進(jìn)行RMST檢驗(yàn)法的發(fā)展與應(yīng)用。本文提出的RMST檢驗(yàn)法既保證了生存率的準(zhǔn)確估計(jì)，又對區(qū)間刪失生存數(shù)據(jù)中組間差異進(jìn)行直觀闡述和比較；不管風(fēng)險(xiǎn)率是否成比例，該檢驗(yàn)法均通過兩組生存曲線下面積差值構(gòu)建統(tǒng)計(jì)量，能夠直觀反映在τ時刻前平均生存時間的差異大小，并做假設(shè)檢驗(yàn)，在事件數(shù)較少時依然適用，為臨床研究者和病人提供平衡風(fēng)險(xiǎn)-成本-效益(risk-cost-benefit)并作出決策的依據(jù)。由上述實(shí)例分析結(jié)果可知，τ的選取對結(jié)果的影響較大。截止時間點(diǎn)τ需要在實(shí)驗(yàn)設(shè)計(jì)階段就明確定義，并且其不同取值對應(yīng)不同的限制性平均生存時間，從而導(dǎo)致檢驗(yàn)結(jié)果有所不同。本研究實(shí)現(xiàn)了在區(qū)間刪失生存數(shù)據(jù)中RMST檢驗(yàn)的拓展，并通過模擬驗(yàn)證其能夠較好地控制I類錯誤同時具有檢驗(yàn)效能，但尚未與其他現(xiàn)存檢驗(yàn)方法進(jìn)行全面比較。此外，本研究只考慮了兩組的情況，涉及多組的組間比較有待進(jìn)一步推導(dǎo)和完善。