胡良平

(1.軍事科學(xué)院研究生院，北京 100850； 2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會，北京 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)見表1[1]。

【對數(shù)據(jù)結(jié)構(gòu)的分析】因變量Y為“計(jì)數(shù)變量”，其均值為7.467、方差為13.016，方差稍大于均值，只能將Y近似視為服從Poisson分布的離散型隨機(jī)變量(理論上，均值應(yīng)等于方差)；擬考察的三個(gè)自變量均是“二值變量”，在回歸分析中是可以接受的，但在統(tǒng)計(jì)學(xué)理論中，默認(rèn)自變量為“計(jì)量變量”。

1.1.2 Poisson分布概率函數(shù)[2]

設(shè)Y是一個(gè)服從Poisson分布的隨機(jī)變量，X=(1,x1,x2,…,xq)'是一個(gè)協(xié)變量向量，β=(β0,β1,β2,…,βq)'是參數(shù)向量。如果Y的數(shù)學(xué)期望的對數(shù)可以表示為協(xié)變量的線性表達(dá)式，即E(Y|X)=exp(X'β)=或ln[E(Y|X)]=X'β，則稱(Y，X)服從Poisson分布回歸模型。對應(yīng)的表達(dá)式見下面的式(1)：

(1)

表1 三個(gè)可能的影響因素(X1-X3)對30例非氣質(zhì)性心臟病且僅有胸悶癥狀就診者24小時(shí)早搏數(shù)Y的觀察結(jié)果

注：Y，24小時(shí)中早搏數(shù)；X1，是否吸煙(1-吸煙、0-不吸煙)；X2，是否喝咖啡(1-喝、0-不喝)；X3，性別(1-男、0-女)

1.1.3 一般計(jì)數(shù)資料Poisson分布回歸模型的表達(dá)式

(2)

值得注意的是：統(tǒng)計(jì)軟件給出的計(jì)算結(jié)果是下面式(3)等號右邊回歸系數(shù)的估計(jì)值。也就是說，Poisson分布回歸模型的“核心部分”是其“期望值(或稱均值)”，用自變量的線性組合來表達(dá)計(jì)數(shù)因變量Y的期望的自然對數(shù)。

ln[λ(X)]=β0+β1x1+β2x2+…+βmxm+ε

(3)

式(3)中的“l(fā)n”代表“自然對數(shù)函數(shù)”，λ(X)代表服從Poisson分布的離散型隨機(jī)變量Y的均值，X代表所有自變量構(gòu)成的“向量”，ε代表模型的誤差。有些統(tǒng)計(jì)學(xué)教科書將式(3)稱為“Poisson分布回歸模型”，這是不夠準(zhǔn)確的。它應(yīng)被視為“Poisson分布回歸模型”的“內(nèi)核”，即“核心內(nèi)容”，外文文獻(xiàn)中稱其為“Poisson分布回歸模型”的“連接函數(shù)”。

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

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

如何求解出式(3)中的回歸系數(shù)呢？首先，應(yīng)該明白：λ(X)是X的函數(shù)，在實(shí)際資料中是無法觀測到的，它隨著X取值的改變而變化，是一個(gè)“隱變量”；回憶求解多重線性回歸模型中的參數(shù)時(shí)所面臨的“處境”，在那里，計(jì)量因變量Y的取值是可以觀測到的，是“顯變量”[3-6]。

顯然，直接根據(jù)式(3)無法求出式中的回歸系數(shù)。需要基于式(2)并采用“最大似然法”[7]構(gòu)造“目標(biāo)函數(shù)”，運(yùn)用高等數(shù)學(xué)方法，將“目標(biāo)函數(shù)”轉(zhuǎn)變成“方程組”，再運(yùn)用“計(jì)算數(shù)學(xué)”或“統(tǒng)計(jì)計(jì)算”等技術(shù)方法求解“方程組”。由于方程組為“非線性方程組”，一般沒有直接解，需用到迭代計(jì)算法，如加權(quán)最小二乘迭代或牛頓-納法生迭代法。從而獲得式(3)的“參數(shù)估計(jì)”結(jié)果。因篇幅所限，參數(shù)估計(jì)公式的詳細(xì)推導(dǎo)過程從略。

1.2.2 Possion回歸系數(shù)的假設(shè)檢驗(yàn)

1.2.2.1 似然比檢驗(yàn)

通過比較兩個(gè)相嵌套模型(如模型P嵌套于模型K內(nèi))的對數(shù)似然函數(shù)統(tǒng)計(jì)量G(又稱Deviance)來進(jìn)行，其統(tǒng)計(jì)量見式(4)：

G=GP-GK=-2× (模型P的對數(shù)似然函數(shù)-模型K的對數(shù)似然函數(shù))

(4)

其中，模型P中的變量是模型K中變量的一部分，另一部分就是需要檢驗(yàn)的變量。這里，G服從自由度為K-P的χ2分布。

1.2.2.2 回歸系數(shù)的Wald檢驗(yàn)

比較估計(jì)系數(shù)與“0”的差別是否有統(tǒng)計(jì)學(xué)意義，其檢驗(yàn)統(tǒng)計(jì)量見式(5)：

(5)

這里z為標(biāo)準(zhǔn)正態(tài)變量。參數(shù)的置信區(qū)間是基于Wald統(tǒng)計(jì)量導(dǎo)出的。β的95%置信區(qū)間見式(6)：

(6)

1.2.3 Possion擬合優(yōu)度檢驗(yàn)

1.2.3.1 Deviance偏差統(tǒng)計(jì)量

設(shè)L(bmax;y)為全模型的似然函數(shù)，L(b;y)為選模型(含部分自變量的模型)的似然函數(shù)，則定義λ統(tǒng)計(jì)量見下面的式(7)：

λ=L(bmax;y)/L(b;y)

(7)

計(jì)算Deviance偏差(D)統(tǒng)計(jì)量，采用公式(8)：

D=2lnλ=2[lnL(bmax;y)-lnL(b;y)]

(8)

且D服從χ2分布。很顯然，D越大，模型的估計(jì)值與觀測值的偏差就越大，擬合效果就越差。

1.2.3.2 廣義χ2統(tǒng)計(jì)量

廣義χ2統(tǒng)計(jì)量見式(9)：

(9)

2 一般計(jì)數(shù)資料Poisson分布回歸模型的SAS實(shí)現(xiàn)

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

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

data maop1006;

input obs X1-X3 Y @@;

cards;

1011111610011 2000 717011 8 3000 318101 9 4101 519000 8 5000 220100 5 61111321011 5 7010 622110 8 8101102311013 9000 424001 810101 725100 611000 126001 412001 927000 613001 6281111314111 1729110 915000 530001 5;run;

2.2 求出因變量Y的均值與方差

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

proc univariate data=maop1006 noprint;

var Y;

output out=aaa mean=ybar var=yvar;

run;

proc print data=aaa;

var ybar yvar;

run;

【SAS輸出結(jié)果】

Obsybaryvar17.4666713.0161

求出Y的均值和方差分別為7.467、13.016。

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

利用下面的SAS過程步構(gòu)建Poisson分布回歸模型：

proc genmod data=maop1006;

model Y=X1-X3/link=log dist=poisson;

run;

【SAS程序說明】“l(fā)ink=log”表明采用“自然對數(shù)”作為連接函數(shù)，即對計(jì)數(shù)因變量取自然對數(shù)變換；“dist=poisson”表明要基于Poisson分布構(gòu)建Poisson分布回歸模型。

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

最大似然參數(shù)估計(jì)值的分析

以上結(jié)果表明：X3(性別)對因變量的影響無統(tǒng)計(jì)學(xué)意義。此結(jié)果是否正確，有待進(jìn)一步研究。

2.4 在構(gòu)建多重Poisson分布回歸模型時(shí)檢驗(yàn)因變量是否存在“過離散”現(xiàn)象

利用下面的SAS過程步程序構(gòu)建Poisson分布回歸模型，重點(diǎn)在于檢驗(yàn)因變量是否存在“過離散”現(xiàn)象：

proc genmod data=maop1006;

model Y=X1-X3/link=log dist=nb noscale;

run;

【SAS程序說明】“dist=nb”說明指定資料中因變量(誤差項(xiàng))的分布為“負(fù)二項(xiàng)分布”；關(guān)鍵問題在于：只有指定為此分布時(shí)，才能發(fā)揮選項(xiàng)“noscale”的作用，它是將冗余參數(shù)k的取值固定為“0”，相當(dāng)于進(jìn)行Poisson分布回歸分析，同時(shí)，采用拉格朗日(Lagrange)乘子檢驗(yàn)“因變量是否存在‘過離散’現(xiàn)象”。

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

第1部分輸出結(jié)果同上，從略。

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

注：*單側(cè)P值

以上結(jié)果表明：因變量確實(shí)存在“過離散”現(xiàn)象，構(gòu)建Poisson分布回歸模型時(shí)應(yīng)對其進(jìn)行校正。

2.5 在構(gòu)建多重Poisson分布回歸模型時(shí)對“過離散”現(xiàn)象進(jìn)行校正

利用下面的SAS過程步程序構(gòu)建Poisson分布回歸模型，對“過離散”現(xiàn)象進(jìn)行校正：

proc genmod data=maop1006;

model Y=X1-X3/link=log dist=poisson scale=deviance;

run;

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

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

最大似然參數(shù)估計(jì)值的分析

以上結(jié)果表明，尺度參數(shù)由“1”調(diào)整為“0.9510”，由此，模型中各參數(shù)的“標(biāo)準(zhǔn)誤”均得到“校正”。故最后兩列的數(shù)值也都得到校正。經(jīng)校正后，原來無統(tǒng)計(jì)學(xué)意義的“X3”變?yōu)橛薪y(tǒng)計(jì)學(xué)意義。由以上結(jié)果可以寫出Poisson分布回歸模型中的“內(nèi)核”：

ln[λ(X)]=1.5066+0.4162X1+0.4012X2+0.2546X3

λ(X)=e1.5066+0.4162X1+0.4012X2+0.2546X3

若將上面的λ(X)代入本文第1.1節(jié)中的式(2)，就可獲得完整的“Poisson分布回歸模型”。

【專業(yè)結(jié)論】因X1、X2和X3前的系數(shù)分別為exp(0.4162)=1.516189、exp(0.4012)=1.493616和exp(0.2546)=1.289946，又因?yàn)閄1代表“是否吸煙(1-吸煙、0-不吸煙)”、X2代表“是否喝咖啡(1-喝、0-不喝)”、X3代表“性別(1-男、0-女)”，說明吸煙者出現(xiàn)早搏的概率是不吸煙者的1.516倍、喝咖啡者出現(xiàn)早搏的概率是不喝咖啡者的1.494倍，而男性受試者出現(xiàn)早搏的概率是女性受試者的1.290倍。