李文靜，鄧文麗，章婷婷

(江西師范大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，江西南昌330022)

0 引言

在臨床實驗或醫(yī)學(xué)研究中，由于客觀因素的限制，失效時間常常不能直接觀測到，而只能知道它在某一個區(qū)間內(nèi)，這類數(shù)據(jù)在統(tǒng)計學(xué)上被稱為區(qū)間刪失數(shù)據(jù)(interval-censored data).如在一些傳染性疾病感染時間的研究中，實驗對象被放入傳染源后，只能知道在某個觀察點實驗對象是否已染上疾病，染上疾病的具體時間卻無法觀測到，所以只能推測出從接觸傳染源到染上傳染病所經(jīng)歷的時間落在某個區(qū)間內(nèi).

區(qū)間刪失數(shù)據(jù)存在于許多應(yīng)用領(lǐng)域中，因此，這引發(fā)了一些統(tǒng)計學(xué)者對相關(guān)問題的研究.Huang Jian等［1］對區(qū)間刪失數(shù)據(jù)的分類及對應(yīng)的統(tǒng)計方法進行了較為詳細地描述.Sun Jian-guo［2］較為全面和系統(tǒng)地概括了區(qū)間刪失數(shù)據(jù)分析中涉及到的基本概念和方法.呂秋萍等［3］運用無偏轉(zhuǎn)換思想構(gòu)造了區(qū)間刪失數(shù)據(jù)函數(shù)的均值估計，并在此基礎(chǔ)上對所構(gòu)造的估計量方差進行了研究.在區(qū)間刪失數(shù)據(jù)的研究中，許多學(xué)者都是基于失效時間變量T和刪失時間變量C相互獨立的假定進行研究的，稱這種刪失情況為獨立刪失或非信息刪失(Independent Censoring，Noninformative Censoring).然而，在實際問題中，這個假定常常會遭到質(zhì)疑.如在對某種疾病的治療中，由于病情惡化或者是已接受的治療方案不奏效，從而導(dǎo)致病人退出治療，這種情況通常預(yù)示著該病人的存活時間會比較短，即刪失的個體對應(yīng)的生存時間更短.相反地，有些病人的退出可能因為病情好轉(zhuǎn)，不需要進一步治療，這種情況的刪失個體的生存時間可能會較長.和獨立刪失相反的是非獨立刪失(Dependent Censoring)，或稱為信息刪失(Informative Censoring).如果對信息刪失數(shù)據(jù)仍采用獨立刪失下的統(tǒng)計分析方法，則可能會得到有偏或者無效的結(jié)論.

在信息刪失數(shù)據(jù)的研究中，對失效時間和刪失時間相依性的假定是至關(guān)重要的.正確的假定可以提高估計的效率，得到更好的統(tǒng)計結(jié)論;不合適的假定可能會導(dǎo)致錯誤的結(jié)論.在實際應(yīng)用中，由于造成信息刪失的應(yīng)用背景和原因的不同，失效時間和刪失時間相依的形式和程度也變得非常復(fù)雜，很難準確估計.敏感性分析可以評價相依關(guān)系的假定對統(tǒng)計分析結(jié)果造成的影響.王純杰［4］基于Copula函數(shù)的一些性質(zhì)，給出了非參數(shù)模型下的信息區(qū)間刪失數(shù)據(jù)分布函數(shù)的相合估計.F.Siannis等［5］對失效時間和刪失時間的相依關(guān)系進行了假定，引入了標示相關(guān)程度的參數(shù)和偏度函數(shù)，且對參數(shù)估計受相依程序的影響進行了敏感性分析.Y.Park等［6］在單個總體和2個總體的情形下，分別對獨立刪失和信息右刪失混合數(shù)據(jù)下的相關(guān)估計問題進行了敏感性分析.Zhang Zhi-gang等［7］在正態(tài)脆弱模型假定下對I型信息區(qū)間刪失數(shù)據(jù)的比例風(fēng)險模型進行了敏感性分析.Huang Xue-lin等［8］基于連接函數(shù)(Copula)對信息右刪失數(shù)據(jù)下的比例風(fēng)險模型的估計問題進行了敏感性分析.

本文擬基于連接函數(shù)對信息區(qū)間刪失數(shù)據(jù)下失效時間的生存函數(shù)進行估計，并在3種不同連接函數(shù)的情形下關(guān)于相依關(guān)系對參數(shù)估計所造成的影響進行敏感性分析.

