穩(wěn)健Poisson模型:兩水平模型與GEE模型在相對(duì)危險(xiǎn)度或患病率比估計(jì)中的應(yīng)用比較*

周舒冬 郜艷暉△ 李麗霞 張 敏 楊 翌 陳 躍

穩(wěn)健Poisson模型:兩水平模型與GEE模型在相對(duì)危險(xiǎn)度或患病率比估計(jì)中的應(yīng)用比較*

周舒冬1郜艷暉1△李麗霞1張 敏1楊 翌 陳 躍2

目的 在處理具有層次結(jié)構(gòu)特征的非罕見結(jié)局事件資料時(shí)，比較基于穩(wěn)健Poisson模型的兩水平模型和GEE模型在估計(jì)RR/PR時(shí)的應(yīng)用。方法 將兩水平穩(wěn)健Poisson模型及穩(wěn)健Poisson-GEE應(yīng)用到2010年歐洲社會(huì)調(diào)查資料，估計(jì)影響居民生活滿意度的各因素相關(guān)的PR及95%CI，以說(shuō)明兩模型在理論和應(yīng)用上的區(qū)別和聯(lián)系。結(jié)果 穩(wěn)健Poisson-GEE模型的PR估計(jì)值與穩(wěn)健Poisson回歸模型相同，但置信區(qū)間較寬;兩水平穩(wěn)健Poisson模型的PR值較GEE模型為低，顯示了隨機(jī)效應(yīng)對(duì)解釋變量的混雜作用。結(jié)論 兩種方法均可處理具有層次結(jié)構(gòu)特征的非罕見結(jié)局事件的RR/PR估計(jì)，但兩水平模型比GEE可提供更多隨機(jī)效應(yīng)的信息，且易于擴(kuò)展至更高水平或隨機(jī)系數(shù)模型。

層次結(jié)構(gòu) 非罕見結(jié)局 穩(wěn)健Poisson回歸 廣義估計(jì)方程 相對(duì)危險(xiǎn)度 患病率比

1．廣東藥學(xué)院公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系，廣東省分子流行病學(xué)重點(diǎn)實(shí)驗(yàn)室(510310)

2．Department of Epidemiology and Community Medicine，University of Ottawa

△通信作者:郜艷暉，E-mail:gao_yanhui@163.com

由于logistic回歸模型已被廣大研究者所熟識(shí)，因此很多文獻(xiàn)習(xí)慣性地將其應(yīng)用于橫斷面研究乃至隊(duì)列研究，計(jì)算優(yōu)勢(shì)比(odds ratio，OR)，作為相對(duì)危險(xiǎn)度(relative risk，RR)或患病率比(prevalence ratio，PR)的估計(jì)并給予同樣解釋。但當(dāng)研究結(jié)局頻率較高時(shí)，OR值作為RR/PR的估計(jì)嚴(yán)重地高估暴露因素對(duì)結(jié)局的影響〔1〕，為此統(tǒng)計(jì)學(xué)者建議宜用Log-binomial模型或穩(wěn)健 Poisson回歸模型等方法直接計(jì)算RR/PR〔2〕。Log-binomial模型屬于廣義線性模型，采用log鏈接函數(shù)，誤差分布為二項(xiàng)分布。但該法的缺陷在于最大似然估計(jì)(maximum likelihood estimate，MLE)當(dāng)參數(shù)落在所限制范圍的邊界，特別是模型中含有連續(xù)型協(xié)變量時(shí)，常得不到似然函數(shù)導(dǎo)數(shù)為零的極大值，導(dǎo)致模型無(wú)法收斂〔3－4〕，后有學(xué)者提出采用 COPY 算法〔5〕解決模型不收斂的問(wèn)題。而穩(wěn)健 Poisson回歸模型〔1，6〕指定誤差分布為Poisson分布，應(yīng)用“三明治”法獲得合理的方差估計(jì)(sandwich variance)，直接估計(jì)RR/PR時(shí)不存在收斂困難。在自變量均為分類變量時(shí)，估計(jì)結(jié)果與 Mantel-Haenszel分層分析法非常近似〔2，6〕。

