復(fù)雜抽樣Poisson回歸分析方法及應(yīng)用*

胡躍華 匡翔宇 金承剛 Hasanat Alamgir 馬林茂 馮國(guó)雙 于石成△

最常見的復(fù)雜抽樣是多階段抽樣，即在抽取樣本時(shí)分為兩個(gè)或兩個(gè)以上階段從總體中抽取樣本的一種抽樣方法。對(duì)這種采用多階段抽樣設(shè)計(jì)得到的數(shù)據(jù)，如果應(yīng)用傳統(tǒng)的統(tǒng)計(jì)方法分析數(shù)據(jù)而沒有考慮由分層、整群及不等概率等因素所造成的設(shè)計(jì)效應(yīng)變化，常常會(huì)導(dǎo)致統(tǒng)計(jì)推斷的錯(cuò)誤。主要體現(xiàn)在參數(shù)估計(jì)的標(biāo)準(zhǔn)誤可能會(huì)被低估，估計(jì)的可信區(qū)間偏窄;假設(shè)檢驗(yàn)的I類錯(cuò)誤概率大于所規(guī)定的檢驗(yàn)水準(zhǔn)α〔1〕。

自1974年Kish和Frankel發(fā)表了針對(duì)復(fù)雜抽樣調(diào)查數(shù)據(jù)的統(tǒng)計(jì)推斷方法以來(lái)〔1〕，近年來(lái)復(fù)雜抽樣調(diào)查數(shù)據(jù)的分析方法越來(lái)越受到關(guān)注〔2-3〕。大多數(shù)統(tǒng)計(jì)分析軟件，如 SAS，SPSS，STATA和 SUDAN，將抽樣權(quán)重納入統(tǒng)計(jì)分析過程進(jìn)行描述性統(tǒng)計(jì)、交叉表、復(fù)雜抽樣一般線性模型和logistic回歸等，使得復(fù)雜抽樣調(diào)查數(shù)據(jù)分析應(yīng)用越來(lái)越普及。

本文針對(duì)美國(guó)德克薩斯州行為危險(xiǎn)因素監(jiān)測(cè)系統(tǒng)(BRFSS)數(shù)據(jù)，探討了45歲及以上的人群體重指數(shù)(BMI)與跌倒性傷害發(fā)生的聯(lián)系;比較了復(fù)雜抽樣Poisson回歸與普通Poisson回歸分析的結(jié)果，旨在闡述復(fù)雜抽樣調(diào)查數(shù)據(jù)應(yīng)用復(fù)雜抽樣統(tǒng)計(jì)推斷方法的必要性。

模型基本原理

復(fù)雜抽樣Poisson回歸也稱調(diào)查加權(quán)Poisson回歸，適用于用多階段抽樣技術(shù)得到的數(shù)據(jù)，如分層隨機(jī)抽樣或整群抽樣。與普通Poisson回歸一樣，復(fù)雜抽樣Poisson回歸的反應(yīng)變量作為一組解釋變量的函數(shù)，要求在一定時(shí)間或空間上事件的發(fā)生獨(dú)立、總體均數(shù)與方差相等以及觀察結(jié)果有可加性;不同的是復(fù)雜抽樣Poisson回歸在統(tǒng)計(jì)推斷時(shí)考慮了使用調(diào)查設(shè)計(jì)的信息來(lái)校正方差估計(jì)值〔4〕。反應(yīng)變量表示在某一特定時(shí)間或空間發(fā)生的獨(dú)立事件數(shù)目，它是一個(gè)非負(fù)的整數(shù)。事件發(fā)生數(shù)目這一隨機(jī)變量服從Poisson分布:

λi是一定時(shí)間ti事件發(fā)生的“預(yù)期數(shù)或率”。Poisson回歸模型應(yīng)用“預(yù)期數(shù)或率”自然對(duì)數(shù)連接函數(shù)來(lái)構(gòu)造與解釋變量x的線性函數(shù)關(guān)系，如下:

