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

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

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

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

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

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

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

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

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

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

模型原理與方法

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

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

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

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

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

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

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

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

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

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

但此時，

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

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

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

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

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

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

實例分析

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

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

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

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

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

討 論

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

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

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

除基于穩(wěn)健Poisson模型的多水平和GEE模型外，在Log-binomial模型基礎(chǔ)上擴展的多水平模型和GEE模型也可處理具有層次結(jié)構(gòu)特征的常見結(jié)局資料。與多水平穩(wěn)健Poisson模型和穩(wěn)健Poisson-GEE模型的比較研究也是我們進一步關(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)系的二分類資料處理方法比較．中國衛(wèi)生統(tǒng)計，2010，10(27):464 －467．

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

附錄:

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

proc genmod;

class ID;/*ID為個體的標識變量*/

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

repeated subject=ID;

run;

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

proc genmod;

class K;/*K為類體的標識變量*/

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為類體的標識變量*/

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é)任編輯:丁海龍)