競爭風(fēng)險數(shù)據(jù)中累積發(fā)生率置信區(qū)間的估計研究*

陳金寶 侯雅文 陳 征△

在疾病的預(yù)后研究中，觀察到的終點事件往往不止一種，并且各終點事件呈相互競爭狀態(tài)，稱為競爭風(fēng)險 (competing risks)[1-4]。例如在SARS流行期間，患者從入院接受治療到研究結(jié)束，可能出現(xiàn)治愈和死亡這兩個競爭的終點事件[5]。累積發(fā)生率 (cumulative incidence function,CIF) 是描述性競爭風(fēng)險分析中最重要的指標之一[6]，CIF置信區(qū)間是可以按預(yù)先給定的概率 (95%、99%等) 確定包含CIF的一個范圍，一般可表示為 (L,U)，分別表示區(qū)間下限 (L≥ 0) 和區(qū)間上限 (U≤ 1)。經(jīng)典CIF區(qū)間是基于假定CIF近似服從正態(tài)分布構(gòu)造得到的，然而卻可能出現(xiàn)區(qū)間下限L<0或區(qū)間上限U> 1越界的異常情況，特別是在小樣本時[7]，而邏輯轉(zhuǎn)換[7]、反正弦平方根轉(zhuǎn)換[10]可以避免區(qū)間出現(xiàn)越界情況。對此本文將基于對CIF進行不同轉(zhuǎn)換，分別為線性轉(zhuǎn)換 (經(jīng)典)、對數(shù)轉(zhuǎn)換、雙對數(shù)轉(zhuǎn)換、反正弦平方根轉(zhuǎn)換以及邏輯轉(zhuǎn)換[8-9]，構(gòu)造5種不同形式的置信區(qū)間，均先構(gòu)造CIF轉(zhuǎn)換形式置信區(qū)間，再轉(zhuǎn)換為其原尺度下的置信區(qū)間。然后基于不同區(qū)間可信度 (90%、95%和99%)、事件發(fā)生率、樣本量、刪失率以及多個時刻點t上，通過模擬研究，綜合評價和比較上述5種置信區(qū)間的錯誤覆蓋率，最后進行一個實例分析。

方法介紹

本文研究只有一個組別，假設(shè)樣本量為n，可能發(fā)生的終點事件有K種，第j個個體的觀測數(shù)據(jù)為(tj,δj),j=1,2,…,n。Tj和Cj分別表示第j個個體的生存時間和刪失時間，并且Tj和Cj相互獨立，則tj=min(Tj,Cj)，δj=I(Tj≤Cj)，I(·)為指示變量，如果觀測到個體發(fā)生終點事件(包括興趣事件和競爭事件)，則定義為對應(yīng)第k種終點事件，k=1,2,…,K，否則為右刪失。競爭風(fēng)險中CIF是指在給定時刻點t之前，發(fā)生第k種終點事件的失效概率，則第k種終點事件CIF估計式為：

其中dj表示在時間tj上合并所有終點事件類型的事件數(shù)。設(shè)定置信區(qū)間可信度為1-α，則時刻點t上CIF的線性轉(zhuǎn)換 (經(jīng)典) 100(1-α)%的雙側(cè)置信區(qū)間 (CIline) 為：

(1)

由置信區(qū)間的一般形式，可以構(gòu)造如下分別基于對數(shù)轉(zhuǎn)換、雙對數(shù)轉(zhuǎn)換、反正弦平方根轉(zhuǎn)換以及邏輯轉(zhuǎn)換的CIF置信區(qū)間，這四種置信區(qū)間均是非對稱區(qū)間。

(2)

(3)

(4)

(5)

模擬研究

通過Monte-Carlo模擬比較上述5種置信區(qū)間的性能。模擬研究只有一組，設(shè)定樣本量為50，100和200，生存時間分布服從參數(shù)為1的指數(shù)分布，刪失時間分布服從均勻分布U(0,a)，改變a值獲得不同刪失率 (0%、10%、20%、30%、40%、50%和60%)。假定競爭風(fēng)險模擬數(shù)據(jù)[12]只有一個興趣事件 (類型1) 和一個競爭事件 (類型2)，則興趣事件和競爭事件理論CIF分別為I1(t)=p1(1-e-t)和I2(t)=(1-p1)(1-e-t)，其中p1和1-p1分別表示興趣事件和競爭事件最大的事件發(fā)生概率，p1取值為0.4、0.6和0.8。通過預(yù)模擬實驗發(fā)現(xiàn)，模擬數(shù)據(jù)時間四分位數(shù)近似為0.25、0.60和0.95，所以模擬比較的三個時間點t為0.25、0.60和0.95。雙側(cè)置信區(qū)間范圍為90%、95%和99%，對應(yīng)的檢驗水準α為0.10、0.05和0.01。通過判斷興趣事件CIF的置信區(qū)間是否包含興趣事件理論CIF，定義評價指標為錯誤覆蓋率，模擬循環(huán)次數(shù)為10000次。

