區(qū)間刪失生存數(shù)據(jù)的統(tǒng)計(jì)分析方法及其應(yīng)用*

李麗賢 湯 茗 曾彥彥 沈羅英 李釗洪 陳慧林 郭勝楠 陳金寶 侯雅文 陳 征△

區(qū)間刪失生存數(shù)據(jù)的統(tǒng)計(jì)分析方法及其應(yīng)用*

李麗賢1＃湯 茗1＃曾彥彥1沈羅英1李釗洪1陳慧林1郭勝楠1陳金寶1侯雅文2陳 征1△

在臨床研究中，當(dāng)只知道事件發(fā)生在某一給定的時(shí)間區(qū)間內(nèi)，而不知道其確切時(shí)間點(diǎn)時(shí)，將這類數(shù)據(jù)稱為區(qū)間刪失數(shù)據(jù)（interval censored data），表示為T∈（L，R］［1］，其中 T表示個(gè)體的生存時(shí)間，L表示刪失區(qū)間的下界，R表示上界。顯而易見，區(qū)間刪失包括左刪失和右刪失，臨床研究中，區(qū)間刪失現(xiàn)象比較常見，特別是在患者進(jìn)行周期性隨訪的臨床試驗(yàn)和隊(duì)列研究中。

在處理區(qū)間刪失數(shù)據(jù)時(shí)，很多研究者往往為簡(jiǎn)便起見，直接用刪失區(qū)間的中點(diǎn)或右端點(diǎn)作為生存時(shí)間的估計(jì)值，再利用類似于處理“右刪失資料”的方法估計(jì)生存率，并據(jù)此進(jìn)行參數(shù)估計(jì)和假設(shè)檢驗(yàn)［2］。Dorey［3］等通過(guò)模擬研究發(fā)現(xiàn)，若將刪失區(qū)間的中點(diǎn)作為觀察時(shí)間，用右刪失方法處理，則會(huì)高估生存率。此外，Rücker［4］等發(fā)現(xiàn)，若將刪失區(qū)間的右端點(diǎn)作為觀察時(shí)間，進(jìn)行Kaplan-Meier估計(jì)，則會(huì)過(guò)低估計(jì)誤差方差，從而易得出假陽(yáng)性的結(jié)果。可見，這些簡(jiǎn)單的處理方式不僅會(huì)降低生存率估計(jì)的準(zhǔn)確性，而且會(huì)影響估計(jì)的精度，是不合理的。因此，采用專門的統(tǒng)計(jì)分析方法來(lái)處理區(qū)間刪失數(shù)據(jù)是非常必要的。本文將介紹有關(guān)區(qū)間刪失數(shù)據(jù)的非參數(shù)方法以及半?yún)?shù)模型，并以“靜脈注射毒品成癮患者的HIV感染情況的生存分析”為例，展示其應(yīng)用。

原理與方法

1.非參數(shù)估計(jì)

設(shè)第 i個(gè)個(gè)體的生存時(shí)間為 Ti（i＝1，…，n），（Li，Ri］為 Ti所屬的刪失區(qū)間，S（t）＝P（T＞t）為生存函數(shù)（即生存率）。令｛τj（即｛0＝τ0＜τ1＜τ2＜…＜τj＜…＜τm｝）等于｛0，Li，Ri；i＝1，…，n｝為唯一順序元素，可見｛τj為研究中能觀察到的全部時(shí)間點(diǎn)（也就是所有觀察區(qū)間的端點(diǎn)），τj為第j＋1個(gè)時(shí)間點(diǎn)。使 αij＝I（（τj－1，τj］∈（Li，Ri］），i＝1，…，n，j＝1，…，m，I（）為指示變量，對(duì)于第 i個(gè)個(gè)體，若區(qū)間（τj－1，τj］包含在區(qū)間（Li，Ri］中，則 αij＝1，否則為 0，由 αij可得知，一個(gè)發(fā)生在區(qū)間（Li，Ri］上的事件，是否發(fā)生在（τj－1，τj］上。定義 pj＝S（τj－1）－S（τj），j＝1，…，m，p＝（p1，…，pj，…，pm）T，pj在這里表示第 j個(gè)時(shí)間點(diǎn)的死亡概率。似然函數(shù)可表示為

