帶有約束條件的零一膨脹Poisson回歸模型

劉 影, 華志強(qiáng)

(1. 吉林師范大學(xué) 數(shù)學(xué)學(xué)院, 吉林 四平 136000; 2. 內(nèi)蒙古民族大學(xué) 數(shù)學(xué)學(xué)院, 內(nèi)蒙古 通遼 028043)

計(jì)數(shù)數(shù)據(jù)是實(shí)際生活中常見的一類數(shù)據(jù)類型, 可描述一個(gè)事件發(fā)生的次數(shù), 通常用取值為自然數(shù)的隨機(jī)變量表示. 該類數(shù)據(jù)在臨床醫(yī)學(xué)、 心理學(xué)、 遺傳學(xué)、 環(huán)境學(xué)和金融學(xué)等領(lǐng)域應(yīng)用廣泛, Poisson模型、 二項(xiàng)模型和負(fù)二項(xiàng)模型是處理計(jì)數(shù)數(shù)據(jù)的傳統(tǒng)方法, Poisson回歸模型、 二項(xiàng)回歸模型和負(fù)二項(xiàng)回歸模型常用來刻畫響應(yīng)變量與解釋變量之間的關(guān)系. 目前, 人們已對(duì)上述經(jīng)典模型進(jìn)行了多方面推廣. 在實(shí)際應(yīng)用中, 通常會(huì)遇到0過多的樣本數(shù)據(jù), 即0觀測(cè)的比例遠(yuǎn)超過了擬合分布的允許范圍, 此時(shí), 傳統(tǒng)的Poisson模型、 二項(xiàng)模型和負(fù)二項(xiàng)模型不再適用. 針對(duì)這種情形, 文獻(xiàn)[1]首次提出了零膨脹數(shù)據(jù)； 文獻(xiàn)[2-3]研究了不考慮協(xié)變量的零膨脹Poisson分布； 文獻(xiàn)[4]建立了零膨脹Poisson回歸模型; 文獻(xiàn)[5]給出了零膨脹模型的Bayes方法; 文獻(xiàn)[6]討論了零膨脹數(shù)據(jù)的Bayes分位回歸模型; 文獻(xiàn)[7]給出了零膨脹分布參數(shù)的固定寬度置信區(qū)間的構(gòu)造方法. 目前, 關(guān)于零膨脹計(jì)數(shù)數(shù)據(jù)的研究已有很多結(jié)果, 但在一些實(shí)際問題中, 樣本不僅會(huì)出現(xiàn)0過多的情況, 也會(huì)出現(xiàn)1過多的情況, 例如： 在一段時(shí)間內(nèi), 某地交通事故發(fā)生的次數(shù)、 極端天氣出現(xiàn)的次數(shù)、 某傳染病的發(fā)病次數(shù)、 汽車保險(xiǎn)的索賠次數(shù)等, 通常0和1出現(xiàn)的次數(shù)都比較多. 針對(duì)這種情形, 文獻(xiàn)[8]提出了零一膨脹Poisson分布; 文獻(xiàn)[9]研究了該分布的性質(zhì); 文獻(xiàn)[10]給出了該分布參數(shù)的極大似然估計(jì)和Bayes估計(jì); 文獻(xiàn)[11]研究了零一膨脹Poisson回歸模型及其統(tǒng)計(jì)推斷方法; 文獻(xiàn)[12]討論了零一膨脹Poisson回歸模型的非參數(shù)統(tǒng)計(jì)分析問題. 此外, 利用統(tǒng)計(jì)模型分析實(shí)際問題時(shí)還發(fā)現(xiàn), 除需要知道樣本信息外, 總體分布也會(huì)受某些條件的限制, 例如： 在考慮產(chǎn)品的合格率時(shí), 由經(jīng)驗(yàn)可知合格率的下限； 由于市場(chǎng)原因或者監(jiān)管法規(guī)的影響, 產(chǎn)品的價(jià)格不能超過某值或者低于某值等. 因此, 模型參數(shù)空間受約束的統(tǒng)計(jì)模型已得到廣泛關(guān)注[13-14]. 本文將該思想引入到零一膨脹模型中, 提出一類帶有一般線性約束條件的零一膨脹Poisson回歸模型, 并討論相應(yīng)的參數(shù)估計(jì)問題.

1 模型簡(jiǎn)介

令Z=(Z0,Z1,Z2)T服從多項(xiàng)分布,X服從Poisson分布, 即Z～Multinomial(1;φ0,φ1,φ2),X～Poisson(λ), 且Z與X相互獨(dú)立. 如果

(1)

其中:P(ξ0=0)=1,P(ξ1=1)=1;φ0∈[0,1)和φ1∈[0,1)分別表示0和1的膨脹率,φ2=1-φ0-φ1∈(0,1]. 顯然, 零一膨脹Poisson分布是3個(gè)隨機(jī)變量的混合： 退化為0的隨機(jī)變量、 退化為1的隨機(jī)變量和Poisson隨機(jī)變量. 特別地, 當(dāng)φ0=0時(shí), 零一膨脹Poisson分布即為一膨脹Poisson分布(OIP)； 當(dāng)φ1=0時(shí), 零一膨脹Poisson分布即為零膨脹Poisson分布(ZIP)； 當(dāng)φ0=φ1=0時(shí), 零一膨脹Poisson分布即為普通的Poisson分布.

