胡良平

(1.軍事科學(xué)院研究生院，北京 100850； 2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專(zhuān)業(yè)委員會(huì)，北京 100029

1 零膨脹計(jì)數(shù)資料及其零膨脹Poisson分布回歸模型構(gòu)建原理

1.1 零膨脹計(jì)數(shù)資料Poisson分布回歸模型的形式

1.1.1適于零膨脹計(jì)數(shù)資料Poisson分布回歸模型的數(shù)據(jù)結(jié)構(gòu)

適于零膨脹計(jì)數(shù)資料Poisson分布模型的數(shù)據(jù)結(jié)構(gòu)見(jiàn)表1。

【對(duì)數(shù)據(jù)結(jié)構(gòu)的分析】因變量count為“計(jì)數(shù)變量”，它是“魚(yú)的條數(shù)”。不難發(fā)現(xiàn)，沒(méi)有釣到魚(yú)的人數(shù)很多，或者說(shuō)，count=0出現(xiàn)的頻數(shù)很多，這種現(xiàn)象被稱(chēng)為“零膨脹”或“零堆積”。顯然，釣魚(yú)者的“性別”“年齡”可被視為兩個(gè)自變量。

1.1.2零膨脹計(jì)數(shù)資料Poisson分布回歸模型的表達(dá)式

在零膨脹Poisson(簡(jiǎn)稱(chēng)ZIP)分布回歸模型中，早期定義的數(shù)據(jù)產(chǎn)生過(guò)程被稱(chēng)為“過(guò)程2”，即Poisson分布概率函數(shù)，由下面的式(1)給出：

(1)

(2)

表1 研究者隨機(jī)調(diào)查在某公園釣魚(yú)的52人的性別、年齡和過(guò)去半年內(nèi)釣魚(yú)條數(shù)

注：sex，性別；age，年齡(歲)；count，魚(yú)的條數(shù)；F，女；M，男；此資料取自SAS 9.3的“幫助信息”，具體來(lái)說(shuō)，取自SAS/STAT的“FMM過(guò)程”中的“getting started”中的第2個(gè)例子

yi的條件期望(或均值)與條件方差分別見(jiàn)下面的式(3)與式(4)：

E(yi|xi,zi)=μi(1-Fi)

(3)

V(yi|xi,zi)=E(yi|xi,zi)(1+μiFi)

(4)

值得注意的是，比較式(4)與式(3)可知，方差大于均值，故ZIP回歸模型(也包括ZINB回歸模型，即零膨脹負(fù)二項(xiàng)分布回歸模型)刻劃的是具有“零堆積”且具有“過(guò)離散”現(xiàn)象的計(jì)數(shù)因變量隨自變量變化而變化的依賴(lài)關(guān)系。

1.2 零膨脹計(jì)數(shù)資料Poisson分布回歸模型的求解

1.2.1 求解參數(shù)的思路

如何求解出式(2)中的概率“Fi”和回歸系數(shù)“βi”呢？與大多數(shù)“廣義線性模型”[1-2]求解思路大同小異，其常規(guī)做法如下：第1步，基于“最大似然原理”并利用“概率函數(shù)(對(duì)離散型隨機(jī)變量而言)或概率密度函數(shù)(對(duì)連續(xù)型隨機(jī)變量而言)”構(gòu)建“目標(biāo)函數(shù)”，即所謂的“似然函數(shù)”；第2步，求取似然函數(shù)的自然對(duì)數(shù)L；第3步，求取L關(guān)于參數(shù)的一階和二階偏導(dǎo)數(shù)；第4步，令二階偏導(dǎo)數(shù)為“0”，得到方程組；第5步，采用某種迭代計(jì)算方法，求出方程組的解，即求得參數(shù)的估計(jì)值。

1.2.2零膨脹Poisson分布回歸模型的自然對(duì)數(shù)似然函數(shù)

概括地說(shuō)，ZIP回歸模型的自然對(duì)數(shù)似然函數(shù)具有下面的形式，見(jiàn)式(5)：

(5)

在一個(gè)被指定的連接函數(shù)(即Logistic分布函數(shù)或標(biāo)準(zhǔn)正態(tài)分布函數(shù))用于概率Fi之后，就可以寫(xiě)出自然對(duì)數(shù)似然函數(shù)與其導(dǎo)數(shù)的精確表達(dá)式。

1.2.3用Logistic分布函數(shù)作為連接函數(shù)的ZIP回歸模型及參數(shù)估計(jì)

首先，考慮在ZIP回歸模型中用Logistic分布函數(shù)作為連接函數(shù)來(lái)表達(dá)概率Fi，這個(gè)函數(shù)就是下面的式(6)：

(6)

此時(shí)，ZIP回歸模型的自然對(duì)數(shù)似然函數(shù)見(jiàn)下面的式(7)：

(7)

式(7)中的“wi”為“權(quán)重系數(shù)”，視具體情形而定，大體可分為以下6種情形：

(1)wi=1

當(dāng)數(shù)據(jù)集中沒(méi)有可作為“權(quán)重”變量或“頻數(shù)”變量及其取值時(shí)，就令“wi=1”。

(2)wi=Wi

當(dāng)用戶(hù)在“WEIGHT語(yǔ)句”中指定了“非正態(tài)化”選項(xiàng)，而且，在“WEIGHT語(yǔ)句”中指定了變量及其非正態(tài)化取值時(shí)，就令“wi=Wi”。

(4)wi=Fi

在“FREQ語(yǔ)句”中，指定了變量及其取值Fi時(shí)，就令“wi=Fi”。

(5)wi=WiFi

當(dāng)用戶(hù)在“WEIGHT語(yǔ)句”中沒(méi)有使用“非正態(tài)化”選項(xiàng)，而且，同時(shí)使用了“WEIGHT語(yǔ)句”和“FREQ語(yǔ)句”，當(dāng)然，Wi和Fi分別為前面提及的兩個(gè)語(yǔ)句中變量的取值，就令“wi=WiFi”。

對(duì)于ZIP回歸模型的偏導(dǎo)數(shù)分別見(jiàn)下面的式(8)和式(9)：

(8)

(9)

1.2.4用標(biāo)準(zhǔn)正態(tài)分布函數(shù)作為連接函數(shù)的ZIP回歸模型及參數(shù)估計(jì)

首先，考慮在ZIP回歸模型中用標(biāo)準(zhǔn)正態(tài)分布函數(shù)作為連接函數(shù)來(lái)表達(dá)概率，這個(gè)函數(shù)就是下面的式(10)：

(10)

此時(shí)，ZIP回歸模型的自然對(duì)數(shù)似然函數(shù)見(jiàn)下面的式(11)：

(11)

式(11)中對(duì)“wi的定義同式(7)，此處從略。

對(duì)ZIP回歸模型的偏導(dǎo)數(shù)分別見(jiàn)下面的式(12)和式(13)：

(12)

(13)

【說(shuō)明】以上計(jì)算公式來(lái)源于SAS/ETS模塊中“COUNTREG過(guò)程”的“details”部分。

1.2.5 二階偏導(dǎo)數(shù)及迭代計(jì)算

通常情況下，需要在一階偏導(dǎo)數(shù)的基礎(chǔ)上，求取二階偏導(dǎo)數(shù)；進(jìn)而，令各參數(shù)的二階偏導(dǎo)數(shù)及二階混合偏導(dǎo)數(shù)為“0”，形成方程組。再利用某種迭代計(jì)算方法，從而求出各參數(shù)的估計(jì)值。因篇幅所限，詳細(xì)做法此處從略。

2 零膨脹Poisson分布回歸模型的SAS實(shí)現(xiàn)

2.1 創(chuàng)建SAS數(shù)據(jù)集

利用下面的SAS數(shù)據(jù)步程序，創(chuàng)建名為“catch”的臨時(shí)SAS數(shù)據(jù)集：

data catch;

input gender $ age count @@;

if gender=' F ' then sex=0;

else if gender=' M ' then sex=1;

datalines;

2.2 編制出因變量count的頻數(shù)分布表并繪制出其頻數(shù)分布直條圖

利用下面的一個(gè)SAS過(guò)程步程序，編制出因變量count的頻數(shù)分布表，并繪制出其頻數(shù)分布直條圖(因?yàn)閏ount是離散型隨機(jī)變量，只取“0”和“非0正整數(shù)”值)：