xi是解釋變量向量，β是參數(shù)向量，p為解釋變量的數(shù)目。Poisson隨機(jī)變量Yi的重要特性之一是其均數(shù)等于方差，即:

將式(1)進(jìn)行對(duì)數(shù)變換，得:

將(3)式代入(2)式，得:

當(dāng)構(gòu)造應(yīng)變量期望值E(yi|xi)的線性模型時(shí)，對(duì)式(4)應(yīng)用自然對(duì)數(shù)連接轉(zhuǎn)換，得到:

log(ti)為第i個(gè)個(gè)體觀察時(shí)間的自然對(duì)數(shù)值，稱為偏移量;上式為標(biāo)準(zhǔn)的Poisson回歸模型。參數(shù)估計(jì)及檢驗(yàn):復(fù)雜抽樣 Poisson回歸參數(shù)估計(jì)構(gòu)造一個(gè)加權(quán)Poisson偽似然函數(shù)(weighted Poisson pseudo-likelihood function)為:

采用牛頓-拉夫遜最大算法(Newton-Raphson maximum algorithm)使上述似然函數(shù)達(dá)到最大值的參數(shù)估計(jì)值，即是復(fù)雜抽樣Poisson回歸模型參數(shù)估計(jì)。復(fù)雜抽樣Poisson回歸參數(shù)假設(shè)檢驗(yàn)用學(xué)生t檢驗(yàn)，Stata提供了設(shè)計(jì)調(diào)整的Wald檢驗(yàn)(design-adjusted Wald tests)進(jìn)行多個(gè)解釋變量參數(shù)的檢驗(yàn)〔5〕。

Poisson回歸整個(gè)模型的擬合評(píng)價(jià)采用Pearson卡方檢驗(yàn)或偏差統(tǒng)計(jì)量及 Deviance殘差圖〔6〕。如果Pearson卡方檢驗(yàn)統(tǒng)計(jì)量或偏差統(tǒng)計(jì)量與其自由度的比值等于1或接近1，P＞0.05，整個(gè)模型擬合優(yōu)度好;如果Pearson卡方檢驗(yàn)統(tǒng)計(jì)量或偏差統(tǒng)計(jì)量與其自由度的比值遠(yuǎn)大于1，說明計(jì)數(shù)數(shù)據(jù)間的變異較大，假定分布的均數(shù)與方差的關(guān)系不正確，Poisson回歸不適合用來(lái)描述該數(shù)據(jù)，也稱為過度離散(over-dispersion)。過度離散的控制可在模型中引入over-dispersion參數(shù)或采用負(fù)二項(xiàng)分布(negative binominal distribution)，負(fù)二項(xiàng)分布可處理資料不獨(dú)立造成的過度離散，如具有傳染性、地方聚集性、家庭聚集性(如乙肝)等疾病〔7-8〕。Deviance殘差圖一般用來(lái)比較兩個(gè)回歸模型的擬合優(yōu)度，如Poisson回歸和負(fù)二項(xiàng)回歸;如果模型擬合較好，其殘差絕對(duì)值在Deviance殘差圖中較小，有向0點(diǎn)收縮的趨勢(shì)〔6〕。

Stata沒有提供復(fù)雜抽樣Poisson回歸或負(fù)二項(xiàng)回歸的擬合優(yōu)度(goodness-of-fit)的Pearson卡方檢驗(yàn)或偏差統(tǒng)計(jì)量及Deviance殘差圖〔9〕;可以用比較計(jì)數(shù)資料的模型分布和實(shí)際分布的圖解技術(shù)方法來(lái)判斷模型擬合優(yōu)度，讀者可參考Long和Freese的Stata程序完成比較〔10〕。但在擬合計(jì)數(shù)數(shù)據(jù)模型的時(shí)候，一個(gè)非常重要的步驟是確定Poisson回歸和負(fù)二項(xiàng)回歸哪個(gè)模型能更好地?cái)M合數(shù)據(jù)。可在擬合一個(gè)復(fù)雜抽樣Poisson回歸模型后，再擬合一個(gè)有相同解釋變量x的負(fù)二項(xiàng)回歸，并計(jì)算負(fù)二項(xiàng)回歸的離散參數(shù)α。Stata沒有給出α的似然比檢驗(yàn)，但輸出α的點(diǎn)估計(jì)值和95%可信區(qū)間;如果α的95%可信區(qū)間包含0，表示可用復(fù)雜抽樣Poisson回歸擬合復(fù)雜抽樣調(diào)查數(shù)據(jù);否則改用復(fù)雜抽樣負(fù)二項(xiàng)回歸估計(jì)模型參數(shù)。

