偏斜正態分布下的ZIP層次回歸模型的貝葉斯方法

呂敏紅，閆奕榮

（1.西安航空學院 理學院，西安 710077；2.西北大學 數學學院，西安 710069；3.西安交通大學 經濟與金融學院，西安 710049）

0 引言

技術數據廣泛存在于醫療、生物學、金融保險以及風險控制，擬合計數數據的單用分布主要有泊松分布，二項分布等。但是在實際問題中零觀測的比例遠超過了擬合分布的允許范圍，即存在零膨脹，故零膨脹模型的研究已成為當今國內外的一個熱點問題。

自從Lambert提出了零點膨脹Psisson回歸模型[1]以來，關于具有零膨脹特征的計數數據已經有了多方面的研究，Greene（1994）[2]在Lambert的思想下提出了零膨脹的可加性負二項回歸模型。Fahrmeir和Echavarria（2006）[3]研究了一類零膨脹的可加模型，Xie（2009）[4]系統研究了廣義的Poisson混合效應模型的統計診斷問題，Ghosh（2006）[5]研究了零膨脹回歸的貝葉斯方法，傳統的零膨脹回歸模型是對隨機效應和隨機誤差作正態的假設，但是在實際中正態假設可能會導致無效的統計結論。本文考慮了隨機誤差和隨機效應服從偏斜正態分布的ZIP層次回歸模型的貝葉斯分析問題，最后用一個實例說明該方法的有效性。

1 零膨脹Poisson回歸模型（ZIP）

ZIP分布的基本思想是取值為零的部分和取值為Poisson的部分各占一定的比例構成ZIP混合分布，即：

其中0＜?＜1為零膨脹系數。顯然當?=0時，ZIP分布變為Poisson分布，λ為泊松分布的均值。（1）式的均值和方差分別為：

2 零膨脹Pission層次回歸模型

在實際問題中，數據可能呈現內在關聯或層次結構，為了刻畫數據的這些關系，本文進一步定義層次回歸模型[6]，層次回歸模型綜合了線性回歸和隨機效應模型的優勢。

假設Yij為本文感興趣的響應變量，yij表示第i個群第 j個樣本的觀察數值。i=1，2，…，m， j=1，2，…，n 相對于傳統模型，層次模型可將傳統模型的誤差項分解到與數據相對應水平上。若Yij~ZIP(φij，λij),針對膨脹參數 ?ij與均值參數λij建立如下混合效應模型：

其中，βij與rij分別是協變量xij與zij的回歸系數，

進一步對上層模型考慮線性回歸，并引入隨機效應：

其中，Wij為協矩陣，β與γ為參數向量，ui與vi為隨機效應。（2）式與（3）式合稱為零膨脹Poisson層次回歸模型。

3 偏斜正態分布下零膨脹Poisson層次回歸模型的貝葉斯方法

經典的零膨脹Poisson回歸模型一般都假設隨機誤差及隨機效應都服從正態分布，但是這種假設過于理想化，現實中很多情況下并不滿足，或者說有些數據按照這種假設建立的模型缺乏穩健性。接下來，本文考慮SN-ZIP層次回歸模型。

3.1 SN-ZIP層次回歸模型

n維隨機變量Y服從n元偏斜正態分布，記作Y～SNn(μ，Σ，Δ)，其概率密度函數為：

其中，μ 為均值，Σ 為尺度矩陣，Δ=diag(δ1，δ2，…，δn)為偏度矩陣，?n和Φn分別為標準正態分布下的概率密度和分布函數。特別當 δ=(δ1，δ2，…，δn)=0 時，分布退化成為多元正態分布。為使用方便，進一步寫出（4）式的層次表示[7]：

假設ZIP層次回歸模型中的隨機誤差和隨機效應都服從SN分布，則ZIP層次回歸模型便成為SN-ZIP層次回歸模型。

首先，ZIP層次回歸模型中的隨機誤差服從SN分布，即（2）式中的：

其中 Δk=diag(δ1(k)，δ2(k)，…，δn(k))，k=1，2 。

其次，ZIP層次回歸模型中的隨機效應也服從SN分布，即（3）式中的：

