基于零膨脹廣義Poisson回歸模型的廣義估計方程方法及其應用

殷明娥， 于 洋

(遼寧師范大學 數(shù)學學院，遼寧 大連 116029)

自Lambert[1](1992)開發(fā)了零膨脹Poisson回歸模型，將Poisson回歸模型與零退化分布相結合分析產品在生產過程中所包含有的瑕疵數(shù)量之后，零膨脹模型在經濟科學、金融保險等領域有了廣泛應用.Gupta等[2](2004)和Famoye等[3](2006)推廣為零膨脹廣義Poisson回歸模型，對過多零、過度分散或壓縮的計數(shù)數(shù)據進行建模，使模型能夠適用于更多的計數(shù)數(shù)據.但是，對于縱向觀測數(shù)據，當來自同一集群的觀測值之間存在相關性時，如果忽略這個因素可能會導致無效的統(tǒng)計推斷.目前，已有多種統(tǒng)計方法被用于縱向數(shù)據分析中，其中，Liang和Zeger[4](1986)提出的廣義估計方程方法是最具有代表性的和應用最廣的方法.它是通過響應變量的一階矩和二階矩刻畫集群水平上響應與協(xié)變量之間的關系，而且只要響應變量的矩被正確地指定，廣義估計方程的估計量具有一致性和漸近正態(tài)性.

對零膨脹特征的縱向計數(shù)數(shù)據建模包括了兩部分，一部分是普通的計數(shù)部分建模，一部分是零膨脹部分的建模，并不能直接應用廣義估計方程方法來建模. Hall和Zhang[5](2004)對一種通用算法的自適應，稱為期望求解(Expectation Solution)算法來估計模型參數(shù). 通過引入隱變量，顯示響應變量是來自零退化分布還是來自標準的計數(shù)的分布，將廣義估計方程的方法推廣到零膨脹計數(shù)模型的應用中.

Kong等[6](2015)考慮了一個零膨脹負二項模型對零膨脹縱向計數(shù)數(shù)據進行統(tǒng)計推斷，模型分別對零膨脹部分與標準負二項模型部分區(qū)分建模，在Hall和Zhang[5](2004)方法基礎上，允許指定了一個工作相關矩陣和邊際關系. 此后，Hyoyoung和Datta[7](2018)將廣義估計方程方法用于零膨脹Conway Maxwell Poisson模型中, 可以靈活處理數(shù)據散度和縱向數(shù)據的集群性特征，并分析了基因數(shù)據.

在本研究中，參考了Hall和Zhang[5](2004)的思路，基于零膨脹廣義Poisson回歸模型擬合具有零膨脹特征的縱向計數(shù)數(shù)據，既考慮了數(shù)據過度分散和零過多的計數(shù)性質, 同時也考慮了縱向數(shù)據個體內部間的相關性，并且分析了一組Iowa氟化物研究數(shù)據, 探討氟化物與齲齒及其影響因素之間的關系.

1 零膨脹廣義Poisson回歸模型與估計方法

1.1 廣義Poisson分布

假設隨機變量Wij～fGP(wij|λij,φ),i=1,…,N,j=1,…,ni，根據Consul和Jain[8](1970)定義的廣義Poisson分布，概率分布律為

P(Wij=wij)=fGP(wij|λij,φ)=

(1)

當參數(shù)φ<1時，m是使λij+m(φ-1)>0的最大正整數(shù)值，且Wij的均值和方差分別為

E(Wij|λij,φ)=λij,

Var(Wij|λij,φ)=φ2λij.

當參數(shù)φ=1時，式(1)退化為Poisson分布；

當參數(shù)φ>1時，其均值要小于方差，此時數(shù)據中會出現(xiàn)離差偏大；

1.2 零膨脹廣義Poisson回歸模型

如果隨機變量Yij服從零膨脹廣義Poisson分布，記為Yij～ZIGP(λij,pij,φ)，λij、pij、φ為參數(shù)，其概率分布律為

則其均值和方差分別為

E(Yij)=(1-pij)λij,