應(yīng)用實(shí)例及結(jié)果解釋

行為危險(xiǎn)因素監(jiān)測(cè)系統(tǒng)是美國(guó)CDC每年開展的行為危險(xiǎn)因素電話調(diào)查，調(diào)查對(duì)象為年滿18周歲及以上成年人，目前全美50個(gè)州都進(jìn)行此項(xiàng)調(diào)查。由于BRFSS在樣本選取時(shí)應(yīng)用了抽取家庭電話號(hào)碼的多階段抽樣，分析時(shí)要考慮電話號(hào)碼選取概率不同、家庭中成年人數(shù)、家庭電話線條數(shù)和年齡性別的事后分層權(quán)重調(diào)整，將估計(jì)值調(diào)整到總體人群。因此，這種調(diào)整解決了每個(gè)人的抽樣概率不等或樣本對(duì)總體特征缺乏代表性的問題。

德克薩斯州BRFSS調(diào)查采用不成比例分層抽樣(disproportionate stratified sample，DSS)設(shè)計(jì)。目標(biāo)人群是德克薩斯州所有年滿18周歲或以上的居民，且它們家庭的電話號(hào)碼包括在用于抽樣的電話號(hào)碼簿里，即抽樣框。將家庭電話號(hào)碼分為兩層，高使用頻度和中等使用頻度分別抽樣，高使用頻度的電話號(hào)碼抽樣概率大。抽到合格家庭后，如果家庭只有一個(gè)人且合格，則調(diào)查該人;如果有兩個(gè)或以上合格的調(diào)查對(duì)象，則最近過生日的人為調(diào)查對(duì)象〔11〕。美國(guó)CDC已對(duì)德克薩斯州BRFSS數(shù)據(jù)計(jì)算出了每個(gè)個(gè)體的最終權(quán)重。

本報(bào)告應(yīng)用2010年美國(guó)德克薩斯州BRFSS數(shù)據(jù)，分析跌倒后造成的傷害與體重指數(shù)的關(guān)系。因變量傷害定義為在過去的三個(gè)月內(nèi)由于跌倒限制了至少一天的日?；顒?dòng)或去看了醫(yī)生。自變量有體重指數(shù)(正常、超重、肥胖)、年齡(歲)、健康狀況(中等和差、好和非常好)、性別、是否患有心血管疾病、婚姻狀況、受教育程度、雇傭狀況，其變量說明見表1。

BRFSS數(shù)據(jù)將在過去的3個(gè)月內(nèi)發(fā)生1次或2次及以上跌倒性傷害作為反應(yīng)變量，記為1，未發(fā)生任何跌倒性傷害為0，反應(yīng)變量為二分類變量，可用復(fù)雜抽樣logistic回歸分析該數(shù)據(jù)，但這樣分析損失了跌倒性傷害次數(shù)的信息。Poisson回歸可處理計(jì)數(shù)的結(jié)果變量，加上本例跌倒性傷害發(fā)生率5.8%，可用復(fù)雜抽樣Poisson回歸分析數(shù)據(jù)。

SAS未提供擬合復(fù)雜抽樣Poisson回歸程序，分析用Stata完成。首先，設(shè)置復(fù)雜抽樣設(shè)計(jì)，其Stata語(yǔ)句為:

