左截斷右刪失數(shù)據(jù)下伽馬分布的參數(shù)推斷

周旭，田茂再，2

(1. 新疆財經(jīng)大學(xué) 統(tǒng)計與數(shù)據(jù)科學(xué)學(xué)院，新疆 烏魯木齊 830012；2. 中國人民大學(xué) 統(tǒng)計學(xué)院，北京100872)

在生存分析的實(shí)際研究中,數(shù)據(jù)常伴隨截斷和刪失。數(shù)據(jù)的不完整給統(tǒng)計推斷過程帶來困難,若直接忽略截斷和刪失數(shù)據(jù)所包含的信息,往往會導(dǎo)致估計結(jié)果不盡人意。近幾年對于截斷和刪失同時存在的數(shù)據(jù)類型進(jìn)行壽命分布的參數(shù)推斷變得越來越普遍,在截斷和刪失同時存在的數(shù)據(jù)中,參數(shù)估計使用MLE(極大似然估計)方法是很困難的,而且由于伽馬函數(shù)的存在,對不完全伽馬函數(shù)的對數(shù)求最大化并不簡單,矩估計是一種數(shù)值方法,且矩估計量具有很好的一致性和漸近正態(tài)性,能夠引入潛在數(shù)據(jù)對參數(shù)進(jìn)行有效估計。伽馬分布是統(tǒng)計學(xué)中一種應(yīng)用較為廣泛的連續(xù)型壽命分布,在概率論中與眾多分布有著密切的聯(lián)系。在此背景下,選擇對左截斷右刪失數(shù)據(jù)下伽馬分布的參數(shù)估計問題進(jìn)行討論具有重要的現(xiàn)實(shí)意義。

1 文獻(xiàn)綜述

生存分析中生存時間數(shù)據(jù)常常是截斷和刪失的,例如,在對某種類型的設(shè)備使用壽命進(jìn)行的調(diào)查中發(fā)現(xiàn),調(diào)查對象通常是使用壽命超過規(guī)定閾值的設(shè)備,而在調(diào)查前結(jié)束使用的設(shè)備和調(diào)查結(jié)束時仍在繼續(xù)使用的設(shè)備并沒有參與調(diào)查,得不到這些設(shè)備的具體使用壽命。因此,得到的樣本是不完全的,即存在左截斷和右刪失,若僅用不完全樣本直接進(jìn)行分析顯然不合適,故對左截斷右刪失數(shù)據(jù)進(jìn)行特殊處理是十分必要的。

近幾年,在截斷和刪失數(shù)據(jù)的基礎(chǔ)上進(jìn)行壽命分布的參數(shù)推斷越來越普遍。例如,對于參數(shù)推斷,Casella等[1]提出極大似然估計法是目前最流行的估計方法;Ahmadi等[2]假設(shè)在給定左截斷右刪失數(shù)據(jù)和生存時間隨機(jī)變量時,比例風(fēng)險模型會被簡化,這種情況下能夠得到未知參數(shù)的MLE和UMVUE(一致最小方差無偏估計),并詳細(xì)給出了指數(shù)分布的點(diǎn)估計;胡江山等[3]研究了左截斷右刪失數(shù)據(jù)下泊松分布的貝葉斯估計,主要比較了極大似然估計和貝葉斯估計的有效性,證明了在小樣本的情況下,貝葉斯估計比極大似然估計的精度要高,但在大樣本時兩者的精度相差不大,這也說明極大似然估計方法面對不完全數(shù)據(jù)的處理還是存在缺陷,故而Dempster等[4]針對不完全數(shù)據(jù)提出了EM(期望最大化)算法,通過引入潛在數(shù)據(jù)有效地進(jìn)行統(tǒng)計推斷。Shang等[5]在左截斷右刪失數(shù)據(jù)下進(jìn)行廣義伽馬分布的估計時,基于EM算法改進(jìn)后提出了EM算法的SEM(隨機(jī)版本)作為計算近似極大似然估計的替代方法,并通過兩種不同的方法獲得迭代估計過程中的初始估計值,用Bootstrap(自助抽樣)方法對廣義伽馬分布進(jìn)行區(qū)間估計,最后運(yùn)用蒙特卡洛模擬評估參數(shù)估計的優(yōu)劣性;胡雋等[6]基于左截斷右刪失數(shù)據(jù)利用EM算法對指數(shù)分布進(jìn)行參數(shù)估計,并通過實(shí)際數(shù)據(jù)證明所建立的似然函數(shù)能夠得到收斂的參數(shù)估計值。此外,矩估計也是一種數(shù)值方法,具有良好的統(tǒng)計性質(zhì),在不完全數(shù)據(jù)類型下也有許多應(yīng)用,如陳超等[7]提出使用矩估計的思想將參數(shù)估計問題轉(zhuǎn)變成方程求解的問題,并證明了矩估計量的一致性和漸近正態(tài)性;Saulo等[8]給出雙參數(shù)伽馬分布的極大似然方法和矩估計量,并對三參數(shù)伽馬分布模型的極大似然估計的二階漸進(jìn)性進(jìn)行討論;梁遠(yuǎn)勝等[9]針對非隨機(jī)缺失的數(shù)據(jù)進(jìn)行研究,使用矩估計方法對兩參數(shù)伽馬分布模型進(jìn)行參數(shù)估計和分位數(shù)估計,也得到了很好的結(jié)果,但沒有對伽馬分布參數(shù)的置信區(qū)間做進(jìn)一步討論。