proc freq data=catch;

tables count/out=aaa plots=freqplot;

run;

【SAS輸出結(jié)果】

count頻數(shù)百分比累積頻數(shù)累積百分比02853.852853.851611.543465.38223.853669.23335.773975.00411.924076.92535.774382.69611.924484.62711.924586.54811.924688.461011.924790.381111.924892.311223.855096.151823.8552100.00

以上是計(jì)數(shù)因變量count的頻數(shù)分布表，count=0出現(xiàn)的頻數(shù)為28，占總頻數(shù)52的53.85%。

以下是“計(jì)數(shù)因變量count的頻數(shù)分布直條圖”，見(jiàn)圖1。

圖1 計(jì)數(shù)因變量count的頻數(shù)分布直條圖

圖1直觀、形象地呈現(xiàn)出計(jì)數(shù)因變量count取“0”的頻數(shù)非常多，這種現(xiàn)象被稱(chēng)為“零膨脹”或“零堆積”。

2.3 求出因變量count的均值和方差

利用下面的兩個(gè)SAS過(guò)程步程序，求出除“0”之外的其他計(jì)數(shù)資料的方差和均值：

proc univariate data=catch(where=(count ne 0)) noprint;

var count;

output out=aaa mean=mfish var=vfish;

run;

proc print data=aaa;