wt代表BRFSS數(shù)據(jù)庫(kù)個(gè)體的最終權(quán)重，a表示數(shù)據(jù)庫(kù)的分層變量，b說明數(shù)據(jù)庫(kù)的群或基本抽樣單位;這3個(gè)變量都需替換為BRFSS數(shù)據(jù)庫(kù)里的相應(yīng)變量名。Stata擬合Poisson回歸、復(fù)雜抽樣Poisson回歸和復(fù)雜抽樣NB回歸語(yǔ)句:

injury為反應(yīng)變量，其他為解釋變量;在 bmi和health等解釋變量前面的i指示該解釋變量為分類變量，以啞變量的形式進(jìn)入回歸模型;irr要求輸出RR值及其95%可信區(qū)間。

模型擬合策略:體重指數(shù)(正常、超重、肥胖)、年齡(歲)、健康狀況(中等或差、好或非常好)、性別、是否患有糖尿病、是否患有心血管疾病、婚姻狀況、受教育程度、家庭收入、雇傭狀況、居住地點(diǎn)和種族12個(gè)變量作單因素復(fù)雜抽樣Poisson回歸分析，變量總的P值顯示性別、居住地點(diǎn)、種族與跌倒性傷害無(wú)聯(lián)系，將不納入多因素分析的模型;從專業(yè)角度考慮，糖尿病與心血管疾病有共線性，將與跌倒性傷害聯(lián)系強(qiáng)的心血管疾病納入多因素分析模型中。分析超重、肥胖與跌倒性傷害的關(guān)聯(lián)時(shí)，將年齡、健康狀況、性別、心血管疾病、婚姻狀況、受教育程度、雇傭狀況作為調(diào)整變量放入模型，擬合普通Poisson回歸、復(fù)雜抽樣Poisson回歸和復(fù)雜抽樣NB回歸模型，結(jié)果見表1。

表1 普通Poisson回歸、復(fù)雜抽樣Poisson回歸和復(fù)雜抽樣NB回歸擬合結(jié)果

對(duì)于復(fù)雜抽樣Poisson回歸或復(fù)雜抽樣NB回歸，目前Stata沒有給出類似普通Poisson回歸的似然比檢驗(yàn);只有應(yīng)用復(fù)雜抽樣NB回歸計(jì)算離散參數(shù)，確定Poisson回歸或NB回歸哪個(gè)更適合描述數(shù)據(jù)。從表1的結(jié)果看，數(shù)據(jù)擬合負(fù)二項(xiàng)回歸得離散參數(shù)為8.15，標(biāo)準(zhǔn)誤1.68，α的95%CI=5.44-12.22，表明跌倒性傷害數(shù)據(jù)的方差是其均數(shù)的8.15倍，因此復(fù)雜抽樣Poisson回歸不適合擬合該數(shù)據(jù)，可應(yīng)用復(fù)雜抽樣負(fù)二項(xiàng)回歸模型。普通Poisson回歸分析與復(fù)雜抽樣NB回歸相比，前者高估了肥胖、患有心血管疾病、小于25000美元的年收入與跌倒性傷害的聯(lián)系;低估了中等和差的健康狀況、未被雇用與跌倒性傷害的聯(lián)系;雖然也低估了已婚和受教育程度低與跌倒性傷害的聯(lián)系，但擬合復(fù)雜抽樣NB回歸時(shí)，其P值小于0.05，而在普通Poisson回歸中，P值小于0.001;年齡效應(yīng)大小沒用改變，其統(tǒng)計(jì)學(xué)假設(shè)檢驗(yàn)結(jié)果也保持未變。

討 論