取對(duì)數(shù)，得對(duì)數(shù)似然函數(shù) log（L（p）），然后對(duì) pj求偏導(dǎo)數(shù)有

這里，ηi表示第 i個(gè)個(gè)體的死亡概率，dj則為在（τj－1，τj］區(qū)間上的所有個(gè)體的ηi倒數(shù)之和。

由公式（3）和（4），可知當(dāng) pj＞0時(shí)，μj＝0，dj＝μ0＝n，而當(dāng) pj＝0時(shí)，μj≥0，dj＝μ0－μj≤n。當(dāng) μj≥0和 dj＋μj－μ0＝0稱之為滿足庫(kù)恩-塔克條件。因此，對(duì)所有的 j，當(dāng) dj＝n（j＝1，…，m）時(shí)，p為所求的 NPMLE。值得注意的是，Peto［7］提出:只有當(dāng) τj－1＝Li和 τj＝Rk（i≠k）時(shí)［7］，有 pj≥0，但不排除滿足以上兩個(gè)條件后，仍有一些pj＝0。

求S（t）和p的過(guò)程中，可通過(guò)不同的迭代計(jì)算方法進(jìn)行求解。其中，Turnbull［8］提出的修正乘積極限估計(jì)算法（Turnbull算法），以及 Gentleman［9］在 EM算法的基礎(chǔ)上提出的修正EM算法應(yīng)用較為廣泛。

（1）Turnbull算法

Turnbull［8］提出了一種類似乘積極限估計(jì)的方法，該生存函數(shù)估計(jì)式可以由自一致性算法來(lái)估計(jì)。迭代的初始值由乘積極限法求得，設(shè)為的第r次迭代結(jié)果，其迭代過(guò)程可表示為:

（2）EM修正算法

修正EM算法為Gentleman和Geyer在EM算法的基礎(chǔ)上提出的一種修正迭代算法［1］，初始值的極大似然估計(jì)的計(jì)算分為兩個(gè)部分——降維與最優(yōu)化。計(jì)算出初始值后，使用EM算法對(duì)其進(jìn)行迭代，計(jì)算出新的生存函數(shù)估計(jì)式與該區(qū)間內(nèi)的死亡概率。當(dāng)某些區(qū)間的死亡概率低于某個(gè)值的時(shí)候，先將區(qū)間的死亡概率歸為0，并驗(yàn)證庫(kù)恩-塔克條件，即驗(yàn)證此區(qū)間的死亡概率是否真正為0。若不滿足，則將區(qū)間內(nèi)的概率密度函數(shù)加上一較小值，再進(jìn)行EM算法迭代。當(dāng)?shù)笞畲蟮淖兓啃∮谀硞€(gè)規(guī)定值，并且滿足庫(kù)恩-塔克條件時(shí)，極大似然估計(jì)值收斂，即求出非參數(shù)估計(jì)結(jié)果。

2.區(qū)間刪失數(shù)據(jù)的非參數(shù)檢驗(yàn)

針對(duì) g組（g≥2）區(qū)間刪失數(shù)據(jù)［1］，設(shè) hl（t）為第 l組在時(shí)間 t的風(fēng)險(xiǎn)函數(shù)，l＝1，2，3，…，g，其中 t（l）為第 l組的生存時(shí)間。

即有如下假設(shè)，

H1∶各組死亡率 hl（t）不等或不全相等，l＝1，2，3，…，g

利用非參數(shù)極大似然估計(jì)函數(shù)公式（1）L（p），可求出得分統(tǒng)計(jì)量 U＝（U1，U2，…，Ug），第 l組的得分統(tǒng)計(jì)量Ul可表示為:為權(quán)重，cjl表示在原假設(shè)下第 l個(gè)組在區(qū)間（τj－1，τj］內(nèi)的死亡人數(shù)的期望值。cj表示第j個(gè)區(qū)間內(nèi)的總死亡人數(shù)的期望值，類似地ajl和aj表示危險(xiǎn)集的期望值，構(gòu)造檢驗(yàn)統(tǒng)計(jì)量χ2＝U∑U′，服從自由度為g-1的卡方分布。其中，∑為U的協(xié)方差矩陣，可通過(guò)置換方法（permutation method）得到［10］。針對(duì) U統(tǒng)計(jì)量的不同權(quán)重 wj，主要有 Sun模型與 Wilcoxon-type模型［2］，其權(quán)重分別為wj＝1和wj＝（tj－1）。