run;

【SAS輸出結(jié)果】

Obsmfishvfish15.8333327.0145

求得除“0”之外的其他計(jì)數(shù)因變量的均值為5.83333、方差為27.0145，方差約為均值的4.63倍。

Obsmfishvfish12.6923120.8054

以上是基于計(jì)數(shù)因變量count的全部52個(gè)觀測(cè)值計(jì)算所得到的均值和方差，方差約為均值的7.73倍。

由此可知，此計(jì)數(shù)資料屬于“過(guò)離散計(jì)數(shù)資料”。

2.4 基于全部自變量對(duì)因變量count構(gòu)建多重零膨脹Poisson分布回歸模型

2.4.1利用下面的SAS過(guò)程步程序嘗試性構(gòu)建多重Poisson分布回歸模型

/*僅擬合多重Poisson分布回歸模型*/

proc fmm data=catch;

class gender;

model count=gender*age/dist=Poisson;

run;

【SAS輸出結(jié)果】

Fit Statistics-2 Log Likelihood 182.7AIC (smaller is better)188.7AICC (smaller is better)189.2BIC (smaller is better)194.6 Pearson Statistic85.9573

關(guān)于模型對(duì)資料的擬合效果，暫不評(píng)價(jià)，需要與其他同類(lèi)模型的擬合結(jié)果進(jìn)行比較，才能得出有說(shuō)服力的判定。

2.4.2利用下面的SAS過(guò)程步程序嘗試性構(gòu)建多重負(fù)二項(xiàng)分布回歸模型

/*僅擬合多重負(fù)二項(xiàng)分布回歸模型*/

proc fmm data=catch;

class gender;

model count=gender*age/dist=nb;

run;

【SAS輸出結(jié)果】

Fit Statistics-2 Log Likelihood162.1AIC (smaller is better)170.1AICC (smaller is better)171.0BIC (smaller is better)177.9 Pearson Statistic38.3597

Parameter Estimates for 'Negative Binomial' Model效應(yīng)gender估計(jì)值標(biāo)準(zhǔn)誤差z值Pr>|z|Intercept-4.38160.8160-5.37<0.0001age*genderF0.13780.019577.04<0.0001age*genderM0.11300.019715.73<0.0001Scale Parameter0.66880.3263

比較以上兩個(gè)統(tǒng)計(jì)模型給出的“擬合統(tǒng)計(jì)量”部分結(jié)果可知，“多重負(fù)二項(xiàng)分布回歸模型”好于“多重Poisson分布回歸模型”，因?yàn)锳IC、AICC、BIC和“Pearson Statistic(皮爾遜卡方統(tǒng)計(jì)量)”數(shù)值均有所下降。