為了能更加清晰綜合評價5種轉(zhuǎn)換置信區(qū)間的表現(xiàn)，使用方差分析技術(shù) (ANOVA)[13]綜合評價錯誤覆蓋率。首先定義錯誤覆蓋率減去檢驗水準作為結(jié)果變量M，即最終評價指標為平均偏差，不同置信區(qū)間范圍有不同檢驗水準。方法表現(xiàn)越好則其期望E(M) 越接近于常數(shù)0。模型中考慮了4個影響因素：Test、P1、Num、Cen和T，分別代表5種置信區(qū)間形式、3種興趣事件最大發(fā)生概率、3種樣本量、7種刪失率和3種時刻點t，擬合不同的不帶截距項模型。模型1表示在控制了其他影響因素下，研究不同興趣事件最大發(fā)生概率對不同置信區(qū)間的影響，模型2～4與1類似，而模型5是控制所有影響因素，研究不同置信區(qū)間的邊際效應(yīng)。

模型1:E(M)=Test×P1+Num+Cen+T;

模型2:E(M)=Test×Num+P1+Cen+T;

模型3:E(M)=Test×Cen+Num+P1+T;

模型4:E(M)=Test×T+Cen+Num+P1;

模型5:E(M)=Test+Cen+Num+P1+T。

表1是不同區(qū)間范圍下，5種CIF置信區(qū)間形式分別擬合模型1～5的平均偏差。不同區(qū)間范圍下，CIline相比于其他置信區(qū)間都有著最大正數(shù)平均偏差，即有著最大的錯誤覆蓋率；CIarcsin次之，且輕微高估錯誤覆蓋率；CIlogit有最小負數(shù)平均偏差，顯得過于保守；CIlog顯得不穩(wěn)健，在90%和95%范圍輕微保守，而99%范圍則輕微高估；CIloglog輕微保守，最接近于常數(shù)0，顯得最為穩(wěn)健可靠。不同模型下，從模型1看出隨著興趣事件最大發(fā)生概率的增大，5種方法基本都呈現(xiàn)越接近0的趨勢；模型2看出隨著樣本量的增加，除了CIarcsin在90%時遠離0，其余情況5種區(qū)間基本都呈現(xiàn)越接近于0的趨勢；模型3看出5種區(qū)間的變化趨勢波動，但依舊可以發(fā)現(xiàn)CIloglog顯得最為穩(wěn)健可靠，特別是在高刪失時；模型4看出5種區(qū)間變化趨勢不一，但CIloglog總是和大多數(shù)的區(qū)間變化趨勢保持一致；模型5結(jié)果基本和不同區(qū)間范圍的結(jié)果一致。

實例分析

一項關(guān)于霍奇金病患者的治療研究[14]，隨訪時間從1968年到1986年。采用外周血液移植治療方式共有49人，興趣事件是移植后發(fā)生慢性移植物抗宿主反應(yīng) (n=45)，移植后死亡或復(fù)發(fā)是競爭事件 (n=3)，其余為右刪失 (n=1)。對外周血液移植治療第0.5、0.75和1年求其CIF值的5種雙側(cè)95%置信區(qū)間。從各時刻點區(qū)間看出 (圖1)，時刻點越大區(qū)間寬度越小，除CIline外其余置信區(qū)間均非對稱。在第0.5和0.75年5種置信區(qū)間上下限均未出現(xiàn)越界情形，然而在第1年時，CIlog上限L=1.007和CIline上限L=1.003均出現(xiàn)了大于1的異常情況，CIarcsin、CIlogit和CIloglog上限L分別為0.982、0.972和0.971，均保持在[0，1]范圍內(nèi)，很好地克服出現(xiàn)越界情況，結(jié)合模擬研究，推薦以CIloglog為準。

