趙遠英,徐登可,冉慶

(1.貴陽學(xué)院數(shù)學(xué)與信息科學(xué)學(xué)院,貴州 貴陽550005;2.浙江農(nóng)林大學(xué)統(tǒng)計系,浙江 杭州311300;3.貴陽市第七中學(xué)數(shù)學(xué)組,貴州 貴陽550001)

1.引言

計數(shù)數(shù)據(jù)[1]廣泛存在于自然科學(xué)領(lǐng)域和社會科學(xué)領(lǐng)域,泊松回歸模型是刻畫計數(shù)型響應(yīng)變量與協(xié)變量之間關(guān)系常用的統(tǒng)計工具,但在實際應(yīng)用時由于計數(shù)數(shù)據(jù)的復(fù)雜性,數(shù)據(jù)和泊松回歸模型經(jīng)常不相吻合,偏大離差或偏小離差等情形時常發(fā)生,這與泊松回歸模型響應(yīng)變量的條件均值與方差相等這一基本要求相違背.因此許多其他的計數(shù)層次回歸模型被提出來解決上述不足.例如:負二項回歸模型[2],泊松對數(shù)正態(tài)模型[3],混合泊松回歸模型[4],廣義泊松回歸模型[5]等.特別地,XIE和WEI[6]基于EM算法研究泊松逆高斯回歸模型的影響診斷問題;解鋒昌等[7]比較系統(tǒng)全面地闡述零膨脹數(shù)據(jù)的統(tǒng)計方法和應(yīng)用價值.有關(guān)計數(shù)數(shù)據(jù)建模及其其它更多統(tǒng)計方法,請參考LIU,TANG和XU[8]等文獻.雖然在計數(shù)層次回歸模型及其方法方面已經(jīng)取得豐碩的研究成果,然而鮮有文獻研究泊松逆高斯回歸模型在貝葉斯框架下的統(tǒng)計推斷問題.因此,本文旨在貝葉斯統(tǒng)計和泊松逆高斯回歸模型的框架背景下:1)估計模型的未知參數(shù);2)評估模型的擬合優(yōu)度.

文章的結(jié)構(gòu)見下:第2部分將泊松逆高斯分布看作逆高斯分布的泊松混合,定義泊松逆高斯回歸模型;第3部分將服從逆高斯分布的隨機變量視為缺失數(shù)據(jù),使用數(shù)據(jù)添加的思想并基于Gibbs抽樣[9]和Metropolis-Hastings算法[10-11]以及Multiple-Try Metropolis算法[12]詳細描述模型未知參數(shù)的貝葉斯估計與貝葉斯擬合優(yōu)度統(tǒng)計量的計算過程,處理算法中包含的高維積分問題;第4部分是模擬研究和實證分析;本文的結(jié)論在第5部分;附錄部分是抽樣算法的具體細節(jié).

2.泊松逆高斯回歸模型及其假設(shè)

與文[13-14]類似,本文關(guān)注泊松逆高斯回歸模型:

其中符號p(νi|λ)～IG(λ)表示隨機變量νi服從逆高斯分布IG(λ),其概率密度函數(shù)為

這里I(·)是示性函數(shù).因此E(νi)= 1,Var(νi)=λ1.而p(yi|νi,μi)～Pois(νiμi)表示在νi和μi給定的條件下,響應(yīng)變量yi服從參數(shù)為νiμi的泊松分布,即E(yi|νi,μi)= Var(yi|νi,μi)=νiμi.在模型(2.1)中,log(μi)=xTi β表明參數(shù)log(μi)與p×1協(xié)變量xi存在線性關(guān)系,β=(β1,··· ,βp)T是未知參數(shù)向量.上述泊松逆高斯分布可視為逆高斯分布的泊松混合[14].下面將在貝葉斯統(tǒng)計和泊松逆高斯回歸模型的框架下研究模型的參數(shù)估計問題和擬合優(yōu)度檢驗問題.

3.模型的貝葉斯統(tǒng)計推斷

記θ=(βT,λ)T,本文假設(shè)θ的先驗分布為:

且

其中N(β0,Hβ)表示期望為β0,協(xié)方差矩陣為Hβ的正態(tài)分布,Γ(a,b)代表參數(shù)為a與b,期望為的伽瑪分布,β0,Hβ,a,b是給定的超參數(shù).