和其他回歸模型一樣，穩(wěn)健Poisson回歸要求觀測(cè)單位間獨(dú)立。但流行病或社會(huì)學(xué)調(diào)查資料中，某些觀測(cè)單位常根據(jù)某些特征聚為一類，甚至具有多水平的特征，如采用多階段抽樣的橫斷面研究中個(gè)體來(lái)自同一街道，街道又來(lái)自同一社區(qū);或縱向研究中同一個(gè)體的多個(gè)部位同一指標(biāo)多次測(cè)量，由此獲得的數(shù)據(jù)具有明顯的層次結(jié)構(gòu)特征，表現(xiàn)為類間獨(dú)立、類內(nèi)相關(guān)的特性，從而違背傳統(tǒng)回歸模型的應(yīng)用條件。近年來(lái)廣泛應(yīng)用的廣義估計(jì)方程(generalized estimating equation，GEE)和多水平模型(multi-level models)是解決層次結(jié)構(gòu)數(shù)據(jù)的有力工具，但基于穩(wěn)健Poisson回歸的相關(guān)理論和應(yīng)用研究尚顯不足。本研究著力解決結(jié)局變量頻率較高且存在類內(nèi)相關(guān)時(shí)的RR/PR估計(jì)，比較穩(wěn)健Poisson-GEE模型和多水平穩(wěn)健Poisson模型在該類特征資料中的應(yīng)用。

模型原理與方法

1.穩(wěn)健Poisson回歸模型

當(dāng)結(jié)局事件頻率較高時(shí)，為直接估計(jì)RR/PR，對(duì)獨(dú)立數(shù)據(jù)，Zou〔3〕建議使用穩(wěn)健Poisson模型。設(shè)yi和Xi=(xi1，xi2，…，xiP)T分別是第i(i=1，2，…，n)個(gè)觀測(cè)的二分類結(jié)局變量和P×1維解釋變量向量，其關(guān)系可通過(guò)Poisson回歸模型表示。

式(1)中pi=Pr(yi=1/Xi)，并假設(shè)誤差分布為Poisson分布。回歸系數(shù)βp表示當(dāng)控制其它自變量后，第p個(gè)自變量xP每變化一個(gè)單位時(shí)log(p)的相應(yīng)變化。因此，與xp相對(duì)應(yīng)的相對(duì)危險(xiǎn)性為:RR(PR)=exp(βp)。

由于Poisson分布方差等于均數(shù)，當(dāng)應(yīng)用到二項(xiàng)分布資料時(shí)，易出現(xiàn)過(guò)度離散(overdispersion)問(wèn)題，導(dǎo)致參數(shù)標(biāo)準(zhǔn)誤的高估，產(chǎn)生過(guò)寬的置信區(qū)間。為此，Cameron〔7〕建議使用穩(wěn)健方差估計(jì)法，如 Huber的“三明治”方差，即:

模型(1)中參數(shù) β =(β0，β1，…，βP)T及“三明治”方差可用準(zhǔn)似然(quasi-likelihood)估計(jì)，在SAS中可用proc genmod實(shí)現(xiàn)，通過(guò)在repeated語(yǔ)句中用“subject=”指定個(gè)體編號(hào)變量(程序見附錄)。

2.穩(wěn)健Poisson-GEE模型

對(duì)非獨(dú)立二分類數(shù)據(jù)，Zou〔8〕提出仍可用穩(wěn)健Poisson回歸估計(jì)RR/PR，這時(shí)采用廣義估計(jì)方程(generalized estimated equation，GEE)的原理進(jìn)行參數(shù)估計(jì)，將類內(nèi)水平的相關(guān)作多余參數(shù)處理。

設(shè)yki和Xki=(xki1，xki2，…，xkiP)T分別為第k(k=1，2，…，K)類內(nèi)第i(i=1，2，…，nk)個(gè)個(gè)體的二分類結(jié)局變量和P×1維解釋變量向量，模型形式同(1):

式(3)中=Pr(yki=1/Xki)，回歸系數(shù)βp的涵義類似(1)。根據(jù)一致性估計(jì)方程理論，對(duì)參數(shù)β=(β0，β1，…，βP)T的“得分(score)”方程為:

方程(4)的解即為參數(shù)β的一致估計(jì)，方差矩陣形式同(2)，Var()=A－1BA－1

但此時(shí)，

式(5)中“三明治”方差的中間項(xiàng)B與(2)不同，先綜合各類內(nèi)個(gè)體的“得分”，再根據(jù)類別匯總。因此(5)可看作(2)在非獨(dú)立數(shù)據(jù)下的自然擴(kuò)展，同時(shí)校正了誤差分布的錯(cuò)誤指定與類內(nèi)響應(yīng)非獨(dú)立對(duì)參數(shù)方差估計(jì)的影響。因此該模型也可利用SAS中的proc genmod，通過(guò)在repeated語(yǔ)句里使用“subject=”指定類別變量來(lái)完成(程序見附錄)。