2.4.3利用下面的SAS過(guò)程步程序構(gòu)建多重零膨脹Poisson分布回歸模型

(1)按性別分層且基于常數(shù)與Poisson分布分別擬合零膨脹與非零膨脹兩部分的零膨脹Poisson分布回歸模型

proc fmm data=catch;

class gender;

model count = gender*age / dist=Poisson;

model + / dist=Constant;

run;

【SAS主要輸出結(jié)果】

Fit Statistics-2 Log Likelihood145.6AIC (smaller is better)153.6AICC (smaller is better)154.5BIC (smaller is better)161.4Pearson Statistic43.4467Effective Parameters4Effective Components2

Parameter Estimates for 'Poisson' Model成分效應(yīng)gender估計(jì)值標(biāo)準(zhǔn)誤差z值Pr>|z|1Intercept-3.52150.6448-5.46<0.00011age*genderF0.12160.013449.04<0.00011age*genderM0.10560.013947.58<0.0001Parameter Estimates for Mixing Probabilities效應(yīng)鏈接尺度估計(jì)值標(biāo)準(zhǔn)誤差z值Pr>|z|概率Intercept0.83420.47681.750.08020.6972

關(guān)于模型對(duì)資料的擬合效果，暫不評(píng)價(jià)，需要與其他同類(lèi)模型的擬合結(jié)果進(jìn)行比較，才能得出有理由的判定。若僅依據(jù)“AIC、AICC、BIC”的數(shù)值來(lái)判斷，要好于前面的“多重負(fù)二項(xiàng)分布回歸模型”的擬合效果。

(2)不按性別分層且基于常數(shù)與Poisson分布分別擬合零膨脹與非零膨脹兩部分的零膨脹Poisson分布回歸模型

proc fmm data=catch;

class gender;

model count = gender age / dist=Poisson;

model + / dist=Constant;

run;

【SAS主要輸出結(jié)果】

Fit Statistics-2 Log Likelihood144.4AIC (smaller is better)152.4AICC (smaller is better)153.3BIC (smaller is better)160.2Pearson Statistic43.7210Effective Parameters4Effective Components2

Parameter Estimates for 'Poisson' Model成分效應(yīng)gender估計(jì)值標(biāo)準(zhǔn)誤差z值Pr>|z|1Intercept-4.12610.6624-6.23<0.00011genderF0.79220.21933.610.00031genderM01age0.11780.01338.86<0.0001Parameter Estimates for Mixing Probabilities效應(yīng)鏈接尺度估計(jì)值標(biāo)準(zhǔn)誤差z值Pr>|z|概率Intercept0.91900.49341.860.06250.7148

若僅依據(jù)“AIC、AICC、BIC”的數(shù)值來(lái)判斷，要略好于前面“按性別分層”的擬合效果，但二者無(wú)差異。

2.5 最終的回歸模型與專(zhuān)業(yè)結(jié)論

2.5.1 最終的回歸模型

依據(jù)上面最后一次選定的回歸模型擬合資料所得到的計(jì)算結(jié)果寫(xiě)出“多重零膨脹Poisson分布回歸模型”，此模型由兩個(gè)表達(dá)式組成，一個(gè)對(duì)應(yīng)于“yi=0”；另一個(gè)對(duì)應(yīng)于“yi>0”。模型表達(dá)式見(jiàn)前面的式(2)，其中，F(xiàn)i=0.7148、μi=e-4.1261+0.7922(sex=F)+0.1178age。將這兩個(gè)參數(shù)的估計(jì)結(jié)果代入式(2)，就可獲得與表1資料對(duì)應(yīng)的“零膨脹Poisson分布回歸模型”。

2.5.2 專(zhuān)業(yè)結(jié)論

任何一人釣不到魚(yú)的概率約為71.5%(因上面第一式“加號(hào)”后面與釣魚(yú)者性別、年齡有關(guān)，但其數(shù)值并不大)；女性較男性、年齡長(zhǎng)者較年齡小者釣魚(yú)數(shù)目多。