在大型流行病學(xué)調(diào)查中多階段抽樣方法較簡(jiǎn)單隨機(jī)抽樣方法更經(jīng)濟(jì)有效，因此常常被使用。如將調(diào)查總體用某一特征(地理位置、經(jīng)濟(jì)狀況、城市農(nóng)村等)分為不同的層，然后從每一層里選取“群”(clusters of subjects)，再?gòu)倪x擇的群里抽取一定數(shù)量的個(gè)體，同一“群”里的個(gè)體更具有相似性。對(duì)這種多階段抽樣得到的數(shù)據(jù)，如果用簡(jiǎn)單隨機(jī)抽樣的常規(guī)分析方法來(lái)分析復(fù)雜抽樣的數(shù)據(jù)，因個(gè)體間的相似程度大，參數(shù)標(biāo)準(zhǔn)誤估計(jì)通常會(huì)偏小?？紤]“群”相似性特征的復(fù)雜抽樣調(diào)查數(shù)據(jù)加權(quán)分析可產(chǎn)生更準(zhǔn)確的標(biāo)準(zhǔn)誤估計(jì)。

對(duì)于德克薩斯州BRFSS數(shù)據(jù)，普通Poisson回歸分析結(jié)果顯示:除超重vs.正常和年收入＜$25000 vs.$75，000+外，其他解釋變量均有統(tǒng)計(jì)學(xué)意義。普通Poisson回歸分析假設(shè)調(diào)查數(shù)據(jù)來(lái)自簡(jiǎn)單隨機(jī)抽樣(SRS)，這導(dǎo)致低估模型參數(shù)的方差和標(biāo)準(zhǔn)誤，進(jìn)而可信區(qū)間的估計(jì)變窄，致使統(tǒng)計(jì)推斷受到影響，將無(wú)統(tǒng)計(jì)學(xué)意義判斷為有統(tǒng)計(jì)學(xué)差異。如參數(shù)檢驗(yàn)，普通Poisson回歸估計(jì)的Z值偏大，P值偏小;作統(tǒng)計(jì)推斷時(shí)，增加了I類錯(cuò)誤的風(fēng)險(xiǎn)。

復(fù)雜抽樣調(diào)查數(shù)據(jù)處理的理論和方法可追溯到100年前〔12-13〕，1945～1975年創(chuàng)立的復(fù)雜抽樣調(diào)查的抽樣設(shè)計(jì)、總體參數(shù)估計(jì)和統(tǒng)計(jì)推斷仍是目前復(fù)雜抽樣調(diào)查數(shù)據(jù)描述分析的基石。上世紀(jì)40年代后期，美國(guó)哥倫比亞大學(xué)社會(huì)學(xué)家Paul Lazarsfeld建立了復(fù)雜抽樣調(diào)查數(shù)據(jù)測(cè)量變量間聯(lián)系的分析方法，而不限于僅是總體特征的描述〔14〕。在1950～1990年隨著新統(tǒng)計(jì)理論和方法的發(fā)展，復(fù)雜抽樣調(diào)查數(shù)據(jù)的分析處理方法得到了更快的發(fā)展，包括列聯(lián)表資料的對(duì)數(shù)線性模型和相關(guān)的方法、廣義線性模型(logistic和Poisson回歸)、生存分析、一般線性混合模型、結(jié)構(gòu)方程模型和隱變量模型〔15〕，這些方法充分考慮了復(fù)雜抽樣調(diào)查設(shè)計(jì)的樣本特性，對(duì)參數(shù)估計(jì)的標(biāo)準(zhǔn)誤和可信區(qū)間有更準(zhǔn)確的估計(jì)。