3.多水平穩(wěn)健Poisson模型

與穩(wěn)健Poisson-GEE模型相比，多水平穩(wěn)健Poisson模型通過(guò)納入高水平層次上的隨機(jī)效應(yīng)來(lái)處理層次數(shù)據(jù)中的類內(nèi)相關(guān)問(wèn)題。以最簡(jiǎn)單的兩水平穩(wěn)健Poisson方差分量模型為例，

式(6)中下標(biāo)和的涵義同(3)，βp為固定效應(yīng)參數(shù)，反映固定效應(yīng)xp對(duì)結(jié)局概率對(duì)數(shù)的影響，且RR(PR)=exp(βp);uk為第k類的隨機(jī)效應(yīng)，假定來(lái)自正態(tài)分布總體，即

如數(shù)據(jù)結(jié)構(gòu)中含有更高水平層次，或解釋變量在高水平上存在隨機(jī)效應(yīng)，模型(6)還可擴(kuò)展為包含更高水平上的隨機(jī)效應(yīng)，或隨機(jī)系數(shù)模型。

模型(6)的參數(shù)估計(jì)可使用SAS中的proc glimmix來(lái)完成，調(diào)用empirical選項(xiàng)用于指定“三明治”方差的結(jié)構(gòu)，默認(rèn)缺省時(shí)的“classical”即為經(jīng)典三明治方差;random語(yǔ)句可指定隨機(jī)截距或其它高水平單位上的隨機(jī)效應(yīng)。(程序見附錄)

實(shí)例分析

實(shí)例來(lái)源于2010年歐洲社會(huì)調(diào)查(European Social Survey)的開放數(shù)據(jù)(www．europeansocialsurvey．org)，該項(xiàng)目是兩年一度的多國(guó)調(diào)查，覆蓋歐洲各國(guó)。本研究選取2010年歐洲26國(guó)49024名居民的“生活滿意度”變量作為結(jié)局變量，研究自我健康評(píng)價(jià)及家庭收支情況對(duì)結(jié)局變量的影響，構(gòu)建模型時(shí)校正了年齡和性別。

49 024名居民中，對(duì)生活持滿意態(tài)度的有21 979人，約占44.83%。表1描述了各因素不同狀態(tài)下居民生活滿意度的分布情況。

表1 各因素不同狀態(tài)下居民生活滿意度的分布情況

考慮到各國(guó)居民的生活滿意度可能存在國(guó)家聚集性，因此建立多水平模型和利用廣義估計(jì)方程的方法來(lái)擬合該數(shù)據(jù)，國(guó)家作為類別指示變量。表2顯示了穩(wěn)健Poisson模型、穩(wěn)健Poisson-GEE模型及兩水平穩(wěn)健Poisson模型估計(jì)的各因素對(duì)居民生活滿意度影響的PR及95%CI。三個(gè)模型結(jié)果均顯示調(diào)整了性別和年齡后，自我健康評(píng)價(jià)和家庭收支對(duì)生活滿意度的影響均有統(tǒng)計(jì)學(xué)意義。但是和穩(wěn)健Poisson模型相比，GEE模型考慮了各國(guó)居民在生活滿意度上的國(guó)內(nèi)相關(guān)，不僅得到穩(wěn)健的PR估計(jì)值，且估計(jì)的PR置信區(qū)間較穩(wěn)健Poisson模型更寬，降低了犯I類錯(cuò)誤的風(fēng)險(xiǎn);而兩水平穩(wěn)健Poisson模型在模型構(gòu)建時(shí)添加隨機(jī)效應(yīng)，估計(jì)的PR和穩(wěn)健Poisson-GEE模型的結(jié)果不同，調(diào)整了性別和年齡后，自我健康評(píng)價(jià)和家庭收支對(duì)生活滿意度影響的PR值均低于GEE模型結(jié)果，反映了隨機(jī)效應(yīng)對(duì)解釋變量可能存在的混雜效應(yīng);從隨機(jī)效應(yīng)的方差估計(jì)值與其標(biāo)準(zhǔn)誤的比值(0.1025/0.0296=3.4628)近似推斷隨機(jī)截距項(xiàng)可能有統(tǒng)計(jì)學(xué)意義，數(shù)據(jù)的層次結(jié)構(gòu)不可忽略。此外，本文也擬合了兩水平logistic方差分量模型，調(diào)整年齡和性別后，得到兩解釋變量的OR值均高于PR值(表2)。

