非獨立數(shù)據(jù)時穩(wěn)健Poisson和Log-binom ial隨機效應(yīng)模型研究*

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

非獨立數(shù)據(jù)時穩(wěn)健Poisson和Log-binom ial隨機效應(yīng)模型研究*

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

目的比較基于穩(wěn)健Poisson回歸及Log-binom ial的隨機效應(yīng)模型對具有層次結(jié)構(gòu)的非罕見結(jié)局事件資料的參數(shù)估計結(jié)果。方法通過計算機模擬生成具有類內(nèi)相關(guān)的數(shù)據(jù)集，設(shè)置了32種參數(shù)情況，分別反映實際資料中的少數(shù)類別數(shù)和中等類別數(shù)，小強度的類內(nèi)聚集和中等強度的類內(nèi)聚集，基線患病率分別為0.15和0.3，關(guān)聯(lián)強度PR值分別是1.5和2.0，暴露率分別為20%和50%的情況；比較各模型在模擬條件下的收斂率、參數(shù)估計結(jié)果以及參數(shù)的95%CI覆蓋率。結(jié)果Log-binomial隨機效應(yīng)模型的收斂效果劣于穩(wěn)健Poisson隨機效應(yīng)模型；兩隨機效應(yīng)模型對各參數(shù)的估計結(jié)果均和真值較為接近；在類內(nèi)聚集性較小時，兩隨機效應(yīng)模型95%CI覆蓋率均較好，但隨類內(nèi)聚集性增大，Log-binom ial的隨機效應(yīng)模型95%CI覆蓋率變動較大，且類別數(shù)增加時覆蓋率普遍較低。結(jié)論基于穩(wěn)健Poisson的隨機效應(yīng)模型收斂率高，參數(shù)估計結(jié)果與真實值接近，對各固定效應(yīng)參數(shù)估計的95%CI的覆蓋率高且穩(wěn)定，優(yōu)于Log-binomial隨機效應(yīng)模型。

非獨立 穩(wěn)健Poisson回歸 Log-binom ial模型 隨機效應(yīng) 患病率比

流行病學封閉性隊列研究和橫斷面研究中，OR值被廣泛應(yīng)用于RR/PR的近似估計，但研究結(jié)局出現(xiàn)頻率較高（如大于10%）時，會嚴重高估暴露對結(jié)局的影響，這一問題早已引起統(tǒng)計學者注意，因此直接估計RR/PR的模型被提出，如Log-binomial模型和穩(wěn)健Poisson回歸模型，其參數(shù)解釋更為合理［1］。但基于穩(wěn)健Poisson回歸和Log-binom ial模型的相關(guān)理論和應(yīng)用研究尚顯不足，因此本研究通過計算機模擬，探討穩(wěn)健Poisson和Log-binomial的隨機效應(yīng)模型參數(shù)估計的有效性及應(yīng)用中的相關(guān)問題。

模型原理與方法

本研究以方差分量模型為基礎(chǔ)，探討基于穩(wěn)健Poisson回歸及Log-binomial的隨機效應(yīng)模型。

1.穩(wěn)健Poisson隨機效應(yīng)模型

設(shè)yki和Xki＝（x1ki，x2ki，…，xpki）T分別為第k（k＝1，2，…，K）類內(nèi)第i（i＝1，2，…，nk）個個體的二分類結(jié)局變量和P×1維解釋變量向量，uk為第k類的隨機效應(yīng)，則構(gòu)建隨機效應(yīng)模型：

式中，pki＝Pr（yki＝1｜Xki），連接函數(shù)采用log函數(shù)，βp為固定效應(yīng)參數(shù)，反映固定效應(yīng)xp對結(jié)局概率對數(shù)的影響，且RR（PR）＝exp（βp）；假定隨機效應(yīng)uk來自正態(tài)分布總體，即隨機效應(yīng)方差的大小反映結(jié)局變量在類內(nèi)的聚集程度。

當誤差分布指定為Poisson分布時，由于Poisson分布的方差無邊界（方差等于均數(shù)），當應(yīng)用到二項分布資料時，易出現(xiàn)過度離散問題，即高估參數(shù)的方差，導致較寬的置信區(qū)間。處理高估方差的一個經(jīng)典方法是引入Huber的穩(wěn)健“三明治”方差：

其中，

從而構(gòu)建穩(wěn)健Poisson隨機效應(yīng)模型。

2.Log-binomial隨機效應(yīng)模型

Log-binom ial隨機效應(yīng)模型同式（1），但誤差分布指定為二項分布。應(yīng)用于獨立數(shù)據(jù)時，基本的Log-binom ial模型最大似然估計時需加一個限制條件，對所有個體，令