3.半?yún)?shù)估計(jì)的區(qū)間刪失ICM算法的Cox模型

基于所有個(gè)體的觀測(cè)時(shí)間t，風(fēng)險(xiǎn)函數(shù)h（t）、基線風(fēng)險(xiǎn)函數(shù)h0（t），當(dāng)引入?yún)f(xié)變量Z及協(xié)變量系數(shù)β時(shí)，Cox比例風(fēng)險(xiǎn)函數(shù)的表達(dá)式為 h（t｜Z）＝h0（t）exp（βTZ），基線生存函數(shù)S0（t）與協(xié)變量Z及協(xié)變量系數(shù)β存在關(guān)系式 S（t｜Z）＝S0（t）exp｜（βTZ），此時(shí)，區(qū)間刪失生存函數(shù)的對(duì)數(shù)極大似然函數(shù)變?yōu)?

針對(duì)對(duì)數(shù)極大似然函數(shù)L（S0，β），可利用ICM（iterative convex minorant）算法進(jìn)行迭代計(jì)算［11］。參數(shù)的初始值轉(zhuǎn)化為右刪失數(shù)據(jù)通過(guò)經(jīng)典比例風(fēng)險(xiǎn)模型計(jì)算得出，并將Breslow估計(jì)值作為基準(zhǔn)風(fēng)險(xiǎn)h0（t）的初始值。在對(duì)β進(jìn)行顯著性檢驗(yàn)，使用輪廓似然法和信息矩陣等方法都能估算出的方差。但由于生存數(shù)據(jù)計(jì)算量大，以上方法迭代較慢，bootstrap方法成為一個(gè)較好的選擇。Efron在1986年的隨機(jī)模擬中發(fā)現(xiàn)［2］，在不同的刪失率下，bootstrap方法計(jì)算速度快，且能保持較小的偏倚。因此綜合考慮采用bootstrap方法進(jìn)行方差估計(jì)。

實(shí)例分析

對(duì)某戒毒中心的881名靜脈注射毒品患者經(jīng)過(guò)戒毒治療后的HIV感染數(shù)據(jù)進(jìn)行分析［6］，此研究的起始時(shí)間為戒毒治療的開始，事件為發(fā)生HIV感染。HIV病毒直接的感染情況可通過(guò)定期的血清檢測(cè)來(lái)確定，以血清檢測(cè)由陰轉(zhuǎn)陽(yáng)的時(shí)間來(lái)確定HIV感染。在該研究中，研究中心以月為單位定期檢測(cè)血清情況。因事件發(fā)生（HIV感染）的確切時(shí)間并不能被直接觀察到，即可將這些患者感染HIV看作區(qū)間刪失型數(shù)據(jù)，其刪失區(qū)間的左端點(diǎn)為最近一次血清檢測(cè)為陰性的時(shí)間，右端點(diǎn)為第一次血清檢測(cè)為陽(yáng)性的時(shí)間。若在觀察期內(nèi)未檢測(cè)到血清抗體呈陽(yáng)性，則該患者為右刪失數(shù)據(jù)。如第694號(hào)個(gè)體，第一次檢測(cè)為陽(yáng)性的時(shí)間為第42月，在此之前最近一次血清檢測(cè)陰性的時(shí)間為第29月，故 t694∈（29，42］。

在觀察隨訪的881個(gè)樣本中，707例為男性，174為女性；患者年齡的中位數(shù)為19歲，其中400例患者小于等于19歲，481例患者大于19歲。576例患者在觀察期內(nèi)感染HIV（事件發(fā)生），305例未感染（右刪失）。另外，按研究的觀察時(shí)間，我們將歷期（calendar period）分為1972至1985年和1986至1997年兩大組，分別為539例和342例。

實(shí)例分析中將對(duì)患者的生存時(shí)間進(jìn)行非參數(shù)估計(jì)，并將性別、歷期、年齡（分為小于等于19歲與大于19歲）作為3個(gè)因素，進(jìn)行組間比較和多因素回歸模型分析。