記y={y1,y2,··· ,yn},x={x1,x2,··· ,xn}和潛變量ν={ν1,ν2,··· ,νn},并記觀測數(shù)據(jù)Z={y,x}.對參數(shù)θ的貝葉斯統(tǒng)計推斷可以基于θ關(guān)于觀測數(shù)據(jù)Z的后驗分布p(θ|Z)∝p(Z|θ)p(θ),其中是觀測數(shù)據(jù)Z的似然函數(shù).顯然上述積分無解析表達式.為了處理上述高維積分問題,借助數(shù)據(jù)添加[15]的思想,本文將潛變量ν視為缺失數(shù)據(jù),考慮參數(shù)θ基于完全數(shù)據(jù){Z,ν}的后驗分布:

(3.3)式的優(yōu)勢在于不涉及高維積分的計算,但是直接從(3.3)式中抽樣是不可能實現(xiàn)的.因此在這里借助Gibbs抽樣[9]方法獲得分布p(θ,ν|Z)的隨機樣本(θ,ν),然后基于此隨機樣本對參數(shù)θ作相應(yīng)的貝葉斯統(tǒng)計推斷.比如,參數(shù)g(θ)的貝葉斯后驗期望估計為

步1初始化β與λ以及ν,使它們的初始值分別為β(0)和α(0)以及ν(0)= (ν1(0),··· ,νn(0))T,并令t=0;

步2從p(β|y,ν(t),x,λ(t))中抽樣β(t+1);這里

步3 從p(λ|y,ν(t),x,β(t+1))中抽樣λ(t+1);其中即

步4?i=1,··· ,n,從p(νi|yi,xi,β(t+1),λ(t+1))中抽樣ν(t+1)i;這里

步5 令t=t+1,循環(huán)執(zhí)行步2到步4直至算法收斂后,獲得需要的隨機樣本.

本文在抽樣的過程中采用Metropolis-Hastings算法[10-11]以及Multiple-Try Metropolis算法[12],細節(jié)見附錄部分.

類似于文[16-17],可得未知參數(shù)θ與潛變量ν的聯(lián)合貝葉斯估計為:

其中{(θ(t),ν(t)):t=1,··· ,T}是Gibbs抽樣方法收斂后從p(θ,ν|Z)中獲得的樣本,且由(3.8)式表示的貝葉斯估計是其相應(yīng)后驗均值的相合估計[18].

為了評估貝葉斯統(tǒng)計框架下所研究模型是否合理,后驗預(yù)測p值[19]和偏后驗預(yù)測p值[20]是評估模型擬合優(yōu)度檢驗的常用統(tǒng)計量.本文考慮以下的假設(shè)檢驗問題H0:數(shù)據(jù){y,x}來自于由(2.1)式定義的泊松逆高斯回歸模型;H1:H0不成立.類似于Gelman,MENG和Stern[21]的研究,定義泊松逆高斯回歸模型的后驗預(yù)測p值為

這里yrep表示y的一個重抽樣值,Pr代表在y,x已知和H0成立的條件下,關(guān)于(θ,ν)的聯(lián)合后驗分布求條件概率,其表達式為

偏差度函數(shù)D(·|·)取為

其中{(θ(t),ν(t)):t= 1,...,T}是產(chǎn)生于P(θ,ν|y,x,H0)的隨機樣本,這些隨機樣本與(3.8)式的樣本相同.值得注意的是,(3.11)式中的散度函數(shù)D(y|θ(t),x,ν(t)))很容易計算,而χ2分布的分位數(shù)可以使用通常的統(tǒng)計軟件獲得,故(3.11)式的計算很容易完成.與文[16]相類似,當(dāng)∈(0.3,0.7)時,逆高斯回歸模型(2.1)對數(shù)據(jù)的擬合效果很好;否則,模型(2.1)對數(shù)據(jù)的擬合效果不理想.

類似于文[20-21]的研究,將模型(2.1)中的ν視為缺失數(shù)據(jù),定義泊松逆高斯回歸模型的偏后驗預(yù)測p值為(3.12)式表示Pr{D(y|θ,x,ν)≥D(yobs|θ,x,ν)}關(guān)于分布p?(θ,ν)求期望,其中yobs為響應(yīng)變量y的觀測值,