由定義易知, 零一膨脹Poisson分布的分布函數(shù)為

E(Y)=φ1+φ2λ, Var(Y)=φ1(1-φ1)+φ2λ(1-2φ1)+φ2λ2(1-φ2).

零一膨脹Poisson分布的其他等價(jià)形式和性質(zhì)參見文獻(xiàn)[9].

在實(shí)際應(yīng)用中, 模型參數(shù)通常會(huì)受一些外在因素的影響, 例如： 某人某段時(shí)期內(nèi)到醫(yī)院就診的次數(shù)可能與其年齡和性別有關(guān), 也可能與其居住地環(huán)境有關(guān), 因此本文將協(xié)變量引入到零一膨脹Poisson分布的參數(shù)中, 得到零一膨脹Poisson回歸模型. 本文只對(duì)參數(shù)λ進(jìn)行回歸分析, 即假設(shè)logλ=WTβ, 其中:W=(1,W1,W2,…,Wp)T為協(xié)變量;β=(β0,β1,…,βp)T為回歸系數(shù).

進(jìn)一步, 統(tǒng)計(jì)模型經(jīng)常有約束條件, 例如： 在醫(yī)療保險(xiǎn)產(chǎn)品定價(jià)中, 某些風(fēng)險(xiǎn)類別的費(fèi)率不能過高, 此時(shí)可假設(shè)回歸系數(shù)滿足如下線性約束：

(2)

2 參數(shù)估計(jì)

本文使用極大似然法進(jìn)行參數(shù)估計(jì), 并借助EM(expectation maximization)算法和Newton-Raphson算法求解, 下面給出迭代過程. 令θ=(φ0,φ1,βT)T為待估參數(shù), 設(shè)Y1,Y2…,Yn為一組來自ZOIP(φ0,φ1,λi)的簡(jiǎn)單隨機(jī)樣本, 即各樣本相互獨(dú)立, 且

其中:i=1,2,…,n;

(3)

為使用EM算法, 引入潛變量Z0i和Z1i(i=1,2,…,n),Z0i=1表示觀測(cè)值yi來自結(jié)構(gòu)零(即式(1)的第一部分),Z0i=0表示觀測(cè)值yi來自分布零(即式(1)的第三部分),Z1i=1表示觀測(cè)值yi來自結(jié)構(gòu)一(即式(1)的第二部分),Z1i=0表示觀測(cè)值yi來自結(jié)構(gòu)一(即式(1)的第三部分). 記y=(y1,…,yn),Z0=(Z01,…,Z0n),Z1=(Z11,…,Z1n), 則完整數(shù)據(jù)下的似然函數(shù)為

易得對(duì)數(shù)似然函數(shù)為

因此, 約束條件下零一膨脹Poisson回歸模型的極大似然估計(jì)可表示為

考慮優(yōu)化如下函數(shù)：

(4)

這里,η是一個(gè)較大的正數(shù). 當(dāng)不滿足約束條件Aβ-c≤0時(shí), 式(4)等號(hào)右邊將減去一個(gè)較大的正數(shù), 相當(dāng)于對(duì)目標(biāo)函數(shù)進(jìn)行了“懲罰”, 從而對(duì)參數(shù)的估計(jì)值進(jìn)行了重新優(yōu)化.

其中:

這里

M-步(β): 令Q1(β)關(guān)于β的一階導(dǎo)數(shù)為0, 即

(5)

其中:S=diag(s1,s2,…,sr),si=max{0,αi0β0+αi1β1+…+αipβp-ci};e=(1,1,…,1)T. 由于方程(5)沒有顯示解, 因此本文采用Newton-Raphson算法得到如下近似解：

其中:

M-步(φ0和φ1): 令Q2(φ0,φ1)分別關(guān)于φ0和φ1的一階導(dǎo)數(shù)為0, 即

(6)

3 數(shù)值模擬與應(yīng)用

假設(shè)φ0=0.5,φ1=0.3, 且lnλi=β1Wi1+β2Wi2(i=1,2,…,n), 其中β1=0.5,β2=-0.3, {Wi1,i≥1}和{Wi2,i≥1}均為獨(dú)立同分布隨機(jī)變量序列, 其共同分布都是Bernoulli分布B(1,0.5). 用R語言產(chǎn)生1 000個(gè)零一膨脹Poisson回歸模型的隨機(jī)數(shù), 圖1為其頻率直方圖. 模擬500次取估計(jì)值的平均數(shù), 無約束條件下各參數(shù)的估計(jì)結(jié)果為:φ0=0.484 8，φ1=0.291 0，β1=0.605 4，β2=-0.234 9.