最大似然估計的解應(yīng)該在此條件限制的參數(shù)空間中。Log-binomial模型在參數(shù)估計過程中，最大似然估計的解如果落在有限制的參數(shù)空間邊界上，易出現(xiàn)不收斂的情形，模型中包含連續(xù)型協(xié)變量時尤甚。實際工作中經(jīng)驗性地設(shè)置截距β0＝-4可能效果較好［2］，Deddens建議也可先對原始數(shù)據(jù)集調(diào)整擴充后再擬合Log-binomial模型，即可得到參數(shù)的近似最大似然估計（maximum-likelihood estimation，簡稱MLE）值，稱為COPY算法［3］，獨立數(shù)據(jù)情況下，其收斂性和近似MLE的唯一性已有嚴格的理論證明［4］。

目前穩(wěn)健Poisson和Log-binomial的隨機效應(yīng)模型都可在SAS中采用Proc Glimmix過程實現(xiàn)，連接函數(shù)選用log連接，誤差分布分別指定為Poisson分布和二項分布，并用random語句指定截距項以擬合方差分量模型。

模擬研究

1.模擬方法及參數(shù)設(shè)置

模擬研究生成具有類內(nèi)相關(guān)的數(shù)據(jù)集，按log-binomial隨機效應(yīng)模型產(chǎn)生。數(shù)據(jù)集中包含兩個自變量X1和X2，假定X1為二分類自變量（1＝暴露；0＝非暴露），X2為連續(xù)型協(xié)變量；Y為二分類因變量（1＝結(jié)局發(fā)生；0＝結(jié)局未發(fā)生）。每個數(shù)據(jù)集共有N類，每類包括Nk個個體，因變量Y具有類內(nèi)相關(guān)，類間獨立的特性。

（1）類別數(shù)N和類內(nèi)個體數(shù)Nk的設(shè)置

假定每個數(shù)據(jù)集的平均樣本量為500，模擬類別數(shù)N為20和50兩種情況。類別數(shù)N＝20時，類內(nèi)個體數(shù)Nk由均勻分布U（0，50）隨機產(chǎn)生，即Nk的均數(shù)為25；類別數(shù)N＝50時，類內(nèi)個體數(shù)Nk由均勻分布U（0，20）隨機產(chǎn)生，即Nk的均數(shù)為10［5-6］。

（2）隨機效應(yīng)及參數(shù)設(shè)置

假設(shè)隨機效應(yīng)uk來自正態(tài)分布，均數(shù)為0，方差為0.1或0.2兩種情況，分別相當于類內(nèi)相關(guān)系數(shù)在0.01和0.15之間，其大小依賴于暴露因素和協(xié)變量的效應(yīng)。研究顯示這種大小在實際資料中常見。

（3）固定效應(yīng)及參數(shù)設(shè)置

固定效應(yīng)x1ki為二分類變量，假定來自二項分布，即x1ki～B（p）；x2ki為連續(xù)變量，假定來自均勻分布，且依賴于x1ki，即x2ki/x1ki～U（-6＋2x1ki，2＋2x1ki）。考慮基線患病率為0.3和0.15兩種情況，暴露因素X1的PR值分別為1.5和2.0，暴露率分別為0.2和0.5，協(xié)變量X2的回歸系數(shù)分別為0.18，0.10，0.36和0.20［7］。共組成8種參數(shù)情況，如表1。

表1 模擬研究固定效應(yīng)參數(shù)設(shè)置情況

（4）因變量Y的設(shè)置

兩分類結(jié)局變量Y（1/0）來自n為1的二項分布，概率Pij為

模擬過程中，如果固定效應(yīng)和隨機效應(yīng)各組合的結(jié)果使Pki≥1，舍棄并模擬產(chǎn)生新的值，直到Pki＜1。

在上述各參數(shù)設(shè)置情況下，本研究共模擬2（類別數(shù)為20和50）×2（隨機效應(yīng)方差0.1和0.2）×8（固定效應(yīng)參數(shù)設(shè)置情況1-8）＝32種情況。每種情況下，改變隨機數(shù)的種子共模擬不相關(guān)的1000個數(shù)據(jù)集。所有模擬數(shù)據(jù)集均在SAS9.2軟件中產(chǎn)生，并擬合穩(wěn)健Poisson隨機效應(yīng)模型和Log-binom ial隨機效應(yīng)模型。

2.評價指標

（1）收斂率

對每種模擬情況，考察基于穩(wěn)健Poisson、Log-binomial的兩類隨機效應(yīng)模型的收斂情況并計算收斂率。

（2）參數(shù)估計值的均值和95%CI覆蓋率

在每種模擬情況下，根據(jù)每種模型的1000個參數(shù)估計結(jié)果，分別計算回歸系數(shù)β0，β1和β2估計值的均值及95%CI對參數(shù)真值的覆蓋率。