Dobs是D(·|·)基于yobs的觀察值.易知在給定θ,x,ν的條件下,Dobs近似服從自由度為n的χ2分布.因此(Dobs|θ,x,ν)條件獨立于θ,ν,從而可得

因為(3.12)式含有難以處理的積分問題,Monte Carlo方法被用來解決這一困難.記{(θ(t),ν(t)):t=1,...,T}是來自p?(θ,ν)的隨機樣本(恰為(3.8)式所述),{y(t):t=1,...,T}是來自p(y|θ(t),x,ν(t))的隨機樣本.則由Monte Carlo積分原理,偏后驗預(yù)測p值pppB的估計為

其中I(.)是示性函數(shù).與文[22]相類似,若的值遠離0.5,則認為H0不真.

4.數(shù)值例子

?i=1,··· ,n,假設(shè)yi產(chǎn)生于(2.1)式定義的泊松逆高斯回歸模型,νi與μi分別由下面的模型產(chǎn)生:

協(xié)變量xi= (xi1,xi2,xi3,xi4,xi5)T的每一分量都分別獨立從均勻分布U(?1,1)中產(chǎn)生,設(shè)置參數(shù)的真值為下面兩種情形:

情形1β=(β1,β2,β3,β4,β5)T=(1,1,1,1,1)T和λ=1;

情形2β=(β1,β2,β3,β4,β5)T=(1,1,1,1,1)T和λ=2.

本文分別研究n=100與n=200兩種情形的模擬數(shù)據(jù)樣本,同時考慮以下兩種類型的先驗分布對貝葉斯估計方法的影響:

Type 1(良好的先驗)β0取為β的真值,Hβ=0.25I5和a=b=1,其中0.25Im表示主對角元素為0.25的m×m對角矩陣.

Type 2(無信息先驗)β0=(0,0,0,0,0)T,Hβ=10I5和a=b=1.

本文從MCMC算法收斂[23]后的Gibbs抽樣中收集4000個隨機樣本用來計算(3.8)式未知參數(shù)和潛變量的聯(lián)合貝葉斯估計.表1和表2是重復(fù)模擬100次的相關(guān)計算結(jié)果,其中Mean,Median,5%th,95%th和SD分別表示100次重復(fù)基礎(chǔ)上得到的貝葉斯估計的均值,中位數(shù),5%和95%樣本分位數(shù)和樣本標(biāo)準(zhǔn)差,RMS代表參數(shù)真值與100次重復(fù)的貝葉斯估計的均方根.根據(jù)表1不難得出:(i)對先驗分布中超參數(shù)的不同選擇,估計結(jié)果都比較穩(wěn)健,且在第一種先驗類型條件下的貝葉斯估計略優(yōu)于第二種類型先驗下的貝葉斯估計;(ii)兩種先驗類型條件下的Mean和Median與參數(shù)真值都非常接近,說明該估計方法具有很高的估計精度;(iii)貝葉斯估計方法的性能伴隨著n的變大越來越穩(wěn)定.

表1 模擬研究的數(shù)值結(jié)果(基于情形1)

實證分析的案例是1987-1988年美國醫(yī)療支出調(diào)查(the National Medical Expenditure Survey)數(shù)據(jù).該數(shù)據(jù)的詳細介紹和描述請參閱文[24].Gómez-Déinz等[25]曾經(jīng)使用下列泊松高斯回歸模型