1 方法

記Ti為第i個個體的失效時間，Ci為第i個個體的刪失時間.試驗中得到的獨立同分布觀測值為{(ci，δi)，i=1，2，…，n}，其中

假設(shè)Ti和Ci的邊際分布函數(shù)分別為F(·)和G(·)，g(·)是Ci的邊際密度函數(shù)，這里i=1，2，…，n.

基于觀測樣本，可以構(gòu)造似然函數(shù):

當(dāng)Ti和Ci相互獨立時，基于觀測樣本的似然函數(shù)為

當(dāng)Ti和Ci不相互獨立時，給定1個有參數(shù)α的連接函數(shù)H(u，v，α)，假設(shè) Ti和 Ci的聯(lián)合分布函數(shù)為J(t，c)=P(Ti≤t，Ci≤c)=H(F(t)，G(c)，α)，聯(lián)合生存函數(shù)為

S(t，c)=P(Ti＞ t，Ci＞ c)=1-F(t)-

G(c)H(F(t)，G(c))，

則第i個個體被刪失的概率為

同理可得，第i個個體失效的概率為

由此可知，

綜上所述，當(dāng)失效時間和刪失時間不相互獨立時，樣本的觀測似然函數(shù)可以表示為

在信息區(qū)間刪失數(shù)據(jù)中，刪失時間能夠完全觀測到，所以G(·)可以直接用它的經(jīng)驗分布函數(shù)代替，其中

似然函數(shù)(1)可以轉(zhuǎn)化為

如果對失效時間的分布形式掌握的信息不多，則通常會考慮用非參數(shù)模型直接估計失效時間的分布函數(shù).類似于文獻［2］給出的獨立區(qū)間刪失數(shù)據(jù)非參數(shù)極大似然估計的方法，可以在(2)式中利用鄧文麗等［9］提出的一類保序最優(yōu)化問題的迭代算法得到分布函數(shù)F的估計.當(dāng)已知一些影響失效時間的協(xié)變量時，比例風(fēng)險模型和加速失效模型是廣泛接受的2類半?yún)?shù)模型.如果假定協(xié)變量的影響滿足比例風(fēng)險模型，則在似然函數(shù)的表達式中，邊際分布函數(shù)F(·)可以用含回歸系數(shù)和基準風(fēng)險函數(shù)的分布函數(shù)表達式代替，然后在似然函數(shù)(2)中，通過迭代的方法得到相關(guān)的估計;如果假定協(xié)變量的影響滿足加速失效模型，則在似然函數(shù)的表達式中，邊際分布函數(shù)F(·)可以用含回歸系數(shù)和隨機誤差項的分布函數(shù)表達式代替，然后在似然函數(shù)(2)中，通過迭代的方法得到相關(guān)的估計.張連增等［10］基于極大似然法研究了Copula的參數(shù)和半?yún)?shù)方法的估計效果.如果失效時間的分布函數(shù)形式已知，而只待估其中包含的參數(shù)，則利用似然函數(shù)(2)就可以得出參數(shù)的極大似然估計.

2 數(shù)值模擬

在實際應(yīng)用中由于失效時間和刪失時間相依的形式和程度非常復(fù)雜，很難準確估計，所以通過敏感性分析來評價相依關(guān)系的假定對統(tǒng)計分析結(jié)果造成的影響.

模擬計算中失效時間T采用威布爾分布隨機生成，因為它的危險率不是常數(shù)，所以，與指數(shù)分布相比，它有較廣闊的應(yīng)用，將其用于調(diào)查深槽輪滾珠軸承的疲勞壽命，或?qū)⑵溆糜诿鑼戨娮庸艿氖?威布爾的分布函數(shù)為 F(t)=1-e-(λt)γ，其中，γ 是分布曲線的形狀參數(shù)，λ是尺度參數(shù).模擬計算中選取了γ=2，λ=0.5.刪失時間C的邊際分布選取的是(0，A)上的均勻分布，調(diào)整A的大小可以改變刪失的比例.

失效時間和刪失時間的相關(guān)性選用阿基米德連接函數(shù)來描述.這里選取了 Clayton、Gumbe-Hougard和Frank 3種連接函數(shù).