1.生存率的非參數(shù)估計(jì)

對(duì)于HIV抗體由陰轉(zhuǎn)陽(yáng)情況，利用中點(diǎn)算法、Turnbull算法、修正EM算法，對(duì)患者的生存時(shí)間進(jìn)行非參數(shù)估計(jì)，其生存率的比較見圖1A。Turnbull算法與修正EM算法計(jì)算出的生存曲線除在尾部有些許差異，其余幾乎重合。此外，類似于 Dorey［3］等人，將刪失區(qū)間的中點(diǎn)當(dāng)作右刪失時(shí)間點(diǎn)，用Kaplan-Meier估計(jì)法進(jìn)行計(jì)算，得到生存曲線，結(jié)果顯示，中點(diǎn)處理后的生存率大體較高于Turnbull方法與修正EM算法，特別是在前期。不同的分組變量的生存率（圖1B、1C、1D）估計(jì)亦顯示中點(diǎn)算法得到的生存率大體上要高于修正EM算法得到的生存率。

圖1 三種方法所得生存率估計(jì)的總比較及分亞組比較

2.以性別、歷期、年齡為分組因素分別做組間比較

采用權(quán)重為1的Sun模型，以歷期為分組因素，進(jìn)行組間比較，可得出檢驗(yàn)統(tǒng)計(jì)量Z＝3.456，P＝0.001，顯示不同的歷期間生存率有統(tǒng)計(jì)學(xué)差異。若對(duì)生存時(shí)間取中點(diǎn)化為右刪失數(shù)據(jù)后，使用log-rank檢驗(yàn)比較組間差異，則得χ2＝0.083，P＝0.773提示差異沒有統(tǒng)計(jì)學(xué)意義。但是，結(jié)合圖1C可看出兩個(gè)歷期的生存率是不同的，從側(cè)面上反映了1986年前進(jìn)入戒毒中心的患者比1985年后進(jìn)入的患者更容易感染HIV病毒。按性別分組，進(jìn)行組間比較，兩種算法均顯示有統(tǒng)計(jì)學(xué)差異，從圖1B中也可以看出，男性的生存率明顯高于女性。由圖1D可知，年齡分組基本沒有差異，特別是中點(diǎn)算法，兩條生存曲線基本重合，兩種算法的組間比較亦顯示沒有統(tǒng)計(jì)學(xué)差異，P值均大于0.05。

表1 SUN模型與取中點(diǎn)的log-rank的組間比較結(jié)果

3.以性別、歷期、年齡作為分組因素做半?yún)?shù)回歸模型

通過(guò)ICM算法，對(duì)患者的性別、歷期、年齡分組進(jìn)行半?yún)?shù)估計(jì)，計(jì)算參數(shù)β，并用bootstrap方法進(jìn)行40000次有放回抽樣，得出β的95%置信區(qū)間。此外，將區(qū)間中點(diǎn)轉(zhuǎn)換為右刪失數(shù)據(jù)，并計(jì)算其Cox模型參數(shù)。如表2所示，兩種算法均顯示不同性別的HIV感染情況有統(tǒng)計(jì)學(xué)差異，女性比男性更容易感染HIV病毒，圖1B亦可看出男性的生存曲線明顯分離高于女性的；年齡之間顯示生存率沒有統(tǒng)計(jì)學(xué)差異，兩種算法的95%置信區(qū)間均包含了0；對(duì)于歷期，雖然取中點(diǎn)的Cox比例風(fēng)險(xiǎn)模型中歷期的置信區(qū)間包含0，提示不同的歷期之間的HIV感染情況沒有統(tǒng)計(jì)學(xué)差異，但是從圖1C中可以看出，不同歷期之間的生存率差異還是比較大的，同時(shí)，ICM算法的Cox模型顯示不同歷期的生存率是有統(tǒng)計(jì)學(xué)差異的，所以在這里取中點(diǎn)的Cox比例風(fēng)險(xiǎn)模型應(yīng)用于區(qū)間刪失數(shù)據(jù)的分析并不是很合理。

表2 ICM算法與取中點(diǎn)Cox模型的參數(shù)估計(jì)及95%置信區(qū)間