表1 ANOVA技術(shù)綜合評價3種不同范圍的5種CIF置信區(qū)間的模擬結(jié)果的平均偏差

P1：擬合模型1；Num：擬合模型2；Cen：擬合模型3；T：擬合模型4；最后一行：擬合模型5。

圖1 外周血液組CIF及其95%CI范圍

討 論

在疾病的預(yù)后研究中，競爭風(fēng)險型數(shù)據(jù)是常見的生存壽命數(shù)據(jù)類型之一，并且其CIF是最重要的描述性指標，因此準確地給出其100(1-α)%雙側(cè)置信區(qū)間也是重要的，然而CIline(經(jīng)典) 上限可能出現(xiàn)越界的異常情形，對此考慮了對CIF進行對數(shù)轉(zhuǎn)換、雙對數(shù)轉(zhuǎn)換、反正弦平方根轉(zhuǎn)換以及邏輯轉(zhuǎn)換來構(gòu)造置信區(qū)間，后三種可以有效地避免越界情形的出現(xiàn)。根據(jù)模擬結(jié)果，CIline和CIarcsin均有較大的正數(shù)偏差，CIlog則易出現(xiàn)波動，CIlogit有最小負數(shù)偏差，只有CIloglog偏差最接近于期望常數(shù)0，顯得最為穩(wěn)健可靠。實例分析中發(fā)現(xiàn)CIline和CIlog無法克服出現(xiàn)越界情況，且CIlog表現(xiàn)不穩(wěn)定，CIarcsin、CIlogit和CIloglog則很好地克服出現(xiàn)越界情況。值得注意，當高刪失率或小樣本時，CIloglog均表現(xiàn)相對穩(wěn)健較佳，另外小樣本下CIarcsin表現(xiàn)也相對穩(wěn)健較佳。在醫(yī)學(xué)研究中，特別是當數(shù)據(jù)是小樣本或高刪失時，要謹慎使用CIline(經(jīng)典)。另外，也可以利用置信區(qū)間上下限值，對兩組間時刻點上CIF進行差異性比較的初步分析。

[1] Lau B,Cole SR,Gange SJ.Competing risk regression models for epidemiologic data.Am J Epidemiol,2009,170(2):244-256．

[2] Austin PC,Lee DS,Fine JP.Introduction to the Analysis of Survival Data in the Presence of Competing Risks.Circulation,2016,133(6):601-609．

[3] 楊召,王少明,粱赫,等.競爭風(fēng)險型數(shù)據(jù)統(tǒng)計分析理論研究進展.中國衛(wèi)生統(tǒng)計,2016,33(6):1088-1091.

[4] 盧梓航,周立志,韓棟,等.競爭風(fēng)險型數(shù)據(jù)的統(tǒng)計推斷處理及應(yīng)用.現(xiàn)代預(yù)防醫(yī)學(xué),2013,40(5):804-807.

[5] 陳征,Nakamura T.基于競爭風(fēng)險理論和概要型數(shù)據(jù)的病死率估計模型.中國衛(wèi)生統(tǒng)計,2010,27(3):249-252.

[6] Kalbfleisch JD,Prentice RL.The Statistical Analysis ofFailure Time Data.New York:Wiley,2002.

[7] Hong Y,Meeher WQ.Confidence interval procedures for system reliabilityand applications to competing risks models.Lifetime Data Anal,2014,20(2):161-184.

[8] 陳金寶,邱李斌,王北琪,等.固定點處組間生存率比較的統(tǒng)計檢驗法.中華流行病學(xué)雜志,2015,36(2):186-188.

[9] 項永兵,高玉堂,金凡,等.生存率置信區(qū)間的五種估計方法.中華流行病學(xué)雜,1995,16(5):306-309.

[10] Choudhury JB.Non-parametric confidence interval estimation for competing risks analysis:application to contraceptive data.Stat Med,2002,21(8):1129-1144.

[11] Aalen O.Nonparametric estimation of partial transition probabilities in multiple decrement models.Ann Stat,1978,6(3):534-545.

[12] Beyersmann J,Latouche A,Buchholz A,et al.Simulating competing risks data in survival analysis.Stat Med,2009,28(6):956-971.

[13] Klein JP,Logan B,Harhoff M,et a1.Analyzing survival curves at a fixed point in time.Stat Med,2007,26(24):4505-4519．

[14] Pintilie M.Competing Risks:A Practical Perspective.England:John Wiley & Sons,2006.