3.結(jié)果

（1）收斂率

表2顯示32種模擬情況下兩類隨機效應(yīng)模型參數(shù)估計的收斂率（%）。32種模擬條件下穩(wěn)健Poisson隨機效應(yīng)模型的收斂率波動范圍為81.4%～ 100.0%，而Log-binomial隨機效應(yīng)模型的收斂率波動范圍為4.7%～98.6%，波動較大。類別數(shù)為20時穩(wěn)健Poisson隨機效應(yīng)模型和Log-binomial隨機效應(yīng)模型的最低收斂率分別為81.4%和5.9%；類別數(shù)為50時穩(wěn)健Poisson隨機效應(yīng)模型和Log-binomial隨機效應(yīng)模型的最低收斂率分別為91.9%和4.7%。

表2 32種模擬情況下兩類隨機效應(yīng)模型參數(shù)估計的收斂率（%）

（2）參數(shù)估計值的均值與95%CI覆蓋率

表3和表4顯示了2類模型在不同隨機效應(yīng)方差水平、類別數(shù)及固定效應(yīng)參數(shù)條件下分別進行1000次參數(shù)估計的β0、β1、β2的均值和95%CI覆蓋率。兩模型對各設(shè)定條件下的參數(shù)估計結(jié)果與真值均較為接近，受類別數(shù)或隨機效應(yīng)方差的影響不大。

穩(wěn)健Poisson隨機效應(yīng)模型對各參數(shù)的覆蓋率均在80%以上，和Log-binomial隨機效應(yīng)模型相比，95%CI覆蓋率變異較??；在隨機效應(yīng)方差較小（0.1）時，穩(wěn)健Poisson和Log-binomial的隨機效應(yīng)模型參數(shù)估計的95%CI覆蓋率較接近，受類別數(shù)目影響較??；當隨機效應(yīng)方差增大和類別數(shù)同時增大時，Log-binomial隨機效應(yīng)模型估計的各參數(shù)95%CI覆蓋率降低，且變異較大。

表3 類別數(shù)為20時，兩種模型參數(shù)估計均值與真值的比較以及95%CI的覆蓋率

*：模型1和模型2分別是穩(wěn)健Poisson-隨機效應(yīng)和Log-binom ial-隨機效應(yīng)（COPY算法）

表4 類別數(shù)為50時，兩種模型參數(shù)估計均值與真值的比較以及95%CI的覆蓋率

討 論

隊列研究或橫斷面研究資料中，穩(wěn)健Poisson回歸和log-binom ial模型常用于估計RR/PR，當用于獨立數(shù)據(jù)時，Deddens［8］顯示中等樣本量情況下，穩(wěn)健Poission回歸對PR估計的偏倚小于Log-binomial模型，但Log-binomial模型所估計的患病概率不會大于1，且可以使用似然比檢驗，相較Wald檢驗可能更好。但運用穩(wěn)健Poisson模型幾乎很少有收斂問題，又可直接對原始數(shù)據(jù)分析，省去了數(shù)據(jù)操作的過程而更容易使用，因此學者建議［9-10］，log-binomial不收斂時可首先考慮穩(wěn)健Poisson回歸。

另一方面，解決非獨立數(shù)據(jù)的問題，很自然地思路是將其推廣到多水平模型，因此本研究主要探討Logbinom ial模型和穩(wěn)健Poisson模型的隨機效應(yīng)模型，采用計算機模擬實驗考察不同組合情況下模型參數(shù)估計的有效性和準確性。

本研究中的模擬實驗共設(shè)置了32種參數(shù)情況，分別反映實際資料中的少數(shù)類別數(shù)和中等類別數(shù)，小強度和中等強度的類內(nèi)聚集，模擬結(jié)果顯示，各類參數(shù)設(shè)置情況下，兩類隨機效應(yīng)模型的參數(shù)估計值與真值差別均不大，估計RR/PR有較高的準確性，且穩(wěn)健Poisson隨機效應(yīng)模型的覆蓋率較高。但隨著類內(nèi)聚集性增大，Log-binom ial的隨機效應(yīng)模型95%CI覆蓋率變動較大，且類別數(shù)增加時覆蓋率普遍較低，表明穩(wěn)健Poisson隨機效應(yīng)模型的穩(wěn)健性較好。此外Log-binom ial隨機效應(yīng)模型的收斂效果劣于穩(wěn)健Poisson隨機效應(yīng)模型，這與獨立數(shù)據(jù)時情況一致，在Log-binomial模型基礎(chǔ)上增加隨機效應(yīng)后可能導致更多參數(shù)估計的不確定性。

