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

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

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

方法介紹

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

(1)

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

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

(2)

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

(3)

S(t)為生存函數(shù)，μ為區(qū)間[0,τ]上生存函數(shù)的積分，即對(duì)應(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為通過(guò)修正EM算法分別估計(jì)兩組的生存曲線圖，圖1B、圖1C陰影部分面積分別為兩組限制性平均生存時(shí)間的估計(jì)值，其中限制性平均生存時(shí)間及其方差的估計(jì)值如下：

(7)

(8)

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

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

(9)

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

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

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

(10)

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

模擬研究

本文采用Monte-Carlo模擬來(lái)探索區(qū)間刪失生存數(shù)據(jù)中RMST檢驗(yàn)法的檢驗(yàn)效能和I類錯(cuò)誤[5]，評(píng)價(jià)其穩(wěn)健性和適用性。生存時(shí)間T：由參數(shù)為(α+βzi)的指數(shù)分布產(chǎn)生，zi是組別指示變量(zi=0表示對(duì)照組，1為試驗(yàn)組)，即對(duì)照組、試驗(yàn)組的生存時(shí)間T分別服從風(fēng)險(xiǎn)率為α、α+β的指數(shù)分布，且β表示試驗(yàn)組與對(duì)照組風(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]；模擬通過(guò)θ1,θ2的不同取值控制左刪失(生存時(shí)間T小于刪失區(qū)間下限)、區(qū)間刪失(生存時(shí)間T位于刪失區(qū)間之間)及右刪失(生存時(shí)間T大于刪失區(qū)間上限)的比例。本研究設(shè)置α=0.2，n1=n2=30,100,200，循環(huán)次數(shù)設(shè)為1000次。

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

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

*：α=0.2。

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

實(shí)例分析

對(duì)某戒毒中心的940名靜脈注射毒品患者經(jīng)戒毒治療后的HIV感染數(shù)據(jù)進(jìn)行分析[14]，該研究起始時(shí)間是戒毒治療的開(kāi)始，終點(diǎn)事件定義為發(fā)生HIV感染，研究中心以月為單位定期檢查血清情況，由此確定在每個(gè)觀察區(qū)間內(nèi)是否發(fā)生終點(diǎn)事件。因發(fā)生HIV感染的確切時(shí)間并不能被直接觀察到，即可將這些患者感染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可見(jiàn)：當(dāng)τ取兩組最大區(qū)間端點(diǎn)的較小值(即τ=166)時(shí)，在男性中限制性平均生存時(shí)間μ0=68.106個(gè)月(圖1B中生存曲線下面積值)，在女性中限制性平均生存時(shí)間是μ1=55.534個(gè)月(圖1C中生存曲線下面積值)，故男性與女性的限制性平均生存時(shí)間差值為12.572個(gè)月，且其對(duì)應(yīng)的95%可信區(qū)間為(2.822,22.322)，RMST檢驗(yàn)統(tǒng)計(jì)量z=2.527，P=0.011，提示在戒毒治療開(kāi)始至第166個(gè)月，男性組、女性組的平均HIV感染時(shí)間有統(tǒng)計(jì)學(xué)差異，且男性較女性長(zhǎng)12.572個(gè)月；假如研究者關(guān)心的時(shí)間點(diǎn)是截止到第50個(gè)月(即τ=50)時(shí)，男性組、女性組的平均感染時(shí)間分別為36.928、 29.859個(gè)月，其差值及95%可信區(qū)間為7.069(3.598,10.540)，z=3.992，P<0.001；同理，假如研究者關(guān)心的時(shí)間點(diǎn)是截止到第100個(gè)月(即τ=100)時(shí)，男性組和女性組的平均感染時(shí)間分別為54.721、45.596個(gè)月，其差值及95%可信區(qū)間為9.125(2.365,15.885)，z=2.646，P=0.008。對(duì)于Sun模型，χ2=3.448，P=0.063，該檢驗(yàn)法只能給出兩組差異的整體性檢驗(yàn)結(jié)果，即提示兩組HIV感染時(shí)間沒(méi)有統(tǒng)計(jì)學(xué)差異。

討 論

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