0-k 膨脹計數回歸模型的參數估計與統計檢驗

李春玉,安博文,李曉天

（1.河北經貿大學數學與統計學學院,河北 石家莊 050061；2.華僑大學經濟與金融學院,福建 泉州 362021；3.新疆理工學院經濟貿易與管理學院,新疆 阿克蘇 843000）

觀測數值為自然數的數據就是計數數據,其在生物醫學、臨床試驗、保險索賠、風險管控以及社會科學等方面都得到了廣泛應用,因此計數數據的回歸預測、統計檢驗等問題一直是數據分析的熱點問題.在實際應用中,0 膨脹現象、0-1 膨脹現象或0-k膨脹現象都較為常見.就0-k膨脹而言,當k取值為0時即為0 膨脹分布,當k取值為1 時即為0-1 膨脹分布.無論是0 膨脹分布還是0-1 膨脹分布均是0-k膨脹分布的特例,因此本文將0-k膨脹計數回歸模型作為研究對象.

0-k膨脹現象是指0 和k同時膨脹,樣本數據表現為觀測值中的0 值和k值頻率都較高.0-k膨脹需要滿足0 和k同時膨脹,這就要求在回歸模型中需要包含兩個膨脹參數,一個膨脹參數控制0 膨脹,另一個膨脹參數控制k膨脹.從現有文獻來看,部分學者對0 膨脹、0-1 膨脹和0-k膨脹的回歸模型都進行了深入研究,例如：Lambert[1]首次基于泊松分布提出0 膨脹回歸模型,同時引入協變量對膨脹數據進行解釋；Melkersson 和Olssen[2]提出0-1 膨脹泊松回歸模型,并結合看牙醫次數的統計數據進行實例分析；Lin 和Tsai[3]針對健康調查數據建立0-k膨脹泊松回歸模型進行統計分析.上述文獻對膨脹模型的建立與參數估計都進行了詳細研究,但并未涉及膨脹類型數據的統計檢驗.如果調查數據僅服從0 膨脹分布或者僅服從k膨脹分布,而在建模過程中依舊采用0-k膨脹分布進行分析,這可能會由于因變量分布設定錯誤而導致回歸結果出現偏差.基于此,本文期望通過統計檢驗對0-k膨脹進行識別以便精準建模.

回顧之前的相關研究,關于回歸系數的估計方法的研究結果如下：Fahrmeir 和Echavarria[4]研究了一類可加的0 膨脹模型,Ghosh 等[5]給出了0 膨脹回歸模型的貝葉斯估計法,Zhang 等[6]討論了0-1 膨脹Poisson 分布的統計性質并給出了回歸系數的區間估計,楊亮和孟生旺[7]建立0 膨脹分位回歸模型并給出了相關參數的貝葉斯估計結果,Tang 等[8]將0-1 膨脹模型重參數化并給出了EM 算法和貝葉斯算法兩種估計方法,夏麗麗和田茂再[9]在0-1 膨脹泊松回歸模型的泊松部分建立非參數模型并采用EM算法給出了估計結果,安博文等[10]基于Fisher 信息矩陣給出了0-1 膨脹貝塔- 二項分布參數的置信區間,李春玉等[11]采用EM算法給出了0-1 膨脹貝塔- 負二項分布參數的置信區間.就膨脹參數的統計檢驗研究結果而言,Broek[12]采用Score統計量對0 膨脹泊松分布的參數進行檢驗,Moghimbeigi 等[13]通過抽樣分布測試了Score統計量的檢驗功效,Xiang 和Teo[14]采用Wald統計量和LR統計量對0 膨脹泊松分布進行了檢驗,Hsu 等[15]建立混合權重0 膨脹回歸模型并采用Wald統計量進行檢驗,黃倩[16]采用Score統計量對0-k膨脹泊松回歸模型進行檢驗.劉娛等[17]分別采用LR統計量、Wald統計量和Score統計量對0-1膨脹泊松回歸模型進行了檢驗.

