胡良平

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

1 過(guò)離散計(jì)數(shù)資料及其負(fù)二項(xiàng)分布模型構(gòu)建原理

1.1 適于過(guò)離散計(jì)數(shù)資料負(fù)二項(xiàng)分布回歸模型的數(shù)據(jù)結(jié)構(gòu)

適于過(guò)離散計(jì)數(shù)資料負(fù)二項(xiàng)分布模型的數(shù)據(jù)結(jié)構(gòu)見(jiàn)表1[1]。

表1數(shù)據(jù)來(lái)自德國(guó)的衛(wèi)生改革項(xiàng)目。改革的目的是減少患者到醫(yī)生處的就診次數(shù)。數(shù)據(jù)共包括2 227例觀(guān)測(cè)，分別屬于改革之前的1996年和改革之后的1998年。響應(yīng)變量為患者在3個(gè)月之內(nèi)的就診次數(shù)，自變量為改革前后(0=改革前，1=改革后)、健康狀況(0=良好，1=不良)、年齡、受教育時(shí)間和家庭收入的對(duì)數(shù)值。由于數(shù)據(jù)量較大，表1中僅列出部分結(jié)果。求出“就診次數(shù)”的均值與方差如下。

均值方差12.5891316.1299

【對(duì)數(shù)據(jù)結(jié)構(gòu)的分析】因變量Y“就診次數(shù)”為“計(jì)數(shù)變量”，其均值為2.589、方差為16.130，方差明顯大于均值，只能將Y視為服從負(fù)二項(xiàng)分布的離散型隨機(jī)變量；擬考察的四個(gè)自變量中，有兩個(gè)是“二值變量”，另外兩個(gè)是“計(jì)量變量”。雖然，在統(tǒng)計(jì)學(xué)的理論上，默認(rèn)自變量全為“計(jì)量變量”，但從實(shí)用性角度出發(fā)，“二值”自變量也是可以接受的。

表1 患者就診次數(shù)及其影響因素?cái)?shù)據(jù)

1.2 過(guò)離散計(jì)數(shù)資料負(fù)二項(xiàng)分布回歸模型的構(gòu)建與求解[2]

1.2.1 Poisson分布回歸模型的推廣

在對(duì)第i個(gè)個(gè)體進(jìn)行觀(guān)測(cè)時(shí)，通過(guò)引入一個(gè)觀(guān)測(cè)不到的非齊性項(xiàng)“τi”，從而使Poisson分布回歸模型得到推廣。因此，假定個(gè)體之間以隨機(jī)的方式呈現(xiàn)彼此的差異，原因在于可觀(guān)測(cè)的協(xié)變量不能完全解釋不同個(gè)體之間實(shí)際存在的差異。此種情況可用如下的式子來(lái)呈現(xiàn)，見(jiàn)式(1)：

(1)

在式(1)中，觀(guān)測(cè)不到的非齊性項(xiàng)τi=eεi獨(dú)立于自變量向量xi。于是，yi關(guān)于xi和τi的條件分布是關(guān)于條件均值和條件方差均為μiτi的Poisson分布，其概率函數(shù)見(jiàn)式(2)：

(2)

令g(τi)為τi的概率密度函數(shù)，于是，通過(guò)對(duì)f(yi|xi,τi)求關(guān)于τi的積分，便可求得f(yi|xi)的概率分布，見(jiàn)式(3)：

(3)

在式(3)中，當(dāng)假定τi服從伽瑪分布時(shí)，其積分的結(jié)果存在解析解，這個(gè)解就是負(fù)二項(xiàng)分布。

為了確定分布的均值，當(dāng)模型中包含常數(shù)項(xiàng)時(shí)，必需假定E(eεi)=E(τi)=1。因此，假定τi為服從gamma(θ,θ)分布且其均值和方差分別為E(τi)=1和V(τi)=1/θ，其概率密度函數(shù)見(jiàn)式(4)：

(4)

在式(4)中，分母中的伽瑪函數(shù)見(jiàn)式(5)：

(5)

而式(4)中的θ是一個(gè)正的參數(shù)。

1.2.2 負(fù)二項(xiàng)分布概率函數(shù)形式之一

于是，在給定了xi的條件下，推導(dǎo)出yi的概率函數(shù)見(jiàn)式(6)：

