部分區間刪失數據下廣義指數分布的參數估計及應用

董小剛, 彭小草, 蔣京京, 王純杰

(長春工業大學 數學與統計學院, 長春 130012)

0 引 言

基于參數模型對各類數據進行統計推斷是生存分析的重要內容之一[1-4]. 廣義指數分布[5]是一個重要的參數模型, 廣泛應用于物理、 壽命試驗和經濟學等領域. 廣義指數分布不僅是指數分布的推廣, 還具有許多與Gamma分布相似的性質, 以及類似于Weibull分布的分布函數, 能很好地規避Gamma分布與Weibull分布的不足, 關于其理論性質[6]和參數估計方法[7-9]的研究目前已有很多結果[10-12]. Gupta等[10]利用極大似然估計(MLE)、 矩估計和最小二乘估計等6種方法估計廣義指數分布的參數, 考慮了模型參數α和β分別固定時的情形, 并分析比較了6種方法的估計效果; Alizadeh等[12]將廣義指數模型應用到飛機空調系統故障次數試驗數據中, 通過Q-Q圖(Quantile Quantile plot)說明MLE方法可以提供更好的擬合效果.

在生存分析和壽命數據研究中, 數據類型多為分組、 刪失或截斷數據. 其中, 基于分組數據、 右刪失數據、 區間刪失數據和截斷數據等在廣義指數分布下的研究已取得了一些成果. 文獻[13-16]用Monte-Carlo模擬、 Newton-Raphson迭代算法、 EM(expectation maximum)算法求解了參數的極大似然估計; 文獻[17]用混合Gibbs算法、 MCMC(Markov chain Monte Carlo)方法求解了參數的Bayes估計； 文獻[18-19]對該類數據同時應用多種參數估計方法并比較了方法的優劣; 文獻[20]對該類數據應用多種模型并比較了模型的擬合效果. 針對部分區間刪失數據, 一些研究者利用Weibull參數模型[21]、 比例風險模型[22-23]、 加性風險模型[24]、 加速失效時間模型[25-26]、 治愈模型[27]等進行半參數模型擬合, 將該類數據與多種模型結合求解參數的極大似然估計和Bayes估計, 并用多種指標評價模型的適用性和有效性. 在該類研究中, 糖尿病數據集[22,28]和AIDS數據集[24,26]是常用的實例數據集. Saeed等[21]應用Weibull參數模型分別在不同刪失比下采用左點和右點填充技術處理數據, 計算參數的點估計和區間估計, 并給出了利用不同情形下的生存曲線圖估計模型的效果.

目前, 關于部分區間刪失數據的廣義指數參數模型研究尚未見文獻報道, 該數據類型包括完全數據、 右刪失和區間刪失數據, 數據類型更靈活. 基于此, 本文擬在部分區間刪失數據下, 在廣義指數參數模型中考慮尺度參數是否受協變量影響建立兩種模型, 并對參數進行極大似然估計. 模擬實驗和實例研究結果驗證了該方法的有效性.

1 數據、 模型及其參數估計

1.1 部分區間刪失數據

下面舉例說明部分區間刪失數據: 在1933—1972年丹麥糖尿病發病研究[22,28]中, 研究人員記錄的患者糖尿病發病時間發生在某兩次臨床檢測之間, 如果能準確觀察到患者患有糖尿病的時間, 則數據為確切觀測數據； 如果患者在第一次檢查前已發現患有糖尿病, 則該數據為左刪失數據； 如果患者糖尿病發病時間發生在兩次臨床檢測之間, 則該數據為區間刪失數據； 如果第二次檢查發現患者仍未患糖尿病, 則該數據為右刪失數據.

1.2 廣義指數分布

本文擬對部分區間刪失數據在廣義指數分布假設下進行參數估計和實例分析.令T=(t1,t2,…,tN)T表示感興趣的失效時間隨機變量, 假設感興趣的失效時間服從形狀參數和尺度參數分別為α和λ的廣義指數分布GE(α,λ), 則其相關函數如下: 對于第i個個體, 其分布函數為

F(ti;α,λ)=[1-exp{-λti}]α,ti>0,α>0,λ>0;

(1)

密度函數為

(2)

生存函數為

S(ti;α,λ)=1-F(ti)=1-[1-exp{-λti}]α,ti>0,α>0,λ>0;

(3)

風險函數為

(4)

1.3 參數估計

在全數據類型下, 廣義指數分布對應的似然函數為

1.3.1 廣義指數參數模型

(8)

其中

(9)

其中

(10)

(11)

將該方差-協方差矩陣求得平方根后取其對角線元素即可得對應的標準誤差.

1.3.2 廣義指數尺度參數回歸模型

(13)

(14)

(15)

將該矩陣求平方根后再取對角線元素可得對應的標準誤差.

