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

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

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

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

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

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

模型原理與方法

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

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

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

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

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

其中，

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

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

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

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

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

模擬研究

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

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

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

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

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

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

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

固定效應(yīng)x1ki為二分類(lèi)變量，假定來(lái)自二項(xiàng)分布，即x1ki～B（p）；x2ki為連續(xù)變量，假定來(lái)自均勻分布，且依賴于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è)置

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

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

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

2.評(píng)價(jià)指標(biāo)

（1）收斂率

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

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

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

3.結(jié)果

（1）收斂率

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

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

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

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

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

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

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

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

討 論

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

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

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

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

（責(zé)任編輯：郭海強(qiáng)）

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

廣東省自然科學(xué)基金（項(xiàng)目號(hào)：10151022401000018）

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

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

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