此外,關(guān)于左截斷右刪失數(shù)據(jù)下構(gòu)造置信區(qū)間估計的研究中,Sauter等[10]提出基于參數(shù)向量MLE的漸近正態(tài)分布計算每個參數(shù)的近似置信區(qū)間;白永昕等[11]探究了左截斷右刪失數(shù)據(jù)中小樣本時剩余壽命分位數(shù)的置信區(qū)間,并與傳統(tǒng)構(gòu)造置信區(qū)間的方法進(jìn)行比較;Chen等[12]提出了伽馬分布參數(shù)的閉式估計,其隨著樣本量的增加呈漸近正態(tài)分布,但這些估計的性能非常接近MLE。

鑒于此,基于左截斷右刪失數(shù)據(jù)選擇矩估計方法對伽馬分布參數(shù)進(jìn)行推斷,通過證明矩估計方法的漸近正態(tài)性質(zhì)得到參數(shù)的漸近置信區(qū)間,并利用蒙特卡洛模擬對置信區(qū)間進(jìn)行評價。

2 左截斷右刪失數(shù)據(jù)標(biāo)識

壽命分析中左截斷和右刪失情況同時存在是常見的數(shù)據(jù)類型,忽略數(shù)據(jù)中的截斷和刪失會導(dǎo)致抽樣誤差,所以對于這類數(shù)據(jù)要進(jìn)行一些特殊處理。左截斷右刪失數(shù)據(jù)的具體含義:在絕大多數(shù)生存分析的實(shí)驗(yàn)研究中,個體在進(jìn)入研究前已經(jīng)被診斷出患病稱為左截斷;在研究結(jié)束時感興趣事件還沒有發(fā)生或者在研究結(jié)束之前由于其他原因退出觀察,這種情況被稱為右刪失。例如在某種疾病對個體生存時間影響的研究中,假設(shè)跟蹤觀察從2011年開始至2022年結(jié)束,共12年,感興趣事件為在跟蹤觀察期間首次診斷出患有該疾病的個體到死亡或研究結(jié)束的生存時間,具體情況如圖1所示。

患者1感興趣事件在第2年初首次診斷出患病,第6年初死亡,其生存時間是可以正常觀測到的;患者2感興趣事件在第3年初首次診斷出患病,第9年初失蹤,存在右刪失;患者3感興趣事件在第7年初首次診斷出患病,第12年初仍活著,存在右刪失;患者4感興趣事件在進(jìn)入研究前就已經(jīng)被診斷出患病,第12年后死亡,存在左截斷右刪失。這4位患者的生存時間分別為4,6+,5+,12+(年)。

3 左截斷右刪失數(shù)據(jù)下伽馬分布的參數(shù)推斷

3.1 基于矩估計左截斷右刪失數(shù)據(jù)下伽馬分布的點(diǎn)估計

有h個服從Gamma(α,β)的樣本Y=(y1,y2,…,yh)T,每個樣本之間獨(dú)立同分布,其密度函數(shù)和分布函數(shù)如下:

(1)

