泊松逆高斯回歸模型的貝葉斯統計推斷

趙遠英,徐登可,冉慶

(1.貴陽學院數學與信息科學學院,貴州 貴陽550005;2.浙江農林大學統計系,浙江 杭州311300;3.貴陽市第七中學數學組,貴州 貴陽550001)

1.引言

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

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

2.泊松逆高斯回歸模型及其假設

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

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

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

3.模型的貝葉斯統計推斷

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

且

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

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

(3.3)式的優勢在于不涉及高維積分的計算,但是直接從(3.3)式中抽樣是不可能實現的.因此在這里借助Gibbs抽樣[9]方法獲得分布p(θ,ν|Z)的隨機樣本(θ,ν),然后基于此隨機樣本對參數θ作相應的貝葉斯統計推斷.比如,參數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,循環執行步2到步4直至算法收斂后,獲得需要的隨機樣本.

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

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

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

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

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

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

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

類似于文[20-21]的研究,將模型(2.1)中的ν視為缺失數據,定義泊松逆高斯回歸模型的偏后驗預測p值為(3.12)式表示Pr{D(y|θ,x,ν)≥D(yobs|θ,x,ν)}關于分布p?(θ,ν)求期望,其中yobs為響應變量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積分原理,偏后驗預測p值pppB的估計為

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

4.數值例子

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

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

情形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兩種情形的模擬數據樣本,同時考慮以下兩種類型的先驗分布對貝葉斯估計方法的影響:

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)式未知參數和潛變量的聯合貝葉斯估計.表1和表2是重復模擬100次的相關計算結果,其中Mean,Median,5%th,95%th和SD分別表示100次重復基礎上得到的貝葉斯估計的均值,中位數,5%和95%樣本分位數和樣本標準差,RMS代表參數真值與100次重復的貝葉斯估計的均方根.根據表1不難得出:(i)對先驗分布中超參數的不同選擇,估計結果都比較穩健,且在第一種先驗類型條件下的貝葉斯估計略優于第二種類型先驗下的貝葉斯估計;(ii)兩種先驗類型條件下的Mean和Median與參數真值都非常接近,說明該估計方法具有很高的估計精度;(iii)貝葉斯估計方法的性能伴隨著n的變大越來越穩定.

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

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

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

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

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

表3 變量及其含義

圖1 美國醫療支出調查數據分析所有未知參數的EPSR值

表4 美國醫療支出調查數據分析未知參數的貝葉斯估計

5.結論

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

6.附錄

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

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

Multiple-Try Metropolis算法[12]是一種重要的MCMC算法,被使用來對(3.7)式的條件分布p(νi|yi,xi,β,λ)進行抽樣ν(t+1)i,其步驟為:1)從建議分布N(ν(t)i ,σ2ν)中獨立地產生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隨機挑選一個潛在轉移樣本ν?,即從離散分布

中隨機產生一個樣本ν?;3)基于ν?,從建議分布N(ν?,σ2ν)中獨立地產生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)中產生一個隨機樣本ε,若