其中 Δu=diag(δu1，δu2，…，δun)，Δv=diag(δv1，δv2，…，δvn)。

式（2）、式（3）、式（6）、式（7）合稱為SN-ZIP層次回歸模型。

3.2 貝葉斯推斷

與似然方法相比，貝葉斯方法綜合了樣本中的先驗信息，對于某些復雜的模型具有特別的靈活性，下面具體研究SN-ZIP層次回歸模型的貝葉斯推斷。

3.2.1 潛變量的數據添加

零膨脹回歸模型中的響應變量Yij可以表示為Yij=Cij(1 -Bij)[5]，其中Bij是具有參數φij的伯努利分布隨機變量，Cij服從參數為λij的Poisson分布，那么給定：

Yij=yij時(Cij，Bij)的聯合條件分布為：

當 yij＞0時，Bij=0，Cij=yij，即：

P(Bij=0，Cij=yij|Yij=yij)=1

當 yij=0時,有兩種情況 Bij=0，Cij=0或 Bij=1，Cij=cij，此時：

3.2.2 先驗分布和參數設定

若 θ=(β，γ，δ(1)，δ(2)，δu，δvσ2(1)，σ2(2)，Σu，Σv)為本文涉及的全體參數，其中 β、γ是本文感興趣的參數，δ(k)=(δ(k)，

1

δ2

(k)，…，δn(k))T，k=1，2。假設 f(θ)為 θ 的先驗密度函數，在后面的貝葉斯推斷中選擇如下的獨立先驗分布，即：

其中Ωk=diag(σ2(k)), β0、γ0δu0、δv0為層次回歸分析的截距項。σ2(k)、Γ(k)、Γu、Γv、ω1(k)、ω2(k)、ψu、ψv，Ru、Rv為超參數，超參數的選取一般通過給定的先驗信息來確定。

3.2.3 模型建立

本文考慮隨機誤差和隨機效應服從偏斜正態分布的ZIP層次回歸模型,利用偏斜正態分布的層次表示方法，即式（5），本文建立如下模型。

第一步：潛變量建模

第二步：回歸系數建模

用貝葉斯的方法，參數θ的后驗分布基于觀測數據是很難直接計算出來的，可以采用Gibbs抽樣和M-H算法[7]，并且借助計算機可以較為簡單的解決上述問題。在抽樣過程中，由于Gibbs抽樣的順序不會影響貝葉斯估計的結果，當樣本收斂后，就會得到感興趣參數的估計值。然后，可以采用Johnson給出的貝葉斯擬合統計量[8]來計算模型對數據的擬合程度。

貝葉斯模型選擇的方法有很多，比如貝葉斯因子，后驗模型概率和后驗預測檢驗等，本文選用BIC作為模型選擇的準則：

4 實例分析

數據來源于Lloyd社記錄的34條船只的5年內發生事故受損的情況。本文對數據進行分析后發現其具有零膨脹特征。為了分析船只種類、建造時間及服務年限對受損情況的影響，本文建立了偏斜正態分布下的ZIP層次回歸模型，計算出參數的后驗均值及MC誤差，具體見表1，其中A1，A2，…,A5表示船舶類型，B1，B2，…,B4表示建造年代，T1、T2表示服務年限。除使用上述模型外，本文還利用一般ZIP回歸模型及ZIP混合效應模型對該數據進行了擬合，通過BIC準則比較了模型的優劣，計算結果見表2。

表1 參數的貝葉斯估計

表2 不同模型的BIC值

由表2可以看出，層次SN-ZIP的BIC值最小，表明用層次SN-ZIP模型對數據進行擬合，擬合程度最高。也是就說對于這組數據，偏態的作用是顯著的，考慮對隨機誤差及隨機效應服從偏斜正態分布比假設兩者服從正態分布要合理。

5 結論

傳統的ZIP層次回歸模型的基本假設是隨機效應和隨機誤差正態分布，然而在實際中，正態假設缺乏穩健性，也可能會導致無效的統計結論。為了精確處理參數估計的問題，本文考慮了隨機誤差和隨機效應服從偏斜正態分布的情況，最后通過實例說明了該方法的有效性。但是缺失數據下的偏斜正態分布還有待進一步研究。