矩估計的主要思想是用樣本的一階和二階原點(diǎn)矩去估計總體的一階和二階原點(diǎn)矩,得到參數(shù)α和β的估計值。將完全樣本X=(x1,x2,…,xn)T中能觀測到的數(shù)據(jù)記為Y=(y1,y2,…,yh)(h

f(xi;α,β)=p·f(xi;α,β|δi=1)+

(1-p)·f(xi;α,β|δi=0),

其中p=P(δi=0)表示刪失數(shù)據(jù)的比例。樣本和總體的一階原點(diǎn)矩和二階原點(diǎn)矩分別為:

(2)

如果數(shù)據(jù)不存在截斷和刪失,則得到參數(shù)數(shù)α和β的估計分別為:

(3)

對來自條件分布的f(xi;α,β|δi=0)的右刪失數(shù)據(jù)添加新的潛在數(shù)據(jù)Z=(z1,z2,…,zl)T,將觀測數(shù)據(jù)和潛在數(shù)據(jù)整理成完整數(shù)據(jù)X=(Y,Z),用完整數(shù)據(jù)的樣本一階和二階原點(diǎn)矩估計總體一階原點(diǎn)矩μ1和二階原點(diǎn)矩μ2:

(4)

假設(shè)觀測到y(tǒng)i的概率(刪失機(jī)制)服從指數(shù)分布f(yi;λ)=1-exp(-yi/λ),λ是指數(shù)分布的參數(shù),刪失機(jī)制也可服從其他分布,但在指數(shù)分布下對于刪失比例的調(diào)整更加靈活,應(yīng)用到實(shí)際數(shù)據(jù)中適用性更強(qiáng)。記沒有觀測到y(tǒng)i的次數(shù)為η(yi),由于η(yi)是個隨機(jī)變量,所以用期望代替。在給定觀測值yi的條件下,η(yi)服從平移的幾何分布:

f(η(yi)=s|yi)=(1-f(yi,λ))sf(yi,λ),

(5)

則有η(yi)期望為

E(η(yi)|yi)=(1-f(yi,λ))/f(yi,λ)。

對于觀測數(shù)據(jù)Y=(y1,y2,…,yh)T(h

f(zi,α,β|δi=0)=

p=(λ/(β+λ))α。

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

3.2 基于矩估計左截斷右刪失數(shù)據(jù)下伽馬分布參數(shù)的漸近置信區(qū)間

矩估計漸近正態(tài)性的證明是基于中心極限定理和Δ方法的。

定理1[13](中心極限定理) 設(shè)X1,X2,…,Xn是一組n維獨(dú)立同分布的隨機(jī)變量,記μ=EXi,Σ=covXiXj,i,j=1,2,…,n,i≠j,且Σ是正定的,則

(14)

進(jìn)一步假設(shè)函數(shù)h(·)在點(diǎn)μ∈k的鄰域內(nèi)具有連續(xù)偏導(dǎo)數(shù),以及Σ正定,用c表示h(·)在μ處的梯度向量:

則

(15)

(16)

(17)

式中:

證明假定θ=(α,β)T,則基于n個服從Gamma(α,β)的左截斷右刪失數(shù)據(jù)X的樣本和總體的一階、二階原點(diǎn)矩分別為:

(18)

(19)

(20)

且Σ正定,通過定理2得證。

(21)

(22)

式中:za/2是標(biāo)準(zhǔn)正態(tài)分布上側(cè)面積為a/2的z值。

4 數(shù)值模擬

為評價使用矩估計對左截斷右刪失數(shù)據(jù)下伽馬分布的參數(shù)進(jìn)行推斷的合理性,基于上述理論依據(jù)式(9)—式(13)和式(21)—式(22),進(jìn)行數(shù)值模擬試驗(yàn)。通過(6)式計算得到成功觀測yi的概率λ,設(shè)定隨機(jī)變量Y服從Gamma(3,0.05),X服從Gamma(2,0.05),刪失比例p為0.1、0.2、0.3,完全數(shù)據(jù)的樣本量n為20、50、100、300,將不同刪失比例和樣本量組合進(jìn)行模擬,對于參數(shù)α和β、λ的初始值分別選擇9、6、1,在迭代時兩次參數(shù)估計值的差值小于設(shè)定的閾值0.000 1則停止迭代。在編寫R語言代碼時主要用到的函數(shù)有rgamma( )以及一些自定義函數(shù)。數(shù)值模擬結(jié)果見表1。

表1 不同樣本量、刪失比例情況的參數(shù)估計模擬結(jié)果

(a)n=20 (b)n=50

5 結(jié)論

針對左截斷右刪失數(shù)據(jù),在固定截斷和刪失時間點(diǎn)的情況下討論了伽馬分布的參數(shù)推斷問題。引入潛在數(shù)據(jù)后利用矩估計方法建立完全樣本一階、二階原點(diǎn)矩與總體一階、二階原點(diǎn)矩之間的關(guān)系,并通過證明矩估計方法的漸近正態(tài)性質(zhì)得到伽馬分布參數(shù)的漸近置信區(qū)間。數(shù)值模擬試驗(yàn)結(jié)果表明:

1) 固定完全數(shù)據(jù)的樣本量時,通過比較不同刪失比例下置信區(qū)間的長度和覆蓋率發(fā)現(xiàn):刪失比例越小,需要引入潛在數(shù)據(jù)的個數(shù)就越少,參數(shù)估計結(jié)果會更接近參數(shù)真值,且置信區(qū)間長度也較短。

2) 固定刪失比例時,通過比較不同完全數(shù)據(jù)的樣本量下置信區(qū)間的長度和覆蓋率發(fā)現(xiàn):樣本量越大,需要引入潛在數(shù)據(jù)的個數(shù)較多,由于矩估計方法適用于大樣本情況,所以參數(shù)估計結(jié)果會更接近參數(shù)真值,且置信區(qū)間長度也較短。