左截斷右刪失數(shù)據(jù)下尺度參數(shù)與協(xié)變量相關時廣義指數(shù)分布的估計及應用

趙 波, 王純杰, 李 群, 佟知真

(長春工業(yè)大學 數(shù)學與統(tǒng)計學院, 長春 130012)

基于參數(shù)模型的統(tǒng)計推斷是生存分析領域中的一個重要分支[1-3]. 廣義指數(shù)分布[4]能在多方面彌補經(jīng)典參數(shù)分布的不足, 在航空、工程、醫(yī)學和生物科學等領域應用廣泛, 因此雙參數(shù)廣義指數(shù)分布在參數(shù)模型的研究中得到廣泛關注[5-9]. 關于不同數(shù)據(jù)類型的生存時間或失效時間特征的分析處理, 也是生存分析領域中的研究熱點. 其中, 基于右刪失數(shù)據(jù)、區(qū)間刪失數(shù)據(jù)、截斷數(shù)據(jù)以及雙刪失數(shù)據(jù)等廣義指數(shù)分布參數(shù)模型的研究已有許多結果[1,10-13]. 但基于左截斷右刪失數(shù)據(jù)的參數(shù)模型研究文獻報道較少[14-18].

本文主要基于左截斷右刪失數(shù)據(jù)和雙參數(shù)廣義指數(shù)模型, 針對尺度參數(shù)是否受到協(xié)變量因素影響建立兩種模型, 采用極大似然估計法估計參數(shù), 并利用Newton-Raphson迭代算法求解, 用隨機模擬分析驗證估計結果的有效性和一致性, 并將兩種模型分別應用于電子變壓器壽命數(shù)據(jù)集[18-19]和Channing house數(shù)據(jù)集[10]中, 給出其生存函數(shù)和風險函數(shù).

1 左截斷右刪失數(shù)據(jù)下廣義指數(shù)的參數(shù)估計

1.1 左截斷右刪失數(shù)據(jù)

圖1 左截斷右刪失數(shù)據(jù)實例Fig.1 Example of left truncated and right censored data

假設實驗開始時間點為a, 實驗結束時間點為b, 則如圖1所示的左截斷右刪失數(shù)據(jù)包含以下4種情形： 1) 完整數(shù)據(jù), 即興趣事件從開始到結束都發(fā)生在區(qū)間[a,b]內(nèi)； 2) 發(fā)生了右刪失, 即到實驗結束時興趣事件仍在繼續(xù)； 3) 發(fā)生了左截斷, 即在實驗開始時興趣事件已經(jīng)發(fā)生； 4) 完全左截斷, 即在實驗開始時興趣事件已經(jīng)結束. 實驗開始時, 情形4)已經(jīng)發(fā)生, 收集不到該情形下的數(shù)據(jù), 所以本文主要考慮情形1)～3).

1.2 廣義指數(shù)分布

廣義指數(shù)分布是由單參數(shù)指數(shù)分布向雙參數(shù)分布的一個推廣, 一般記作GE(α,λ), 與單參數(shù)指數(shù)分布同為連續(xù)型分布, 其分布函數(shù)為

F(t;λ,α)=(1-exp{-λt})α,t>0,α>0,λ>0,

概率密度函數(shù)為

生存函數(shù)為

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

風險函數(shù)為

其中α和λ分別表示雙參數(shù)廣義指數(shù)分布的形狀參數(shù)和尺度參數(shù). 當α=1時, 風險函數(shù)可表示為h(t;λ,1)=λ, 與時間變量無關, 密度函數(shù)表示為f(t;λ,1)=λexp{-λt}, 即廣義指數(shù)分布退化為指數(shù)分布.

1.3 參數(shù)估計