表2 不同模型估計(jì)各因素對(duì)居民生活滿意度影響的PR和OR及95%CI*

討 論

當(dāng)研究結(jié)局出現(xiàn)頻率較高時(shí)，將OR習(xí)慣性地解釋為RR/PR將嚴(yán)重高估暴露因素對(duì)結(jié)局的影響，這一問(wèn)題早已引起統(tǒng)計(jì)學(xué)者的注意，因此提出各種直接估計(jì)RR/PR的模型和方法，如log-binomial模型和穩(wěn)健Poisson回歸模型，其點(diǎn)估計(jì)和區(qū)間估計(jì)均比logistic回歸模型的OR解釋起來(lái)更為合理〔9〕。本文實(shí)例使用兩水平logistic模型估計(jì)自我健康評(píng)價(jià)和家庭收支對(duì)生活滿意度影響的OR值均高于穩(wěn)健Poisson-GEE和兩水平穩(wěn)健Poisson模型估計(jì)的PR值。

很多大型的流行病學(xué)調(diào)查都具有層次結(jié)構(gòu)的特征，廣義估計(jì)方程或者多水平模型是分析該類資料的兩種相對(duì)成熟的方法〔10－11〕。本文實(shí)例歐洲社會(huì)調(diào)查項(xiàng)目中，各國(guó)居民由于擁有一些共同的屬性或國(guó)內(nèi)居民間的相互影響，個(gè)體間的研究結(jié)局并不獨(dú)立。本實(shí)例擬合兩水平穩(wěn)健Poisson模型時(shí)隨機(jī)效應(yīng)參數(shù)估計(jì)結(jié)果顯示數(shù)據(jù)的層次結(jié)構(gòu)不可忽略，采用穩(wěn)健Poisson-GEE模型估計(jì)的PR值雖與穩(wěn)健Poisson模型相等，但置信區(qū)間有更寬的變化，除反映GEE模型在處理非獨(dú)立數(shù)據(jù)時(shí)能有效降低I類錯(cuò)誤的能力，某種程度上也體現(xiàn)出該數(shù)據(jù)具有較為明顯的類內(nèi)相關(guān)特征。

多水平穩(wěn)健Poisson模型和穩(wěn)健Poisson-GEE模型均可用于非獨(dú)立數(shù)據(jù)估計(jì)解釋變量對(duì)常見結(jié)局影響的RR/PR。從回歸系數(shù)的性質(zhì)看，穩(wěn)健Poisson-GEE模型中將類內(nèi)相關(guān)作多余參數(shù)，主要考慮固定效應(yīng)，因此回歸系數(shù)是群體效應(yīng)(平均效應(yīng))的回歸系數(shù);而多水平穩(wěn)健Poisson模型在隨機(jī)效應(yīng)條件下，估計(jì)的回歸系數(shù)是個(gè)體效應(yīng)的回歸系數(shù)，所以兩模型參數(shù)估計(jì)結(jié)果有時(shí)可能不同。如本研究中兩水平模型得到的各因素PR估計(jì)值均較GEE模型為低，特別是家庭收支變量，體現(xiàn)了隨機(jī)效應(yīng)可能對(duì)解釋變量的混雜作用，提示引入隨機(jī)效應(yīng)項(xiàng)后回歸系數(shù)變化較大的解釋變量也可能存在類內(nèi)相關(guān)，并且和結(jié)局變量的隨機(jī)效應(yīng)存在某種程度上的關(guān)聯(lián)。本例中根據(jù)目前結(jié)果提示每個(gè)國(guó)家的居民家庭收支狀況可能存在國(guó)家聚集性并對(duì)研究的關(guān)聯(lián)產(chǎn)生影響。在進(jìn)一步分析中可建立穩(wěn)健Poisson隨機(jī)系數(shù)模型進(jìn)行考察和比較(我們將另文研究)。此外，從模型的擴(kuò)展角度來(lái)說(shuō)，穩(wěn)健Poisson-GEE模型只能處理兩水平的資料;而多水平模型理論上可以處理更高水平結(jié)構(gòu)的資料，且可將隨機(jī)效應(yīng)分解到解釋變量上，進(jìn)一步構(gòu)建隨機(jī)系數(shù)模型，而這種分析是GEE模型目前無(wú)法完成的。

除基于穩(wěn)健Poisson模型的多水平和GEE模型外，在Log-binomial模型基礎(chǔ)上擴(kuò)展的多水平模型和GEE模型也可處理具有層次結(jié)構(gòu)特征的常見結(jié)局資料。與多水平穩(wěn)健Poisson模型和穩(wěn)健Poisson-GEE模型的比較研究也是我們進(jìn)一步關(guān)注的方向。