(6)

1.2.3 負(fù)二項(xiàng)分布概率函數(shù)形式之二

yi=0,1,2…

(7)

因此，負(fù)二項(xiàng)分布概率函數(shù)是由服從Poisson分布與伽瑪分布的兩類(lèi)隨機(jī)變量混合而成的。其條件均值和條件方差分別見(jiàn)式(8)和式(9)：

(8)

(9)

由式(9)可知，負(fù)二項(xiàng)分布的條件方差明顯大于其條件均值，這種“過(guò)離散”的結(jié)果是忽視未觀(guān)測(cè)到的非齊性項(xiàng)所致。

式(9)顯示出：負(fù)二項(xiàng)分布的方差函數(shù)是其均值的二次方，這個(gè)結(jié)果可被稱(chēng)為第2種定義下的負(fù)二項(xiàng)分布回歸模型(即“Cameron and Trivedi 1986”)，簡(jiǎn)稱(chēng)為“2型負(fù)二項(xiàng)分布回歸模型”。當(dāng)α=0時(shí)，負(fù)二項(xiàng)分布就退化成為Poisson分布了。檢驗(yàn)“α=0”這個(gè)假設(shè)是否成立，可使用WALD檢驗(yàn)。

1.2.4 2型負(fù)二項(xiàng)分布回歸模型的自然對(duì)數(shù)似然函數(shù)及一階偏導(dǎo)數(shù)

2型負(fù)二項(xiàng)分布回歸模型的自然對(duì)數(shù)似然函數(shù)可由下面的式(10)給出：

(10)

在式(10)中，wi的定義很繁瑣，參見(jiàn)有關(guān)文獻(xiàn)，此處從略；在式(10)等號(hào)右邊第二個(gè)求和號(hào)之后的第一項(xiàng)是怎么得到的？它是由如下的原因所致：在式(7)中，兩個(gè)伽瑪函數(shù)相除可以改寫(xiě)成乘積形式，見(jiàn)下面的式(11)：

(11)

基于式(10)，求L關(guān)于參數(shù)β的偏導(dǎo)數(shù)，見(jiàn)下面的式(12)：

(12)

基于式(10)，求L關(guān)于參數(shù)α的偏導(dǎo)數(shù)，見(jiàn)下面的式(13)：

(13)

1.2.5 1型負(fù)二項(xiàng)分布回歸模型的自然對(duì)數(shù)似然函數(shù)及一階偏導(dǎo)數(shù)

V(yi|xi)=μi+αμi

(14)

為了估計(jì)這樣的模型，在“model語(yǔ)句”中指定“DIST=NEGBIN(p=1)”這樣的選項(xiàng)。

于是，1型負(fù)二項(xiàng)分布回歸模型的自然對(duì)數(shù)似然函數(shù)由下面的式(15)給出：

(15)

在式(15)中，wi為權(quán)重系數(shù)，其定義視具體情況而定，因篇幅所限，此處從略。

基于式(15)，求L關(guān)于參數(shù)β的偏導(dǎo)數(shù)，見(jiàn)下面的式(16)：

α-1ln(1+α)μixi}

(16)

基于式(15)，求L關(guān)于參數(shù)β的偏導(dǎo)數(shù)，見(jiàn)下面的式(17)：

(17)