盡管本研究不能窮盡各種可能的實際情況，但模擬條件下結(jié)果顯示，和Log-binomial隨機效應(yīng)模型相比，運用穩(wěn)健Poisson隨機效應(yīng)模型幾乎很少有收斂問題，又可直接對原始數(shù)據(jù)分析，省去了COPY算法中數(shù)據(jù)操作過程而更容易使用；參數(shù)估計值與真值的絕對偏差均較小，對各固定效應(yīng)參數(shù)估計的95%CI的覆蓋率高且穩(wěn)定，故在結(jié)局頻率較高的非獨立數(shù)據(jù)中估計RR/PR時推薦使用。不過Kauermann＆Carroll在實驗中發(fā)現(xiàn)三明治方差的估計效率較基于模型的方差估計方法的差，尤其是當數(shù)據(jù)來源于小樣本的實驗單位而且不同實驗單位間協(xié)變量不同時，主要表現(xiàn)為參數(shù)估計的置信區(qū)間覆蓋率較低，此時建議對置信區(qū)間進行校正［12］。

1.Str?mberg U.Prevalence odds ratio v prevalence ratio.Occup Environ Med，1994，51（2）：143-144.

2.Deddens JA，Petersen MR，Lei X.Estimation of prevalence ratioswhen PROC GENMOD does not converge.Proceedings of the 28th Annual SASUsers Group International Conference，Seattle，Washington，2003：270-28.

3.Deddens JA，Petersen MR.Approaches for estimating prevalence ratios. Occup Environ Med，2008，65（7）：501-506.

4.Savu A，Liu Q，Yasui Y.Estimation of relative risk and prevalence ratio.Stat Med，2010，29（22）：2269-2281.

5.Yelland LN，Salter AB，Ryan P.Relative Risk Estimation in Cluster Randomized Trials：A Comparison of Generalized Estimating Equation Methods.International Jof BIOST，2011，7（1）：1-26.

6.Yu B，Wang Z.Estimating relative risks for common outcome using PROC NLP.Comput Methods Programs Biomed，2008，90（2）：179-186.

7.Zou GY，Donner A.Extension of themodified Poisson regressionmodel to prospective studies w ith correlated binary data.Stat Methods Med Res，2013，22（6）：661-670.

8.Petersen MR，Deddens JA.A comparison of twomethods for estimating prevalence ratios.BMC Med Res Methodol，2008，8：9.

9.Lee J，Tan C S，Chia K S.A practical guide formultivariate analysis of dichotomous outcomes.Ann Acad Med Singapore，2009，38（8）：714-719.

10.Spiegelman D，Hertzmark E.Easy SAS calculations for risk or prevalence ratios and differences.Am JEpidem iol，2005，162（3）：199-200.

11.Goran K，Raymond JC.A note on the efficiency of sandw ich covariancematrix estimation.JAm Stat Assoc，2001，12：1387-1396.

（責任編輯：郭海強）

A Com parative Study between Robust Poisson and Log-binom ial M odel w ith Random Effects M odels to Estimate RR/PR in Non-independent Data

Zhou Shudong，Gao Yanhui，Li Lixia，et al（DepartmentofEpidemiologyandBiostatistics，SchoolofPublicHealth，Guangdong PharmaceuticalUniversity，GuangdongKeyLaboratoryofMolecularEpidemiology（510310），Guangzhou）

ObjectiveTo compare the parameter estimation results based on robust Poisson and Log-binomialmodels w ith random effects in the hierarchical and uncommon outcome datasets.MethodsThe simulation datasets were generated by computer，including 32 kinds of parameters conditions，respectively，reflecting the small and the medium number of categories，low and medium aggregation w ithin the classes，0.15 and 0.3 of the baseline prevalence，1.5 and 2.0 of the PR value，20%and 50%of the cases’exposure rate.The convergence rates，parameter estimation results and the 95%CIcoverage parameters in all the simulations were compared.ResultsThe convergence effects of Log-binom ialmodel w ith random effects were worse than the robust Poisson.The estimation results of two random effectsmodel were both closer to the true value.Two random effects models’95%CIcoveragewere good when the classes’aggregation was low，butw ith the addition in intra-class clustering，95% CI coverage of Log-binom ial random effectsmodelwas unstable，and the coverage was generally low when the number of categories increased.ConclusionThe convergence rates of robust Poisson random effectsmodel were good，the parameter estimation resultswere close to the true value，and the 95%CIcoverageswere high and stable，superior to Log-binom ial random effects model.

Non-independent；Robust Poisson regression；Log-binomialmodel；Random effects；Prevalence ratio

廣東省自然科學基金（項目號：10151022401000018）

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

2.Department of Epidem iology and Community Medicine，University of Ottawa.

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