Var(Yij)=(1-pij)λij(φ2+pijλij).

零膨脹廣義Poisson回歸模型與廣義線性模型的構成類似，模型結構如下：

(ⅰ) 模型的響應變量為Yij，Yij～ZIGP(λij,pij,φ),且Yij相互獨立，i=1,…,N,j=1,…,ni.

1.3 模型參數(shù)的估計

零膨脹廣義Poisson回歸模型中回歸系數(shù)β和γ滿足的廣義估計方程為

(2)

其中,Yi=(Yi1,Yi2,…,Yini)T，μi=E(Yi)=(μi1,…,μini)T，

由Hall和Zhang[5](2004)知，由于參數(shù)β和γ分別出現(xiàn)在零膨脹模型的兩個部分，不能直接用廣義估計方程(2)估計，只能通過引進隱變量uij,i=1,…,N,j=1,…,ni分別估計回歸系數(shù)β和γ.

γ的廣義估計方程為

(3)

則

同理，β的廣義估計方程為

(4)

注意到，由廣義估計方程(4)中矩陣diag(1-ui)可知，只是來自廣義Poisson分布中的觀測值Yij參與廣義估計方程(4)的計算，而來自零退化分布的觀測值Yij不參與其中的計算.

若隱變量uij給定，γ和β的估計值可用常規(guī)的Fisher得分法通過迭代求解估計方程(3)和方程(4)獲得.但是在每次迭代時，需要估計uij的值. 在假設給定y,γ,β,φ的條件下，用EM算法獲得uij的估計.

假設用γ(b)，β(b)，φ(b)分別表示γ,β,φ第b次迭代時的估計值，uij可由如下式子迭代，得到估計值.

(5)

其中,

γ和β可由如下迭代式子得到

(6)

(7)

散度參數(shù)φ的估計:

λijφ2=vij.

故當Yij服從ZIGP分布時，參數(shù)φ只與Yij的方差函數(shù)有關，與均值函數(shù)無關，因此，通過求解以下估計方程估計參數(shù)φ，

(8)

其中,

Hi=diag((1-ui1)2,…,(1-uini)2).

由式(8)可得如下方程

(9)

由vij(φ)=φ2λij可得

則

(10)

相關參數(shù)α1，α2的估計：

估計α1的值，令

(11)

且

Uγi=(Ui12,Ui13,…,Uini-1,ni)T，ργi(α1)=E(Uγi)=(ρi12,ρi13,…,ρini-1,ni)T.

因此，α1可由如下估計方程估計

(12)

其中,

Wγi是Uγi的協(xié)方差矩陣.如果將Wγi設定為單位矩陣，并且假設Rγi(α1)是等相關可交換矩陣，α1常見的一個估計可以通過下式估計，

(13)

(14)

同理，對于α2的估計值，令

ρβi(α2)=E(Uβi)=(ρi12,ρi13,…,ρini-1,ni)T.

由于Yi服從廣義Poisson分布，則α2可由如下方程估計

(15)

其中,

且

Hβi=diag{(1-ui1)(1-ui2),…,(1-uini-1)(1-uini)}.

若假設Wβi是單位矩陣，Rβi(α2)是等相關可交換矩陣，α2的估計為

(16)

標準化α2的估計為

(17)

本文中對α1，α2的估計分別用的是式(14)和式(17).

為獲得參數(shù)β，γ，φ，α1，α2的估計值，需要迭代求解.給定參數(shù)φ，α1，α2的估計值，由式(5)、式(6)和式(7)迭代，分別得到隱變量和參數(shù)β和γ的值. 再作為給定值，由式(10)、式(14)和式(17)得到參數(shù)φ，α1，α2的值. 重復上述過程，直到收斂為止.

1.4 模型參數(shù)的方差估計