1.2.6 二階偏導(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 過(guò)離散計(jì)數(shù)資料負(fù)二項(xiàng)分布回歸模型的SAS實(shí)現(xiàn)[3]

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

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

data nbd1;

infile 'd:SASTJFXmdvisit.dat';

input id numvisit reform badh age edu loginc;

run;

【SAS程序說(shuō)明】將本例中全部2 227行數(shù)據(jù)(注意：表1中僅列出了前9行和最后一行)錄入計(jì)算機(jī)，以“文本格式”且以“mdvisit.dat”為文件名，存儲(chǔ)在名為“SASTJFX”的D盤(pán)文件夾內(nèi)。每行有7個(gè)數(shù)據(jù)，分別為“編號(hào)(id)”“三個(gè)月內(nèi)的就診次數(shù)(numvisit)”“改革前后(reform)(0=改革前，1=改革后)”“健康狀況(badh)(0=良好，1=不良)”“年齡(age)”“受教育時(shí)間(edu)”和“家庭收入的對(duì)數(shù)值(loginc)”。

2.2 求出因變量numvisit的均值和方差

利用下面的兩個(gè)SAS過(guò)程步程序，求出因變量Y的均值和方差：

proc univariate data=nbd1 noprint;

var numvisit;

output out=aaa mean=ybar var=yvar;

run;

proc print data=aaa;

var ybar yvar;

run;

【SAS輸出結(jié)果】

Obsybaryvar12.5891316.1299

求出3個(gè)月內(nèi)就診次數(shù)(numvisit)的均值和方差分別為2.589和16.130。此結(jié)果表明，方差約為均值的6.23倍，屬于較嚴(yán)重的“過(guò)離散”。

2.3 基于全部自變量對(duì)因變量Y構(gòu)建多重Poisson分布回歸模型

利用下面的SAS過(guò)程步，嘗試構(gòu)建Poisson分布回歸模型，目的是為了使用選項(xiàng)“noscale”進(jìn)行拉格朗日乘子檢驗(yàn)，以明確是否存在明顯的“過(guò)離散”現(xiàn)象。

proc genmod data=nbd1;

model numvisit= reform badh age edu loginc/link=log dist=nb noscale;

run;

【SAS程序說(shuō)明】“l(fā)ink=log”表明采用“自然對(duì)數(shù)”作為聯(lián)接函數(shù)，即對(duì)計(jì)數(shù)因變量取自然對(duì)數(shù)變換；“dist=nb”表明要基于負(fù)二項(xiàng)分布構(gòu)建負(fù)二項(xiàng)分布回歸模型。但是，選項(xiàng)“noscale”的作用是將冗余參數(shù)k的取值固定為“0”，此時(shí)，相當(dāng)于構(gòu)建Poisson回歸模型。與此同時(shí)，該選項(xiàng)的主要用途是進(jìn)行拉格朗日乘子檢驗(yàn)，以明確因變量是否存在明顯的“過(guò)離散”現(xiàn)象。

【SAS主要輸出結(jié)果】基于Poisson分布回歸模型的參數(shù)估計(jì)結(jié)果沒(méi)有實(shí)用價(jià)值，此處從略。

拉格朗日乘數(shù)統(tǒng)計(jì)量參數(shù)卡方Pr > 卡方離散度551.7077<.0001*

注：*單側(cè)P值

以上結(jié)果表明，因變量存在極嚴(yán)重的“過(guò)離散”現(xiàn)象。

2.4 嘗試對(duì)“過(guò)離散”現(xiàn)象進(jìn)行校正后再構(gòu)建Poisson分布回歸模型

利用下面的SAS過(guò)程步程序?qū)Α斑^(guò)離散”現(xiàn)象進(jìn)行校正后再構(gòu)建Poisson分布回歸模型：

proc genmod data=nbd1;

model numvisit= reform badh age edu loginc/link=log dist=poisson scale=deviance;

run;

【SAS程序說(shuō)明】選項(xiàng)“scale=deviance”是要求對(duì)“過(guò)離散”進(jìn)行校正后，再構(gòu)建Poisson分布回歸模型。

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

評(píng)估擬合優(yōu)度的準(zhǔn)則準(zhǔn)則自由度值值/自由度偏差22217422.12443.3418調(diào)整后的偏差22212221.00001.0000Pearson 卡方22219681.69204.3592調(diào)整后的 Pearson X222212897.15411.3044對(duì)數(shù)似然129.4790完全對(duì)數(shù)似然-5943.8280AIC(越小越好)11899.6561AICC(越小越好)11899.6939BIC(越小越好)11933.9066

以上是采用Poisson分布回歸模型擬合此資料所對(duì)應(yīng)的“擬合優(yōu)度”評(píng)價(jià)指標(biāo)及其取值，這些結(jié)果需要與其他同類(lèi)模型比較，才有價(jià)值。

最大似然參數(shù)估計(jì)值的分析參數(shù)自由度估計(jì)值標(biāo)準(zhǔn)誤差Wald 95%置信限Wald 卡方Pr > 卡方Intercept1-0.42150.4917-1.38520.54220.730.3913reform1-0.13990.0485-0.2350-0.04488.310.0039badh11.13260.05541.02411.2412418.17<.0001age10.00490.00230.00040.00944.550.0329edu1-0.01180.0109-0.03320.00951.180.2781loginc10.15200.06580.02310.28105.340.0208尺度01.82810.00001.82811.8281