目前，SAS 9.2已包含復(fù)雜抽樣調(diào)查數(shù)據(jù)的描述(PROC SURVEYMEANS和 PROC SURVEYFREQ)、多元線性回歸(PROC SURVEYREG)和logistic回歸分析(PROC SURVEYLOGISTIC);但仍沒有復(fù)雜抽樣Poisson回歸分析模塊。Stata Version 13包含了描述分析、二分類結(jié)果(logistic回歸和Probit回歸)、有序結(jié)果(有序logistic回歸和有序Probit回歸)，多分類結(jié)果(多項(xiàng)式 logistic回歸和 Probit回歸)，計(jì)數(shù)結(jié)果(Poisson回歸，負(fù)二項(xiàng)回歸，零膨脹Poisson回歸和零膨脹負(fù)二項(xiàng)回歸)和生存分析等。SPSS最新版本除有基本的描述分析外，還有多元線性回歸、二分類logistic回歸、有序logistic回歸和Cox回歸分析等。SUDAN Version 9.0是一個(gè)很好的處理復(fù)雜抽樣調(diào)查數(shù)據(jù)的統(tǒng)計(jì)分析軟件，它包括描述性統(tǒng)計(jì)(PROC DESCRIPT，PROC CROSSTAB，PROC RATIO)、多元線性回歸(PROC REGRESS)、二分類 logistic回歸(PROC LOGISTIC/RLOGIST)、有序和多分類logistic回歸(PROC MULTILOG)、Poisson回歸(PROC LOGLINK)和生存分析(PROC KAPMEIER，PROC SURVIVAL)。除了上述分析復(fù)雜抽樣調(diào)查數(shù)據(jù)四大軟件外，其他統(tǒng)計(jì)軟件如M-Plus，R，IVEware和WesVar等軟件也具有復(fù)雜抽樣調(diào)查數(shù)據(jù)統(tǒng)計(jì)描述和統(tǒng)計(jì)分析的功能，感興趣的讀者可參考有關(guān)書籍。

1.Kish L，F(xiàn)rankel MR.Inference from complex samples.Journal of the Royal Statistical Society，series B，1974，36:1-37.

2.Rao JNK.Interplay between sample survey theory and practice;an appraisal.Survey Methodology，2005，31:117-138.

3.劉建華，金水高.復(fù)雜抽樣調(diào)查總體特征量及其方差的估計(jì).中國(guó)衛(wèi)生統(tǒng)計(jì)，2008，25(4):377-379.

4.Binder DA，On the variances of asymptotically normal estimators from complex surveys.Survey Methodology，1981，7:157-170.

5.Heeringga S，O'Muicheartaigh C，Survey Methods in Multinational Contexts.247-263.

6.楚慧珠，郜艷輝，鄒宇華，等.負(fù)二項(xiàng)回歸和Poisson回歸在改水降氟效果中的對(duì)比分析.數(shù)理醫(yī)藥學(xué)雜志，2008，21(6):655-657.

7.劉亞杰，李海波，潘萍，等.朝陽(yáng)市某高級(jí)中學(xué)結(jié)核病爆發(fā)的班級(jí)聚集性分析.中國(guó)衛(wèi)生統(tǒng)計(jì)，2010，27(4):371.

8.薛付忠，王潔貞，馬希蘭.疾病空間分布狀態(tài)的負(fù)二項(xiàng)分布概率生成模型的討論.中國(guó)衛(wèi)生統(tǒng)計(jì)，2000，17(6):366-368.

9.Heeringa SG，Alcser KH，Doerr K，et al.Potential selection bias in a community-based study of PSA Levels in African-American men.Journal of Clinical Epidemiology，2001，54(2):142-148.

10.Long JS，F(xiàn)reese J.Regression Models for Categorical Dependent Variables Using Stata.2nd ed，Stat Press，College Station，Texas，2006.

11.http://www.cdc.gov/brfss/.

12.Bowley AL.Address to the Economic Science and Statistics Section of the British Association for the Advancement of Science.Journal of the Royal Statistical Society，1906，69:548-557.

13.Fisher RA.Statistical Methods for Research Workers.Oliver and Boyd，Edinburgh，1925.

14.Kendall PL，Lazarsfeld PE.problems of survey analysis，in R.K.Merton and PF Lazaarsfeld(Eds.)，Continuities in Social Research:Studies in the Scope and Method of“The American Solders，”.Free Press，Chicago，1950.

15.Binder DA.On the variances of asymptotically normal estimation from complex surveys.Survey Methodology，1981，7:157-170.