1．Barros AJ，Hirakata VN．Alternatives for logistic regression in cross-sectional studies:an empirical comparison of models that directly estimate the prevalence ratio．BMC Med Res Methodol，2003，3:21．

2．Petersen MR，Deddens JA．A comparison of two methods for estimating prevalence ratios．BMC Med Res Methodol，2008，8:9．

3．Lumley T，Kronmal R，Ma S:Relative risk regression in medical research:models，contrasts，estimators，and algorithms．UW Biostatistics working Paper Series．2006:293．http://www．bepress．com/uwbiostat/paper293．

4．Deddens JA，Petersen MR，Lei X．Estimation of prevalence ratios when proc genmod does not converge．In:Proceedings of the 28th Annual SAS Users Group International Conference，Paper 270 － 28．Cary，NC:SAS Institute Inc 2003．

5．Deddens JA，Petersen MR．Re:“Estimating the relative risk in cohort studies and clinical trials of common outcomes”．Am J Epidemiol，2004，159(2):213 －4;author reply 214 －5．

6．Zou G．A modified poisson regression approach to prospective studies with binary data．Am J Epidemiol，2004，159(7):702 －706．

7．Badi H，Baltagi．A companion to theoretical econometrics，Blackwell，Oxford(U．K．)，2001:331 －348．

8．Zou G，Allan D．Extension of the modified Poisson regression model to prospective studies with correlated binary data．Stat Methods Med Res，2011，11，8．［Epub ahead of print］

9．Lee J，Chia KS．Estimation of prevalence rate ratios for cross sectional data:an example in occupational epidemiology．Br J Ind Med，1993，50(9):861－862．

10．張春茂，李嬋娟，蔣志偉，等．具有相關(guān)關(guān)系的二分類資料處理方法比較．中國(guó)衛(wèi)生統(tǒng)計(jì)，2010，10(27):464 －467．

11．張華君，閔捷．廣義估計(jì)方程與多水平模型在相關(guān)資料中的比較研究．中國(guó)衛(wèi)生統(tǒng)計(jì)，2012，4(29):214 －217．

附錄:

程度1:穩(wěn)健Poisson模型SAS程序

proc genmod;

class ID;/*ID為個(gè)體的標(biāo)識(shí)變量*/

model y=X1 X2 X3/d=poisson link=1og;

repeated subject=ID;

run;

程度2:穩(wěn)健Poisson-GEE模型SAS程序

proc genmod;

class K;/*K為類體的標(biāo)識(shí)變量*/

model y=X1 X2 X3/d=poisson link=1og;

repeated subject=K;

run;

程度3:兩水平穩(wěn)健Poisson方差分量模型SAS程序

proc glimmix;

MPIRICAL=CLASSICAL;

class K;/*K為類體的標(biāo)識(shí)變量*/

model Y=X1 X2 X3/d=poisson link=1og solution;

random int/subject=K;

run;

A Comparison between Two-level and GEE Based on Robust Poisson Regression Models in the Estimation of Relative Risk or Prevalenc Rti

Zhou Shudong，Gao Yanhui，Li Lixia，et al．Department of Epidemiology and Biostatistics，School of Public Health，Guangdong Pharmaceutical University，Guangdong Key Laboratory f Mlecul Epidmig 510310，Gnzh

Objective To compare two-level and GEE based robust Poisson regression models in estimation of relative risk(RR)or prevalence ratio(PR)for common outcome data with intra-class correlation.MethodsTwo-level and GEE based robust Poisson regression models were compared by examing factors associated with life satisfaction using data from the 2010 European Social Survey．Prevalence ratios and 95%confidence intervals(95%CIs)were estimated.ResultsCompared to results from regular robust Poisson model，the GEE based robust Poisson model provided the same PR point estimates but wider 95%CIs．The two level robust Poisson model revealed lower point estimates，indicating potential confounding effects caused by random effects on the assocation of interest.ConclusionBoth two-level and GEE based methods are suitable for estimating relative risk or prevalence ratio for common outcomes with the hierarchical structure．The two-level model is superior when there are random effects，and can be easily extended for higher hierarchical structures．

Hierarchical structure;Non-rare outcome;Robust Poisson regression;Generalized estimating equations;Relative risk;Prevalence ratio

2010年廣東省自然科學(xué)基金資助(10151022401000018)

book=686,ebook=333

(責(zé)任編輯:丁海龍)