Newton-Raphson迭代算法步驟如下：

1) 給出參數迭代的初始值θ0=(λ0,α0)T；

2) 根據對數似然函數計算未知參數的一階偏導數U(θ0)和二階偏導數G(θ0)；

2 模擬實驗

下面對本文提出的參數估計方法進行模擬實驗, 生成不同樣本量的隨機數.感興趣的失效時間T由廣義指數分布中產生, 為產生區間刪失數據, 假設相鄰檢查之間的時間間隔服從均勻分布U(0,a), 其中a為可調節參數, 通過調節參數a改變刪失比例, 每個個體的觀測時間數據為1加均值為10的Poisson隨機數.若觀測時間大于感興趣的失效時間, 則大于失效時間的最小觀測數據為左刪失數據； 若觀測時間小于感興趣的失效時間, 則小于失效時間的最大觀測數據為右刪失數據； 否則為區間刪失數據, 即包含失效時間的最小區間為區間刪失數據.

為考察低刪失比率和高刪失比率下估計方法的有效性, 本文將刪失比分別設置為0.2,0.5,0.8, 并考慮樣本量分別為200,400,800,1 000情形下的模擬效果, 模擬循環500次, 分別得到不同刪失比及不同樣本量下的模擬結果, 其中BIAS表示估計參數的平均偏差, SE表示估計參數的標準差, SEE表示估計標準誤差的均值, CP表示95%置信區間的覆蓋率.

模擬Ⅰ： 當尺度參數不受協變量影響時,T服從廣義指數分布, 參數真值α=1,λ=1, 模擬結果列于表1.

表1 模擬Ⅰ的參數估計結果

由表1和表2可見, 參數估計值均較接近真實值, 樣本量越大, 偏差越小, 因此符合大樣本性質, 由SE和SEE較接近可再一次驗證理論結果的正確性, 且其覆蓋率接近95%置信區間. 隨著確切觀測數據的增加, 估計的效果越來越好, 表明對部分區間刪失數據建立廣義指數模型具有穩定性和有效性, 同時也驗證了Newton-Raphson算法在極大似然估計中能準確有效地得到廣義指數分布下基于部分區間刪失數據的參數估計值.

表2 模擬Ⅱ的參數估計結果

3 實證分析

3.1 丹麥糖尿病數據集

糖尿病是一種以正常胰島素分泌紊亂為特征的慢性疾病, 主要分為兩種類型： Ⅰ型, 即胰島素依賴型糖尿病, 是最嚴重的類型, 主要發生在生命前期； Ⅱ型, 即非胰島素依賴型糖尿病, 是一種病情較輕的類型, 主要發生在生命后期. 在一次檢測過程中, 如果患者連續4次尿樣分析中(時間至少間隔1個月)每24 h含蛋白質超過0.5 g, 則確定為糖尿病腎病. 在丹麥糖尿病實例(Diabetes)數據[22,28]中, 所有病人在進入研究或者在研究結束時都已患有糖尿病腎病, 表明數據中無右刪失數據, 有595個確切觀測值, 136個區間刪失觀測值, 這731名病人中, 有男性454名, 女性277名, 222名患者的年齡小于10歲, 509名患者年齡在10～30歲之間.

下面考慮部分區間刪失數據, 當尺度參數不受協變量影響時, 應用廣義指數模型及其參數估計方法, 對上述731名患者的觀測數據建立模型, 對性別分組分析, 估計結果列于表3. 由表3可見, 男性和女性生存函數的形狀參數和尺度參數較接近.

表3 當尺度參數不受協變量影響時, Diabetes數據估計結果

根據性別分組得到男性和女性的Fisher信息矩陣分別為

(16)

方差-協方差矩陣分別為

(17)

表4 當尺度參數受協變量影響時, Diabetes數據估計結果

該模型下的Fisher信息陣為

(18)

方差-協方差矩陣為

(19)

由于丹麥糖尿病實例數據集中的性別協變量因素對生存時間的影響不顯著, 所以本文又考慮了AIDS數據集.

3.2 AIDS數據集

表5 當尺度參數受協變量影響時, AIDS數據估計結果

該數據下的Fisher信息矩陣為

(20)

方差-協方差矩陣為

(21)

綜上所述, 本文在部分區間刪失數據下根據尺度參數是否受協變量的影響建立了廣義指數模型, 研究了模型下的極大似然估計, 并利用Newton-Raphson算法給出了參數的估計值, 利用Fisher信息矩陣求得了參數的方差-協方差矩陣, 對不同刪失比及不同樣本量下的模型進行了模擬研究, 最后將該參數估計方法應用到丹麥糖尿病實例數據集和AIDS數據集中. 模擬實驗和實例分析結果驗證了模型的有效性.