討 論

區(qū)間刪失數(shù)據(jù)在醫(yī)學(xué)領(lǐng)域中是一種常見的數(shù)據(jù)類型，但是醫(yī)務(wù)工作者常將其簡(jiǎn)化成右刪失的形式，再采用Kaplan-Meier估計(jì)、Log-rank方法、Cox比例風(fēng)險(xiǎn)模型進(jìn)行統(tǒng)計(jì)分析，從本文實(shí)例分析以及 Dorey［3］與Rücker［4］等人所做相關(guān)模擬可知，區(qū)間刪失數(shù)據(jù)若簡(jiǎn)化成右刪失的形式，會(huì)造成不合理的結(jié)果，因此對(duì)于區(qū)間刪失數(shù)據(jù)處理，采用專門的分析方法是十分必要的。

本文所介紹的區(qū)間刪失數(shù)據(jù)屬于II型區(qū)間刪失數(shù)據(jù)［12］。抽涉及的非參數(shù)估計(jì)、組間比較、半?yún)?shù)回歸模型是較為主流的方法，幾乎都能通過(guò)現(xiàn)行的SAS軟件、R軟件進(jìn)行實(shí)現(xiàn)，例如本文使用R軟件的interval，intcox等程序包，能夠?yàn)樾袠I(yè)內(nèi)相關(guān)人士在具體的研究當(dāng)中提供一定的幫助，得出更合理的估計(jì)、以及恰如其分的統(tǒng)計(jì)學(xué)結(jié)論。

［1］Fay MP，Shaw PA.Exact and asymptotic weighted log-rank tests for interval censored data:the interval R package.Journal of Statistical Software，2010，36（2）:i2.

［2］Mongoué-TchokotéS，Kim J.New statistical software for the proportional hazards model with current status data.Computational Statistics＆Data Analysis，2008，52（9）:4272-4286.

［3］Dorey FJ，Little RJA，Schenker N.Multiple imputation for threshold crossing data with interval censoring.Statistics in Medicine，1993，12（17）:1589-1603.

［4］Rücker G，Messerer D.Rem ission duration:an example of interval censored observations.Statistics in Medicine，1988，7（11）:1139-1145.

［5］Gómez G，Luz Calle M，Egea JM，et al.Risk of HIV infection as a function of the duration of intravenous drug use:a non-parametric Bayesian approach.Statistics in Medicine，2000，19（19）:2641-2656.

［6］Hanson MA.Invexity and the Kuhn-Tucker Theorem.Journal of Mathematical Analysis and Applications，1999，236（2）:594-604.

［7］Peto R.Experimental survival curves for interval-censored data.Journal of the Royal Statistical Society.Series C（Applied Statistics），1973，22（1）:86-91.

［8］Turnbull BW.Nonparametric estimation of a survivorship function with doubly censored data.Journalof the American Statistical Association，1974，69（345）:169-173.

［9］Gentleman R，Geyer CJ.Maximum likelihood for interval censored data:consistency and computation.Biometrika，1994，81（3）:618-623.

［10］Heinze G，Gnant M，Schemper M.Exact log-rank tests for unequal follow-up.Biometrics，2003，59（4）:1151-1157.

［11］Pan W.Extending the iterative convex m inorant algorithm to the Cox model for interval-censored data.Journal of Computational＆Graphical Statistics，1999，8（1）:109-120.

［12］梁潔，王彤，崔燕.II型區(qū)間刪失數(shù)據(jù)的生存分析.中國(guó)衛(wèi)生統(tǒng)計(jì)，2016，33（2）:357-361.

國(guó)家自然科學(xué)基金（81202288）；廣州市科技計(jì)劃（2012J5100023）；廣東省科技計(jì)劃（2010B031600100）

1.南方醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計(jì)學(xué)系、廣東省熱帶病研究重點(diǎn)實(shí)驗(yàn)室（510515）

2.暨南大學(xué)經(jīng)濟(jì)學(xué)院統(tǒng)計(jì)學(xué)系

＃共同第一作者

△通信作者:陳征，E-mail:zchen@smu.edu.cn

（責(zé)任編輯:郭海強(qiáng)）