D.G.Clayton［11］給出在 τ=1/(1+2α)下的Copula函數(shù):

E.J.Gumbel等［12］給出 τ=1-1/α 下的Copula函數(shù):

H(u，v，α)=exp{- ［(-log u)α+

(-log v)α］1/α}(α ≥ 1).

2.1.1 30 g/L甲基二磺隆對麥田禾本科雜草的防效經(jīng)過藥前雜草基數(shù)調(diào)查可知，試驗藥劑處理區(qū)、對照藥劑處理區(qū)及空白處理區(qū)野燕麥發(fā)生基數(shù)為：17.92，18.83，20.42，20.83，20.25，19.83 株；雀麥的發(fā)生基數(shù)為：18.92，19.08，16.92，20.00，20.08，20.08 株。

M.J.Frank［13］給出的 Copula 函數(shù):

H(u，v，α)=log{1+(αu-1)(αv-1)

(α -1)}(α ＞ 0，α ≠1)，其中τ=1+4γ-1［D1(γ)-1］，γ=-log α，D1(γ)=

R.B.Nelsen［14］對于連接函數(shù)的相關(guān)性質(zhì)和特殊的連接函數(shù)進行了詳細介紹.

下面主要是通過數(shù)值計算分析連接函數(shù)的選取對參數(shù)γ和λ的估計產(chǎn)生的影響.這里的穩(wěn)健性分析包括參數(shù)敏感性分析和連接函數(shù)敏感性分析.

2.1 參數(shù)的敏感性分析

分別取 τ=0.8、τ=0.5、τ=0.2的 Frank Copula作為連接函數(shù)，T服從γ=2，λ=0.5的威布爾分布，C服從(0，37)上的均勻分布，生成容量為200的樣本，刪失比例P(T＜C)為0.5.在本文方法中，選取 Frank連接函數(shù)，τ=0.8.模擬次數(shù)為1000，得到上述情況下λ～和γ～的均值、標準差和偏差的估計值(見表1).

表1 總體的參數(shù)τ變化下參數(shù)γ和λ的估計

由表1可以看出:如果生成樣本的Frank連接函數(shù)的參數(shù)τ為0.2、0.5、0.8，采用獨立刪失的估計方法，得到的γ和λ估計量的偏差都較大，特別是γ估計值的偏差很大.而采用本文方法(選取Frank連接函數(shù)，參數(shù)τ=0.8)得到的估計量都比較理想，其估計值的偏差遠遠小于獨立刪失下估計值的偏差.由此可見，在失效時間和刪失時間不相互獨立的情況下，采用獨立刪失方法進行估計會得到不理想的估計結(jié)果，因此，應(yīng)該采用帶相關(guān)性假定的模型進行分析.

2.2 改變連接函數(shù)下參數(shù)估計的敏感性分析

分別取Clayton Copula和Gumbel Copula作為連接函數(shù)，τ=0.8，T服從γ=2，λ=0.5的威布爾分布，C服從(0，37)上的均勻分布，生成容量為200的樣本，刪失比例P(T＜C)為0.5.在估計方法中采用Frank連接函數(shù)，τ=0.8.2種數(shù)據(jù)集下λ～和γ～的均值、標準差如表2所示.

表2 總體的連接函數(shù)形式變化時參數(shù)γ和λ的估計

由表2可以看出:當(dāng)τ=0.8時，如果生成樣本的連接函數(shù)分別選取Clayton和Gumbel連接函數(shù)，采用獨立刪失的估計方法，得到的γ和λ估計量的偏差都很大.而采用本文方法(選取Frank連接函數(shù)，參數(shù)τ=0.8)得到的估計量都比較理想，估計量的偏差比較小.由此可見，在失效時間和刪失時間不相互獨立的情況下，采用獨立刪失方法進行估計可能會得到不理想的估計，所以應(yīng)該采用帶相關(guān)性假定的模型進行分析.

其次，在連接函數(shù)的選取對估計量的影響方面，當(dāng)連接函數(shù)的假定和總體不一致時本文方法能夠得到較穩(wěn)健的估計量.

2.3 不同刪失比例下參數(shù)估計的敏感性分析

選取C服從均勻U(0，3.3)和U(0，4.0)，刪失比例P(T＜C)分別為0.3、0.7.T服從γ=2，λ=0.5的威布爾分布，連接函數(shù)為Frank，τ=0.8，生成容量為200的數(shù)據(jù)集，估計不同數(shù)據(jù)集下λ～和γ～的均值、標準差，估計結(jié)果如表3所示.