假設T表示觀測時間變量,ti=min{yi,Ci}(i=1,2,…,n)表示第i個個體從興趣事件發(fā)生到興趣事件結束的時長,yi表示觀測的興趣事件從開始到結束的時長(右刪失數(shù)據(jù)下的完全數(shù)據(jù));Ci表示右刪失時間. 當存在截斷時,τi表示左截斷時間. 設δi表示刪失指示函數(shù), 當δi=0時, 表示觀測數(shù)據(jù)發(fā)生右刪失, 當δi=1時, 表示觀測數(shù)據(jù)未發(fā)生右刪失;νi表示截斷指示函數(shù), 當νi=0時, 表示觀測數(shù)據(jù)發(fā)生左截斷, 當νi=1時, 表示觀測數(shù)據(jù)未發(fā)生左截斷. 當不考慮協(xié)變量存在時, 數(shù)據(jù)結構為D={ti,δi,νi,τi;i=1,2,…,n}； 當考慮協(xié)變量矩陣X存在時, 觀測數(shù)據(jù)結構為D={ti,δi,νi,τi,Xi;i=1,2,…,n}. 取集合S,S1和S2分別表示全集、左截斷數(shù)據(jù)集及S1關于S的補集(即未發(fā)生截斷, 右刪失數(shù)據(jù)和完整數(shù)據(jù)組成的集合), 記作S1={i:νi=0},S2={i:νi=1}.

在左截斷右刪失數(shù)據(jù)下, 參數(shù)的似然函數(shù)為

(1)

其中:θ表示未知參數(shù);G(τi)表示截斷部分的分布函數(shù). 整理可得

(2)

根據(jù)左截斷右刪失數(shù)據(jù)下常見分布形式截斷部分分布函數(shù)的取法[16-18,20], 取截斷點的分布與興趣時間一致為廣義指數(shù)分布, 即G(τi)=(1-exp{-λτi})α, 所以在左截斷右刪失數(shù)據(jù)下, 廣義指數(shù)分布的似然函數(shù)可表示為

因此, 對數(shù)似然函數(shù)為

1.3.1 廣義指數(shù)參數(shù)模型 由于不考慮分布的尺度參數(shù)受協(xié)變量因素的影響, 因此可直接根據(jù)該數(shù)據(jù)類型下廣義指數(shù)分布的對數(shù)似然函數(shù), 用極大似然估計法和Newton-Raphson迭代算法極大化求解, 估計參數(shù)θ=(λ,α). 對似然函數(shù)LL(λ,α;D)分別關于參數(shù)λ和α求偏導, 并令其偏導數(shù)為0：

(5)

因此關于參數(shù)θ=(α,β0,β1,…,βm)的一階導數(shù)為

(7)

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

1) 給定參數(shù)迭代的初始值θ0=(λ0,α0);

2 數(shù)值模擬

根據(jù)反解法生成廣義指數(shù)分布的隨機數(shù), 即F(t;θ)=1-S(t;θ)=U,U服從均勻分布U(0,1), 反解出分布函數(shù)中的事件發(fā)生時間t. 所以當尺度參數(shù)不受協(xié)變量因素影響時,t=-λ-1log(1-U1/α)； 當尺度參數(shù)受協(xié)變量因素影響時,

文獻[14]中生成左截斷右刪失隨機數(shù)的方法, 用R軟件實現(xiàn)隨機模擬. 考察不同樣本量n=100,300,500, 不同刪失比p=0.4,0.6和截斷比q=0.3下的參數(shù)估計效果.

當尺度參數(shù)不受協(xié)變量因素影響時, 假設興趣時間變量T服從參數(shù)為θ=(α,λ)的廣義指數(shù)分布, 截斷部分時間的分布與興趣時間變量T的分布一致, 但其截斷數(shù)受截斷比q控制. 刪失時間點分布服從參數(shù)為k的指數(shù)分布, 參數(shù)k的取值由刪失比p的取值決定.

模擬ⅠT服從參數(shù)為α=1.2,λ=0.6的廣義指數(shù)分布. 實驗模擬在不同樣本量n=100,300,500下, 根據(jù)刪失截斷分成兩組, 一組刪失比和截斷比分別為0.4,0.3, 另一組刪失比和截斷比分別為0.6,0.3, 循環(huán)1 000次, 求出參數(shù)估計值, 并與真實參數(shù)值進行比較分析, 模擬結果列于表1.

表1 模擬Ⅰ(α=1.2, λ=0.6)的模擬結果*

*EST表示樣本估計值; BIAS表示估計的樣本偏差; SE表示標準誤差; SEE表示標準差; RMS表示標準均方誤差; CP表示置信水平為95%置信區(qū)間的覆蓋概率.

當尺度參數(shù)受協(xié)變量因素影響時, 設興趣時間變量T服從參數(shù)為θ=(α,λi=exp{β0+β1x1i})的廣義指數(shù)分布, 協(xié)變量x1服從概率為0.5的二項分布, 即x1～B(n,1,0.5).

