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

(長(zhǎng)春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 長(zhǎng)春 130012)

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

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

1 左截?cái)嘤覄h失數(shù)據(jù)下廣義指數(shù)的參數(shù)估計(jì)

1.1 左截?cái)嘤覄h失數(shù)據(jù)

圖1 左截?cái)嘤覄h失數(shù)據(jù)實(shí)例Fig.1 Example of left truncated and right censored data

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

1.2 廣義指數(shù)分布

廣義指數(shù)分布是由單參數(shù)指數(shù)分布向雙參數(shù)分布的一個(gè)推廣, 一般記作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,

風(fēng)險(xiǎn)函數(shù)為

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

1.3 參數(shù)估計(jì)

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

在左截?cái)嘤覄h失數(shù)據(jù)下, 參數(shù)的似然函數(shù)為

(1)

其中:θ表示未知參數(shù);G(τi)表示截?cái)嗖糠值姆植己瘮?shù). 整理可得

(2)

根據(jù)左截?cái)嘤覄h失數(shù)據(jù)下常見(jiàn)分布形式截?cái)嗖糠址植己瘮?shù)的取法[16-18,20], 取截?cái)帱c(diǎn)的分布與興趣時(shí)間一致為廣義指數(shù)分布, 即G(τi)=(1-exp{-λτi})α, 所以在左截?cái)嘤覄h失數(shù)據(jù)下, 廣義指數(shù)分布的似然函數(shù)可表示為

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

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

(5)

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

(7)

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

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

2 數(shù)值模擬

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

文獻(xiàn)[14]中生成左截?cái)嘤覄h失隨機(jī)數(shù)的方法, 用R軟件實(shí)現(xiàn)隨機(jī)模擬. 考察不同樣本量n=100,300,500, 不同刪失比p=0.4,0.6和截?cái)啾萹=0.3下的參數(shù)估計(jì)效果.

當(dāng)尺度參數(shù)不受協(xié)變量因素影響時(shí), 假設(shè)興趣時(shí)間變量T服從參數(shù)為θ=(α,λ)的廣義指數(shù)分布, 截?cái)嗖糠謺r(shí)間的分布與興趣時(shí)間變量T的分布一致, 但其截?cái)鄶?shù)受截?cái)啾萹控制. 刪失時(shí)間點(diǎn)分布服從參數(shù)為k的指數(shù)分布, 參數(shù)k的取值由刪失比p的取值決定.

模擬ⅠT服從參數(shù)為α=1.2,λ=0.6的廣義指數(shù)分布. 實(shí)驗(yàn)?zāi)M在不同樣本量n=100,300,500下, 根據(jù)刪失截?cái)喾殖蓛山M, 一組刪失比和截?cái)啾确謩e為0.4,0.3, 另一組刪失比和截?cái)啾确謩e為0.6,0.3, 循環(huán)1 000次, 求出參數(shù)估計(jì)值, 并與真實(shí)參數(shù)值進(jìn)行比較分析, 模擬結(jié)果列于表1.

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

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

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

模擬Ⅱ 形狀參數(shù)α=1.2, 決定尺度參數(shù)的兩個(gè)回歸參數(shù)為β0=0.8,β1=0.6. 實(shí)驗(yàn)?zāi)M在樣本量n=100,300,500、刪失比和截?cái)啾确謩e為0.4,0.3的條件下循環(huán)模擬1 000次, 模擬結(jié)果列于表2.

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

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

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

3 實(shí)例分析

例1電子變壓器壽命數(shù)據(jù)集[18-19]. 該數(shù)據(jù)集共有286個(gè)觀測(cè)值, 正常失效數(shù)據(jù)39個(gè), 發(fā)生刪失的數(shù)據(jù)247個(gè), 發(fā)生截?cái)嗟臄?shù)據(jù)167個(gè), 所以其刪失比和截?cái)啾确謩e為86.36%,58.39%. 用t表示跟蹤的壽命時(shí)長(zhǎng),τ表示觀測(cè)的截?cái)嗖糠謺r(shí)間,δ,v分別表示刪失和截?cái)嗟氖拘宰兞? 發(fā)生刪失或者截?cái)? 則δ或v取0, 反之取1.

考慮在左截?cái)嘤覄h失數(shù)據(jù)下, 廣義指數(shù)分布的參數(shù)模型, 可得變壓器的生存(可靠性)函數(shù)為

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

危險(xiǎn)率函數(shù)為

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

例2美國(guó)Channing house 數(shù)據(jù)集[1,10]. 該數(shù)據(jù)集, 即1964-01至1975-07美國(guó)加州帕洛阿爾托市退休中心462名老年居民的死亡時(shí)間數(shù)據(jù)集. 其中數(shù)據(jù)只有實(shí)際壽命超過(guò)實(shí)驗(yàn)起始點(diǎn)時(shí)才能觀測(cè)、記錄, 因此退休居民進(jìn)入實(shí)驗(yàn)的年齡即為所需截?cái)帱c(diǎn). 進(jìn)入年齡小于840個(gè)月(70 a)的個(gè)體組成的數(shù)據(jù)子集表示發(fā)生截?cái)? 截?cái)啾葹?8.0%, 截?cái)嗖糠钟忙颖硎? 示性函數(shù)v取0表示截?cái)? 取1表示未發(fā)生截?cái)啵?存活超過(guò)實(shí)驗(yàn)截止時(shí)間表示刪失, 刪失比為61.8%, 示性函數(shù)δ取0表示刪失, 取1表示未發(fā)生刪失. 該數(shù)據(jù)集中有女性365名, 男性97名, 用協(xié)變量x表示,x=1表示女性,x=2表示男性.

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

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

風(fēng)險(xiǎn)函數(shù)為

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

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

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

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