Note: The scale parameter was estimated by the square root of DEVIANCE/DOF

最后一行“注釋”的含義是：尺度參數(shù)是“離差/自由度”的平方根，即(7422.124/2221)的平方根，其數(shù)值為1.8281。

以上結(jié)果表明：尺度參數(shù)由原先的“1”調(diào)整為“1.8281”，由此，模型中各參數(shù)的“標(biāo)準(zhǔn)誤”都得到“校正”。從而，最后兩列的數(shù)值也都得到校正。

2.5 在“過(guò)離散”現(xiàn)象嚴(yán)重的情況下構(gòu)建多重負(fù)二項(xiàng)分布回歸模型

proc genmod data=nbd1;

model numvisit= reform badh age edu loginc/link=log dist=nb;

run;

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

評(píng)估擬合優(yōu)度的準(zhǔn)則準(zhǔn)則自由度值值/自由度偏差22212412.34641.0862調(diào)整后的偏差22212412.34641.0862Pearson 卡方22212693.55141.2128調(diào)整后的 Pearson X222212693.55141.2128對(duì)數(shù)似然1813.1299完全對(duì)數(shù)似然-4563.3905AIC(越小越好)9140.7809AICC(越小越好)9140.8314BIC(越小越好)9180.7398

以上“擬合優(yōu)度評(píng)價(jià)指標(biāo)的取值”與前面基于“Poisson分布回歸模型”的相應(yīng)評(píng)價(jià)指標(biāo)的取值相比，現(xiàn)在的模型，即“負(fù)二項(xiàng)分布回歸模型”對(duì)此資料的擬合效果好多了(偏差變小了，AIC、AICC和BIC都變小了很多)。

以上結(jié)果表明，截距項(xiàng)無(wú)統(tǒng)計(jì)學(xué)意義，去掉截距項(xiàng)后，可使模型精簡(jiǎn)一些。

2.6 構(gòu)建不含截距項(xiàng)且刪除無(wú)統(tǒng)計(jì)學(xué)意義的自變量edu的負(fù)二項(xiàng)分布回歸模型

proc genmod data=nbd1;

model numvisit= reform badh age loginc/noint link=log dist=nb;

run;

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

與前面的結(jié)果相比，擬合優(yōu)度評(píng)價(jià)指標(biāo)的取值又有所改善。

最大似然參數(shù)估計(jì)值的分析參數(shù)自由度估計(jì)值標(biāo)準(zhǔn)誤差Wald 95%置信限Wald卡方Pr >卡方Intercept00.00000.00000.00000.0000..reform1-0.13970.0510-0.2398-0.03977.490.0062badh11.13090.07450.98481.2769230.29<0.0001age10.00560.00240.00100.01035.710.0169loginc10.07640.01190.05310.099741.17<0.0001離散度11.00290.04770.91371.1008

基于以上的結(jié)果可以寫(xiě)出最終的負(fù)二項(xiàng)分布回歸模型的“內(nèi)核”，如下：

ln[μ(X)]=-0.1397reform+1.1309badh+0.0056age+0.0764loginc

μ(X)=e-0.1397reform+1.1309badh+0.0056age+0.0764loginc

若將上面“μ(X)”等號(hào)后面的內(nèi)容代入式(6)或式(7)，便可獲得多重負(fù)二項(xiàng)分布回歸模型的完整表達(dá)式。

【專(zhuān)業(yè)結(jié)論】改革與否、健康狀況和年齡這三個(gè)因素對(duì)患者就診次數(shù)的影響具有統(tǒng)計(jì)學(xué)意義，結(jié)合參數(shù)估計(jì)的結(jié)果來(lái)看：改革與否的系數(shù)為-0.1397，說(shuō)明改革之后，患者的就診次數(shù)下降了。具體地說(shuō)，因exp(-0.1397)=0.869619，即就診次數(shù)下降了約1-0.869619=13.04%；相對(duì)于良好健康狀況、年齡小的患者，具有不良健康狀況、年齡大均會(huì)使就診次數(shù)增加(因?yàn)樗鼈兊幕貧w系數(shù)為正值)。