區間刪失生存數據的統計分析方法及其應用*

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

區間刪失生存數據的統計分析方法及其應用*

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

在臨床研究中，當只知道事件發生在某一給定的時間區間內，而不知道其確切時間點時，將這類數據稱為區間刪失數據（interval censored data），表示為T∈（L，R］［1］，其中 T表示個體的生存時間，L表示刪失區間的下界，R表示上界。顯而易見，區間刪失包括左刪失和右刪失，臨床研究中，區間刪失現象比較常見，特別是在患者進行周期性隨訪的臨床試驗和隊列研究中。

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

原理與方法

1.非參數估計

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

取對數，得對數似然函數 log（L（p）），然后對 pj求偏導數有

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

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

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

（1）Turnbull算法

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

（2）EM修正算法

修正EM算法為Gentleman和Geyer在EM算法的基礎上提出的一種修正迭代算法［1］，初始值的極大似然估計的計算分為兩個部分——降維與最優化。計算出初始值后，使用EM算法對其進行迭代，計算出新的生存函數估計式與該區間內的死亡概率。當某些區間的死亡概率低于某個值的時候，先將區間的死亡概率歸為0，并驗證庫恩-塔克條件，即驗證此區間的死亡概率是否真正為0。若不滿足，則將區間內的概率密度函數加上一較小值，再進行EM算法迭代。當迭代之后最大的變化量小于某個規定值，并且滿足庫恩-塔克條件時，極大似然估計值收斂，即求出非參數估計結果。

2.區間刪失數據的非參數檢驗

針對 g組（g≥2）區間刪失數據［1］，設 hl（t）為第 l組在時間 t的風險函數，l＝1，2，3，…，g，其中 t（l）為第 l組的生存時間。

即有如下假設，

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

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

3.半參數估計的區間刪失ICM算法的Cox模型

基于所有個體的觀測時間t，風險函數h（t）、基線風險函數h0（t），當引入協變量Z及協變量系數β時，Cox比例風險函數的表達式為 h（t｜Z）＝h0（t）exp（βTZ），基線生存函數S0（t）與協變量Z及協變量系數β存在關系式 S（t｜Z）＝S0（t）exp｜（βTZ），此時，區間刪失生存函數的對數極大似然函數變為:

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

實例分析

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

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

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

1.生存率的非參數估計

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

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

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

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

表1 SUN模型與取中點的log-rank的組間比較結果

3.以性別、歷期、年齡作為分組因素做半參數回歸模型

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

表2 ICM算法與取中點Cox模型的參數估計及95%置信區間

討 論

區間刪失數據在醫學領域中是一種常見的數據類型，但是醫務工作者常將其簡化成右刪失的形式，再采用Kaplan-Meier估計、Log-rank方法、Cox比例風險模型進行統計分析，從本文實例分析以及 Dorey［3］與Rücker［4］等人所做相關模擬可知，區間刪失數據若簡化成右刪失的形式，會造成不合理的結果，因此對于區間刪失數據處理，采用專門的分析方法是十分必要的。

本文所介紹的區間刪失數據屬于II型區間刪失數據［12］。抽涉及的非參數估計、組間比較、半參數回歸模型是較為主流的方法，幾乎都能通過現行的SAS軟件、R軟件進行實現，例如本文使用R軟件的interval，intcox等程序包，能夠為行業內相關人士在具體的研究當中提供一定的幫助，得出更合理的估計、以及恰如其分的統計學結論。

［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型區間刪失數據的生存分析.中國衛生統計，2016，33（2）:357-361.

國家自然科學基金（81202288）；廣州市科技計劃（2012J5100023）；廣東省科技計劃（2010B031600100）

1.南方醫科大學公共衛生學院生物統計學系、廣東省熱帶病研究重點實驗室（510515）

2.暨南大學經濟學院統計學系

＃共同第一作者

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

（責任編輯:郭海強）