表3 不同刪失比例下參數(shù)γ和λ的估計

由表3可以看出:如果連接函數(shù)及其參數(shù)τ的假定都是正確的，則在不同的刪失比例下，本文方法都能夠得到較好的估計，但采用獨立刪失方法得到的估計卻不理想.

綜合上述模擬計算的結(jié)果，可以得出:1)本文方法的參數(shù)估計效果比獨立刪失方法的參數(shù)估計效果更好;2)由連接函數(shù)不能準確識別或者參數(shù)τ不能正確識別所導(dǎo)致的偏差小于由獨立刪失錯誤假設(shè)所引起的偏差;3)當(dāng)連接函數(shù)或參數(shù)τ的假定發(fā)生偏差時，本文方法依然能夠較穩(wěn)健.

3 討論

通過上面的敏感性分析，可以看出在帶信息的刪失數(shù)據(jù)下估計參數(shù)的標準差小于在獨立情況下估計參數(shù)的標準差.且對于不同連接函數(shù)、相關(guān)系數(shù)以及刪失比例，效果都較穩(wěn)健.另外，除了假設(shè)分布函數(shù)為威布爾分布外，還關(guān)于對數(shù)正態(tài)、指數(shù)分布等分布函數(shù)作了估計，效果也較好.因此，文本提供的方法具有一定的實用價值.

本文的工作還有許多地方可以進一步深入地研究，如在本文方法的框架下，繼續(xù)解決半?yún)?shù)模型的分布函數(shù)估計［15］;考慮有協(xié)變量影響下基于連接函數(shù)的帶信息刪失的非參數(shù)估計等.

［1］ Huang Jian，Wellner JA.Interval censored survival data:a review of recent progress［M］.New York:Springer-Verlag，1997:123-169.

［2］Sun Jianguo.The statistical analysis of interval-censored failure time data［M］.New York:Springer-Verlag，2006.

［3］呂秋萍，鄧文麗.區(qū)間刪失數(shù)據(jù)函數(shù)的均值估計［J］.江西師范大學(xué)學(xué)報:自然科學(xué)版，2011，35(1):96-100.

［4］王純杰.基于Copula函數(shù)的相依刪失數(shù)據(jù)的非參數(shù)統(tǒng)計推斷［D］.長春:吉林大學(xué)，2012.

［5］Siannis F，Copas J，Lu Baobing.Sensitivity analysis for informative censoring in parametric survivalmodels［J］.Biostatistics，2005，6(1):77-91.

［6］ Park Y，Lee Jenwei.One-and two-sample nonparametric inference procedures in the presence of amixture of independent and dependent censoring［J］.Biostatistics，2006，7(2):252-267.

［7］Zhang Zhigang，Sun Liuquan，Sun Jianguo，et al.Regression analysis of failure time data with informative interval censoring［J］.Statistics in Medicine，2007，26(12):2533-2546.

［8］ Huang Xuelin，Zhang Nan.Regression survival analysis with an assumed copula for dependent censoring:a sensitivity analysis approach ［J］.Biometrics，2008，64(4):1090-1099.

［9］鄧文麗，朱瑩瑩.一類保序最優(yōu)化問題的迭代算法［J］.統(tǒng)計與決策，2011(14):10-11.

［10］張連增，胡祥.Copula的參數(shù)與半?yún)?shù)估計方法的比較［J］.統(tǒng)計研究，2012，31(2):91-95.

［11］Clayton D G.Amodel for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence ［J］.Biometrika，1978，65(1):141-151.

［12］ Hougaard P.A class ofmultivariate failure time distributions［J］.Biometrika，1986，73(3):671-678.

［13］Frank M J.On the sumultaneous association of F(x，y)and x+y-F(X，Y)［J］.Aequationes Mathematicae，1979，21(41):37-38.

［14］Nelsen R B.An introduction to copulas［M］.2nd ed.New York:Springer-Verlag，2006.

［15］楊金英，趙培信.缺失數(shù)據(jù)下ρ～混合誤差線性模型的參數(shù)估計［J］.西南大學(xué)學(xué)報:自然科學(xué)版，2012，34(9):35-37.