綜上所述,本研究將對0-k膨脹回歸模型進行參數估計和統計檢驗,借鑒Tang 等[8]對0-1 膨脹泊松回歸模型的設定方法引入隱變量建立0-k膨脹回歸模型并基于EM 算法給出回歸系數與膨脹參數的估計結果；基于Fisher 信息矩陣構造各個參數的置信區間；采用Wald統計量和LR統計量對0-k膨脹分布是否會出現退化進行檢驗；以0-k膨脹泊松回歸模型為例結合相關數據進行實例分析.

1 模型設定

假設隨機變量V服從非膨脹的基礎分布,記為V~f(θ|v)（其中θ 表示分布f(θ|v)中原有的參數,和分布膨脹與否無關）；引入隱變量B和隱變量C分別控制0 膨脹部分和k膨脹部分,且B~B(1,p)、C~B(1 ,q)；要求隨機變量V、B和C相互獨立.令Y=V(1?B)+k×B(1?C),則隨機變量Y所服從的分布即為0-k膨脹分布（簡記為ZKI分布）.Y與(B,C,V)的對應關系如下

研究樣本Y1,???,Yn獨立同分布于0-k膨脹分布,用y=(y1, ???,yn)表示樣本的觀測數據,則0-k膨脹分布的概率質量函數如式（1）.

式（1）中：θi表示非膨脹分布的原有參數,pi為0 膨脹參數,qi為k膨脹參數.經過計算,式（1）的數學期望表示為.當pi＞0 且qi=0 時,式（1）退化為k膨脹分布（簡記為KI分布）；當pi＞0 且qi=1 時,式（1）退化為0 膨脹分布（簡記為ZI分布）；當pi=0 時,式（1）退化為非膨脹的基礎分布（簡記為NI分布）.

將參數θi、pi和qi作為被解釋變量,引入協變量zi、xi和wi,假定鏈接函數為

式（2）中：zi、xi和wi為可觀測數據, β、γ 和α 為待估系數,觀測數據與待估系數的具體形式如下

2 參數估計

這部分基于EM算法和Fisher 信息矩陣給出0-k膨脹分布回歸系數的點估計值與區間估計值.

2.1 點估計

被解釋變量y以及協變量z、x和w均可以被觀測,隱變量B和隱變量C不可被觀測,這里假定其樣本數據形式依次為b=(b1, ???,bn)和c=(c1, ???,cn).由此可以得到擴充數據的對數似然函數為

將式（2）代入式（3）中即可得到含有協變量的對數似然函數.

以下采用EM算法求解式（4）的最大化問題.EM 算法分為E 步和M 步,E 步為計算隱變量B和隱變量C的數學期望,M步為最大化對數似然函數式（4）.

E 步：根據文獻[8]可以分別計算出隱變量B和隱變量C的數學期望,如式（5）

M步：采用極大似然估計法求解式（4）中的參數.具體過程如式（7）

這里采用Newton-Raphso 迭代計算方程組（7）的數值解,記,則β、γ 和α 的極大似然估計結果為

式（8）中：t表示Newton-Raphso 算法的迭代次數.

2.2 區間估計

用空間A表示樣本y的所有觀測值,空間A0表示樣本y中所有值為0 的觀測值,空間Ak表示樣本y中所有值為k的觀測值.記A中元素個數為n,A0中元素個數為n0,Ak中元素個數為nk.對式（3）和式（4）進行改寫,寫出不含隱變量的對數似然函數,具體見式（9）

根據式（10）可以計算回歸系數β、γ 和α 的Fisher 信息矩陣,記為J.由于回歸模型的參數個數為(k+r+m+3),則J為(k+r+m+3)階方陣.考慮第（ii= 1, ???,k+r+m+3）個回歸系數o（ioi∈{0, ???,},0, ???,r,0, ???,m}）,其對應的Fisher 信息矩陣為,其中

由此可得回歸系數oi的方差為

從而,回歸系數oi的(1?)置信區間為

3 膨脹參數檢驗

采用Wald統計量和LR統計量對0-k膨脹分布的膨脹參數進行假設檢驗.0-k膨脹是指0 和k同時膨脹,因此需要考慮3 種情況：①是否退化為0 膨脹分布；②是否退化為k膨脹分布；③是否退化為非膨脹的基礎分布.對此提出3 個假設檢驗：①是否退化為0 膨脹分布的檢驗,原假設表述為0-k膨脹分布會退化為0 膨脹分布,則H0為p＞0 且q=1,拒絕H0說明該分布不會退化為0 膨脹分布；②是否退化為k膨脹分布的檢驗,原假設表述為0-k膨脹分布會退化為k膨脹分布,則H0為p＜0 且q=0,拒絕H0說明該分布不會退化為k膨脹分布；③是否退化為非膨脹分布的檢驗,原假設表述為0-k膨脹分布會退化為非膨脹分布,則H0為p=0,拒絕H0說明該分布不會退化為非膨脹分布.