Liang和Zeger[4](1986)表明，在正確指定回歸模型的前提下，對于任意選擇的工作相關矩陣，廣義估計方程估計量都具有一致性和漸近正態(tài)性.然而，標準的廣義估計方程估計不含有隱變量，如果忽略隱變量所引起的變化，則回歸參數(shù)估計量的方差會小于真實方差.所以,當估計方程中含有隱變量時，基于響應變量的前二階矩以及回歸參數(shù)的當前估計，期望求解算法將隱變量替換為它的條件均值.

參數(shù)θ的漸近協(xié)方差由以下sandwich形式給出

2 實例分析

本文使用的數(shù)據集是從Iowa氟化物研究數(shù)據中收集的一部分數(shù)據(http://www1.dentistry.uiowa.edu/preventiing-fluori-study)，這是一項研究兒童齲齒的數(shù)量與氟暴露關系的縱向數(shù)據，目的是從飲食和非飲食兩方面說明氟化物來源，識別兒童齲病的危險因素.該研究是在美國Iowa州，從1992—1995年間，每隔1.5—6個月時間，通過牙科醫(yī)院向兒童父母定期發(fā)送問卷，獲取健康和飲食方面的數(shù)據，包括氟化物攝入量數(shù)據. 再通過牙醫(yī)對所調查兒童牙齒健康和齲齒的程度的診斷，對每顆牙齒進行0、1或2的評分，然后將這些分數(shù)相加得到每個兒童齲齒數(shù)量評分(CES)，將CES作為計數(shù)數(shù)據進行分析.顯然，0在數(shù)據集中占主導地位，過多的0提供了零膨脹的證據. 由于兒童不同牙齒的CES有潛在的相關性，一個基于廣義估計方程方法分析齲齒是否與不同的暴露和治療有關是適合的. 模型的響應變量是齲齒數(shù)量評分(CES)，協(xié)變量為受試者牙齒檢查時的年齡(Age)、受試者從飲食中攝入的氟化物(AUC)、受試者每次觀察時間點前6個月牙醫(yī)就診次數(shù)比例(Dental Visit)、受試者每次觀察時間點前6個月牙齒接受專業(yè)氟化物治療平均次數(shù)比例(Flouride Treatment)、受試者平均每天刷牙頻率(Brushing)，去除變量缺失的觀察結果后，研究子集包含352名兒童，總計4 353個觀測值.

通過廣義估計方程方法，對零膨脹廣義Poisson回歸模型中廣義Poisson計數(shù)模型和零膨脹的logit模型分別建模，兩個回歸模型的協(xié)變量分別都使用上述的5個協(xié)變量，回歸系數(shù)(β，γ)的估計值、標準誤及Wald檢驗的P值見表1的上半部分.模型的其他參數(shù)(φ，α1，α2)估計結果見表1的下半部分.

表1 基于估計方程方法分析數(shù)據的結果

在表1中計數(shù)模型部分的結果顯示，協(xié)變量中牙齒檢查時的年齡(Age)、受試者從飲食中攝入的氟化物(AUC)，受試者就診次數(shù)比例(Dental Visit)、受試者接受治療次數(shù)比例(Flouride Treatment)4個因素是顯著的.其中，受試者從飲食中攝入的氟化物(AUC)與齲齒數(shù)量呈負相關，其他因素與牙與齲齒的數(shù)量呈正相關.

在表1中零膨脹部分的結果顯示，受試者每天平均刷牙頻率(Brushing)因素是顯著的，而且與零膨脹部分呈正相關，說明刷牙頻率越大，零計數(shù)就越多，齲齒的數(shù)量就越少.

3 討 論

通過對Iowa州的氟化物研究數(shù)據的實例分析，得到了兒童齲齒的影響因素，對實施兒童的預防保健措施提供了依據.本文所得結果與Kong等[6](2015)研究得到結果相仿，說明數(shù)據擬合的效果較好，說明本文提出模型與估計方法是合適的. 事實上，具有零膨脹特征的縱向計數(shù)數(shù)據很多，采用哪種模型擬合，哪種估計方法，還需要進一步深入探討.本文嘗試的基于零膨脹廣義Poisson模型的廣義估計方程方法，分析實際數(shù)據，為今后類似的數(shù)據建模提供一個新的路徑.