假設(shè)根據(jù)實(shí)際情形, 需對(duì)W1=1,W2=0類別的個(gè)體進(jìn)行限制, 例如： 取A=(1,0),c=0.55, 考慮約束條件Aβ≤c, 即β1≤0.55, 則可得修正后的參數(shù)估計(jì)結(jié)果為:φ0=0.479 4,φ1=0.285 4,β1=0.459 1,β2=-0.197 4.

綜上可見, 本文方法使得參數(shù)的估計(jì)結(jié)果滿足約束條件, 具有一定的可行性和有效性.

選擇一組澳大利亞新南威爾士州水資源委員會(huì)沖浪/健康研究機(jī)構(gòu)提供的耳病發(fā)生次數(shù)數(shù)據(jù)(http://www.statsci.org/data/oz/earinf.html), 該數(shù)據(jù)集給出了1990年287組個(gè)體觀測(cè)數(shù)據(jù). 響應(yīng)變量是每個(gè)觀測(cè)者的耳病發(fā)生次數(shù), 其中發(fā)生次數(shù)為0的人數(shù)占總?cè)藬?shù)的52.6%, 發(fā)生次數(shù)為1的人數(shù)占總?cè)藬?shù)的13.9%. 解釋變量包括游泳頻率、 游泳地點(diǎn)、 觀測(cè)者的年齡和性別. 此時(shí), 式(3)的回歸模型可寫為

logλi=β0+β1Wi1+β2Wi2+β3Wi3+β4Wi4,

其中:Wi1表示游泳頻率,Wi1=1表示觀測(cè)者經(jīng)常游泳,Wi1=0表示觀測(cè)者不經(jīng)常游泳;Wi2表示游泳地點(diǎn),Wi2=1表示觀測(cè)者在海濱游泳,Wi2=0表示觀測(cè)者不在海濱游泳;Wi3表示年齡,Wi3=1表示觀測(cè)者年齡為15～24歲,Wi3=0表示觀測(cè)者年齡為25～29歲;Wi4表示性別,Wi4=1表示觀測(cè)者為女性,Wi4=0表示觀測(cè)者為男性.

耳病發(fā)生次數(shù)的分布如圖2所示. 由圖2粗略可見, 超過1/2觀測(cè)者的耳病發(fā)生次數(shù)為0, 也有較多觀測(cè)者的耳病發(fā)生次數(shù)為1, 表明該數(shù)據(jù)集存在零一膨脹特點(diǎn), 適用于零一膨脹模型. 無約束條件下各參數(shù)的估計(jì)結(jié)果為:φ0=0.505 2，φ1=0.121 4，β0=1.519 3，β1=-0.541 8，β2=-0.076 3，β3=-0.147 7，β4=-0.049 0.

圖1 隨機(jī)數(shù)的頻數(shù)直方圖

圖2 耳病發(fā)生次數(shù)的分布

為說明約束條件的影響, 考慮因?yàn)棣?E(Y|Z3=1), 因而λ可理解為觀測(cè)者發(fā)生耳病的高風(fēng)險(xiǎn)部分, 在審核相應(yīng)的保險(xiǎn)產(chǎn)品費(fèi)率時(shí), 其具有重要作用. 假設(shè)在市場(chǎng)條件下, 行業(yè)政策規(guī)定針對(duì)這部分的費(fèi)率不能過高, 如不允許超過4. 下面考慮在該約束條件下的模型和參數(shù)估計(jì). 記

此時(shí), 可在模型中加入約束條件Aβ≤c, 并得到約束條件下各參數(shù)的估計(jì)結(jié)果:φ0=0.505 5,φ1=0.120 8,β0=1.386 5,β1=-0.455 4,β2=-0.015 1,β3=-0.085 3,β4=-0.000 2.

按照解釋變量W1,W2,W3,W4的不同取值, 觀測(cè)者可被分為16組, 其中第i組解釋變量的取值對(duì)應(yīng)矩陣A第i行中第2列～第5列元素, 例如： 第一組對(duì)應(yīng)W1=0,W2=0,W3=0,W4=0, 第二組對(duì)應(yīng)W1=0,W2=1,W3=0,W4=0, 以此類推. 表1列出了不同組在約束前后的λ預(yù)測(cè)值. 由表1可見, 約束前有4組λ的預(yù)測(cè)值超過4, 而約束條件改變了參數(shù)的估計(jì)結(jié)果, 使得所有組的λ預(yù)測(cè)值都不超過4, 且所有組的λ預(yù)測(cè)值總和變化不大, 表明相關(guān)醫(yī)療保險(xiǎn)產(chǎn)品的定價(jià)受市場(chǎng)原因或者政策法規(guī)的影響較小, 但模型在約束前后得到的總價(jià)格保持不變, 兼顧了保險(xiǎn)公司的利益, 因而本文的模型和方法能滿足預(yù)期目標(biāo), 有一定的實(shí)用價(jià)值.

表1 有約束條件和無約束條件下λ的預(yù)測(cè)值比較