根據式（9）可以計算關于膨脹參數的Fisher 信息矩陣.記

其中

Wald檢驗統計量和LR檢驗統計量分別為：

是否退化為0 膨脹分布的檢驗統計量

是否退化為k膨脹分布的檢驗統計量

是否退化為非膨脹分布的檢驗統計量

4 實例應用

這部分以0-k膨脹泊松回歸模型為例進行實際應用.取基礎分布（其中=）,取鏈接函數(·)=log (·)、(·)=logit (·)和(·)=logit(·),則0-k膨脹泊松回歸模型的具體形式為

其中

4.1 回歸系數估計

實例數據來源于pscl 包中提供的bioChemists 數據集,bioChemists 數據集以915 名生物化學博士為調查對象,記錄了他們在博士后三年期間發表論文數量（ART）、博士后性別（FEM）、是否婚配（MAR）、所擁有5 歲及以下的子女數量（KID）、博士在讀學校的聲望（PHD）以及他們導師近三年撰寫文章的數量（MENT）.經過篩選,本文選取ART作為被解釋變量,將MENT和KID分別作為泊松部分、0 膨脹部分和k膨脹部分的協變量,具體模型如下

于ART變量而言,觀測樣本量為915,其中,觀測數值為0 的樣本量是275,占比為30.05%；觀測數值為1 的樣本量是246,占比為26.89%；觀測數值為2 的樣本量是178,占比為19.45%.因此這里依次將k取值為1 和2,分別建立0-1 膨脹泊松回歸模型和0-2 膨脹泊松回歸模型,具體估計結果如表1、表2所示.

表1 0-1 膨脹泊松回歸模型系數估計結果Tab.1 0-1 Coefficient estimation results of expansion Poisson regression model

表2 0-2 膨脹泊松回歸模型系數估計結果Tab.2 0-2 Coefficient estimation results of expansion Poisson regression model

回歸系數的估計結果顯示,于0-1 膨脹泊松回歸模型而言,泊松部分的ART與MENT呈正向變動關系,即博士后導師發表論文數量越多則博士后發表論文數量也越多,說明博士后導師發表論文會對博士后發表論文產生積極影響；ART與KID呈反向變動關系,5 歲及以下子女數量越多則博士后發表論文數量越少,說明隨著5 歲及以下子女數量增加博士后的時間精力被逐漸分散,從而發表論文數量減少.0 膨脹部分的ART與MENT呈反向變動關系,隨著博士后導師發表論文數量增加博士后不發表論文的可能性逐漸減??；ART與KID呈正向變動關系,隨著5 歲及以下子女數量增加博士后不發表論文的可能性逐漸增大.1 膨脹部分的ART與MENT呈反向變動關系,博士后導師發表論文數量增加會減小博士后發表1 篇論文的可能性；ART與KID呈正向變動關系,5 歲及以下子女數量增加會增大博士后發表1篇論文的可能性.于0-2 膨脹泊松回歸模型而言,泊松部分、0 膨脹部分都與0-1 膨脹泊松回歸模型的估計結果相近,回歸系數的正負號和數量級均未發生明顯變化,通過對比兩種回歸模型也驗證了泊松部分、0 膨脹部分估計結果的穩健性.2 膨脹部分,ART與MENT呈反向變動關系,1 膨脹部分MENT的系數為-0.073 8,2 膨脹部分MENT的系數為-0.116 0,隨著博士后導師發表論文數量增加博士后發表1篇論文可能性的減弱程度要低于發表兩篇論文可能性的減弱程度；ART與KID呈正向變動關系,1 膨脹部分KID的系數為0.211 8,2 膨脹部分KID的系數為0.389 3,隨著5 歲及以下子女數量增加博士后發表2 篇論文可能性的加強程度要高于發表1 篇論文可能性的加強程度.