擬合該數(shù)據(jù),其中響應(yīng)變量yi= HOSPi,協(xié)變量xi1= 1,xi2= EXCLHLTHi,xi3=POORLHLTHi,xi4= PRIVINSi,xi5= MEDICAIDi,變量的具體含義由表3給出.相關(guān)參數(shù)的極大似然估計結(jié)果為= (?1.476,?0.950,0.878,0.138,0.264)T和= 0.139.在貝葉斯統(tǒng)計框架下,本文依然使用模型(4.3)擬合該數(shù)據(jù)并在參數(shù)的先驗分布中選擇如下的超參數(shù):β0= (?1.476,?0.950,0.878,0.138,0.264)T,Hβ= 0.25I5和a=b= 1.在進行Gibbs抽樣時,取σ2β= 3.6和σ2ν= 138,得到相應(yīng)的β和ν的接受概率約為27.8%和32.5%.基于模型(4.3)和上述超參數(shù)設(shè)置,使用參數(shù)不同的初始值產(chǎn)生馬爾科夫鏈計算參數(shù)的EPSR值[23].圖1描繪了相應(yīng)EPSR的計算結(jié)果.根據(jù)圖1可知,當(dāng)?shù)螖?shù)iter numbers=500時,所有參數(shù)的EPSR值都未超過1.2.為了確保收斂,丟棄Gibbs抽樣參數(shù)前1000次迭代的迭代樣本,使用1000次迭代之后的3000個隨機樣本,按照(3.8)式計算模型未知參數(shù)的貝葉斯估計,表4是最終的估計結(jié)果.從表4可得,未知參數(shù)的貝葉斯估計為= (?1.537,?0.948,1.056,0.168,0.275)T和= 0.401.這和Gómez-Déniz等[25]的估計結(jié)果是基本相合的.另外,根據(jù)(3.11)式與(3.14)式并借助Gibbs抽

樣收集的3000個隨機樣本計算擬合優(yōu)度統(tǒng)計量,得到后驗預(yù)測p值和偏后驗預(yù)測p值的估計分別為= 0.513和= 0.595,結(jié)果一致表明所建立的模型(4.3)對美國醫(yī)療支出調(diào)查(the National Medical Expenditure Survey)數(shù)據(jù)的擬合效果比較理想.

表2 模擬研究的數(shù)值結(jié)果(基于情形2)

表3 變量及其含義

圖1 美國醫(yī)療支出調(diào)查數(shù)據(jù)分析所有未知參數(shù)的EPSR值

表4 美國醫(yī)療支出調(diào)查數(shù)據(jù)分析未知參數(shù)的貝葉斯估計

5.結(jié)論

文章在泊松逆高斯回歸模型的框架下,借助Gibbs抽樣與Metropolis-Hastings算法以及Multiple-Try Metropolis算法等MCMC方法進行貝葉斯統(tǒng)計推斷研究.數(shù)值例子的結(jié)果表明,所描述的方法能有效地估計模型的未知參數(shù)和評價模型的擬合優(yōu)度.

6.附錄

根據(jù)3.2節(jié)的敘述,(3.6)式的分布是可以直接進行抽樣的伽瑪分布,但是(3.5)式與(3.7)式分布的抽樣就十分棘手和困難.

對如何從(3.5)式的條件分布抽樣β(t+1),本文采用Metropolis-Hastings算法[10-11].首先選取建議分布為N(0,σ2βΣβ),其中假設(shè)β在當(dāng)前的迭代值為β(t),然后從N(β(t),σ2βΣβ)中隨機產(chǎn)生潛在的轉(zhuǎn)移樣本β?,并同時從均勻分布U(0,1)中獨立地產(chǎn)生一個隨機樣本u,若則令β(t+1)=β?,否則令β(t+1)=β(t).在具體使用MH算法時,經(jīng)常選擇σ2β使得β的接受概率位于區(qū)間[0.25,0.45].

Multiple-Try Metropolis算法[12]是一種重要的MCMC算法,被使用來對(3.7)式的條件分布p(νi|yi,xi,β,λ)進行抽樣ν(t+1)i,其步驟為:1)從建議分布N(ν(t)i ,σ2ν)中獨立地產(chǎn)生k個試驗樣本c1,c2,··· ,ck(比如k= 10,50等,本文取k= 10已足夠);2)對任意j= 1,2,··· ,k,用式子πj=p(νi=cj|yi,xi,β,λ)計算πj,并做正則化處理,即令然后從試驗樣本c1,c2,··· ,ck中以概率p1,p2,··· ,pk隨機挑選一個潛在轉(zhuǎn)移樣本ν?,即從離散分布

中隨機產(chǎn)生一個樣本ν?;3)基于ν?,從建議分布N(ν?,σ2ν)中獨立地產(chǎn)生k?1個參考樣本c?1,c?2,··· ,c?k?1,并令c?k=ν(t)i.對任意j= 1,2,··· ,k,用式子π?j=p(νi=c?j|yi,xi,β,λ)計算π?j;4)獨立地從均勻分布U(0,1)中產(chǎn)生一個隨機樣本ε,若