模擬Ⅱ 形狀參數(shù)α=1.2, 決定尺度參數(shù)的兩個回歸參數(shù)為β0=0.8,β1=0.6. 實驗模擬在樣本量n=100,300,500、刪失比和截斷比分別為0.4,0.3的條件下循環(huán)模擬1 000次, 模擬結果列于表2.

表2 模擬Ⅱ(α=1.2, β0=0.8, β1=0.6)的模擬結果*

*EST表示樣本估計值; BIAS表示估計的樣本偏差; SE表示標準誤差; SEE表示標準差; RMS表示標準均方誤差; CP表示置信水平為95%置信區(qū)間的覆蓋概率.

由表1和表2可見, 估計值基本接近真實值, SE,SEE,RMS的取值較接近, 估計穩(wěn)定, 且隨樣本量的增大, 估計的偏差和標準差均減小, 因此在該種數(shù)據(jù)類型下, 估計的結果有很好的穩(wěn)定性和一致性, 表明Newton-Raphson迭代算法在極大似然估計中能精確、有效地估計模型參數(shù).

3 實例分析

例1電子變壓器壽命數(shù)據(jù)集[18-19]. 該數(shù)據(jù)集共有286個觀測值, 正常失效數(shù)據(jù)39個, 發(fā)生刪失的數(shù)據(jù)247個, 發(fā)生截斷的數(shù)據(jù)167個, 所以其刪失比和截斷比分別為86.36%,58.39%. 用t表示跟蹤的壽命時長,τ表示觀測的截斷部分時間,δ,v分別表示刪失和截斷的示性變量, 發(fā)生刪失或者截斷, 則δ或v取0, 反之取1.

考慮在左截斷右刪失數(shù)據(jù)下, 廣義指數(shù)分布的參數(shù)模型, 可得變壓器的生存(可靠性)函數(shù)為

S(t;λ,α)=1-(1-exp{-λt})α,

危險率函數(shù)為

表3 變壓器壽命數(shù)據(jù)分析結果

例2美國Channing house 數(shù)據(jù)集[1,10]. 該數(shù)據(jù)集, 即1964-01至1975-07美國加州帕洛阿爾托市退休中心462名老年居民的死亡時間數(shù)據(jù)集. 其中數(shù)據(jù)只有實際壽命超過實驗起始點時才能觀測、記錄, 因此退休居民進入實驗的年齡即為所需截斷點. 進入年齡小于840個月(70 a)的個體組成的數(shù)據(jù)子集表示發(fā)生截斷, 截斷比為18.0%, 截斷部分用τ表示, 示性函數(shù)v取0表示截斷, 取1表示未發(fā)生截斷； 存活超過實驗截止時間表示刪失, 刪失比為61.8%, 示性函數(shù)δ取0表示刪失, 取1表示未發(fā)生刪失. 該數(shù)據(jù)集中有女性365名, 男性97名, 用協(xié)變量x表示,x=1表示女性,x=2表示男性.

在帶有尺度參數(shù)回歸時, 可得生存函數(shù)為

S(t;α,β0,β1)=1-(1-exp{-exp{β0+β1x}t})α,

風險函數(shù)為

當尺度參數(shù)不受協(xié)變量影響時, 利用廣義指數(shù)參數(shù)模型, 并根據(jù)性別分組分析, 結果列于表4. 由表4可見, 男性和女性的生存函數(shù)形狀參數(shù)、尺度參數(shù)較接近, 男性和女性的生存函數(shù)基本相同.

表4 當尺度參數(shù)不受協(xié)變量影響時, Channing house數(shù)據(jù)分析結果

當尺度參數(shù)受協(xié)變量影響時, 利用廣義指數(shù)尺度參數(shù)回歸模型將性別因素作為協(xié)變量處理, 結果列于表5. 由表5可見, 截距項β0=1.865 37, 協(xié)變量性別的回歸系數(shù)β1=0.012 94, 計算協(xié)變量β1對應的p值為0.206 34, 表明性別因素對生存概率的影響不顯著, 性別因素對加州帕洛阿爾托市退休中心老年居民生存概率的影響不明顯.

表5 當尺度參數(shù)受協(xié)變量影響時, Channing house數(shù)據(jù)分析結果