4.2 分布參數估計與膨脹參數檢驗

這部分對0-k膨脹分布是否會出現退化的現象進行假設檢驗.在假設檢驗之前需要計算相關的膨脹參數,表3 和表4 依次給出了0-1 膨脹泊松分布參數的估計結果和0-2 膨脹泊松分布參數的估計結果.于0-1 膨脹泊松分布而言,0 膨脹參數的估計值為0.417 2,1 膨脹參數的估計值為0.619 9,泊松部分的參數估計值為2.632 8,經過計算得到該分布的數學期望為1.693 0；于0-2 膨脹泊松分布而言,0 膨脹參數的估計值為0.336 9,2 膨脹參數的估計值為0.714 8,泊松部分的參數估計值為2.408 2,經過計算得到該分布的數學期望為1.789 0.從樣本數據來看,ART的平均值為1.692 9（或稱為真實值）,基于0-1 膨脹泊松分布計算的數學期望與真實值十分接近,相對誤差在±0.01%之內；基于0-2 膨脹泊松分布計算的數學期望與真實值的相對誤差為5.68%,說明0-2 膨脹泊松分布會高估原始數據的均值.由此可見,采用0-1 膨脹泊松回歸模型進行建模分析會更加可靠.

表3 0-1 膨脹泊松分布的參數估計結果Tab.3 0-1 Parameter estimation results of expansion Poisson distribution

表4 0-2 膨脹泊松分布的參數估計結果Tab.4 0-2 Parameter estimation results of expansion Poisson distribution

從統計檢驗的角度而言,表5 和表6 分別給出了0-1 膨脹泊松分布是否會退化和0-2 膨脹泊松分布是否會退化的檢驗結果.通過查表可得自由度為2、置信水平為0.05 的卡方分布臨界值是5.99.從表5來看,Wald統計量與LR統計量均拒絕了退化為0 膨脹泊松分布、退化為1 膨脹泊松分布以及退化為非膨脹的泊松分布的原假設,說明原始數據存在0-1 膨脹現象,可以采用0-1 膨脹泊松回歸模型進行建模分析.從表6 來看,Wald統計量均拒絕了退化為0 膨脹泊松分布、退化為1 膨脹泊松分布以及退化為非膨脹的泊松分布的原假設,而LR統計量僅拒絕了退化為2 膨脹泊松分布和退化為非膨脹的泊松分布的原假設,但并未拒絕退化為0 膨脹泊松分布的原假設,說明從LR檢驗來看0-2 膨脹泊松分布會退化為0 膨脹泊松分布,這意味著原始數據不存在2 膨脹的現象,因此不應該采用0-2 膨脹泊松回歸模型進行建模分析.綜上所述,對于樣本數據更適合建立0-1 膨脹泊松回歸模型而非0-2 膨脹泊松回歸模型,這一結論與前文所述的“0-2 膨脹泊松分布會高估原始數據均值”相互佐證.

表5 0-1 膨脹泊松分布的膨脹參數檢驗結果Tab.5 0-1 Expansion parameter test results of Poisson distribution of expansion

表6 0-2 膨脹泊松分布的膨脹參數檢驗結果Tab.6 0-1 Expansion parameter test results of Poisson distribution of expansion

5 結語

本研究通過引入隱變量設定0-k膨脹回歸模型,采用EM算法給出了各個參數的點估計值并基于Fisher 信息矩陣構造了置信區間；對0-k膨脹分布是否會退化的現象進行假設檢驗,通過Wald統計量和LR統計量給出了檢驗方法；以0-k膨脹泊松回歸模型為例進行實際應用,包括回歸系數的估計、分布參數的估計以及膨脹參數的統計檢驗.由此得到以下主要結論：采用EM 算法對參數進行點估計時,回歸系數與膨脹參數估計值的迭代收斂速度較快并且估計結果穩??；基于Fisher 信息矩陣給出了0-k膨脹回歸模型的區間估計結果,為參數的置信區間構造提供了一種思路；于0-k 膨脹泊松分布而言,當k取不同值時可能會對原始數據的擬合出現偏差,在統計檢驗時發現LR 統計量的檢驗效果要優于Wald統計量的檢驗效果.