重慶醫(yī)科大學(xué)公共衛(wèi)生與管理學(xué)院衛(wèi)生統(tǒng)計學(xué)教研室（400016） 文 雯 文小焱 胡 珊 彭 斌

貝葉斯層次模型在嵌套結(jié)構(gòu)調(diào)查數(shù)據(jù)中的應(yīng)用研究*

重慶醫(yī)科大學(xué)公共衛(wèi)生與管理學(xué)院衛(wèi)生統(tǒng)計學(xué)教研室（400016） 文 雯 文小焱 胡 珊 彭 斌△

目的針對分層抽樣流行病調(diào)查數(shù)據(jù)的結(jié)構(gòu)特點，構(gòu)建兩種基于分層嵌套思想的貝葉斯層次模型，并探討其優(yōu)缺點。方法以貝葉斯層次模型為基礎(chǔ)，利用嵌套結(jié)構(gòu)中的層級關(guān)系構(gòu)建模型，其中，模型一以嵌套層效應(yīng)分解為特點構(gòu)建，模型二以嵌套層效應(yīng)逐級傳遞為特點構(gòu)建。以重慶市出生缺陷調(diào)查數(shù)據(jù)為例，采用OpenBUGS軟件進行模型擬合及分析。結(jié)果以偏差信息準則（deviance information criterion，DIC）作為擬合優(yōu)度評價，模型一和模型二的DIC值分別為101.8和101.6，大致相等；敏感性分析顯示，在總體率的超參數(shù)μ設(shè)置不同先驗信息下，模型一和模型二對總效應(yīng)估計的變異性分別為（用標準差度量，10－4）：后驗均數(shù)1.191和27.546；后驗中位數(shù)1.038和7.617，模型一的變異性比模型二小。結(jié)論模型一和模型二均可用于嵌套結(jié)構(gòu)的調(diào)查數(shù)據(jù)建模分析及預(yù)測，擬合效果相當；但模型一比模型二受先驗信息影響小，穩(wěn)健性更好，更適合先驗信息欠缺時的數(shù)據(jù)分析。

嵌套結(jié)構(gòu)數(shù)據(jù) 分層抽樣 貝葉斯層次模型 OpenBUGS

隨著馬爾科夫鏈蒙特卡洛（Markov chain Monte carlo，MCMC）方法的不斷成熟以及計算機軟件的快速發(fā)展，貝葉斯統(tǒng)計方法被越來越多地用于科學(xué)研究和數(shù)據(jù)處理［1］。與經(jīng)典統(tǒng)計方法不同，貝葉斯統(tǒng)計方法將未知參數(shù)看作隨機變量，用一個概率分布去描述，稱為先驗分布。傳統(tǒng)的統(tǒng)計方法應(yīng)用樣本信息和總體信息進行統(tǒng)計推斷，而貝葉斯方法則利用樣本信息、總體信息及先驗信息進行統(tǒng)計推斷［2］。當有先驗信息可以利用時，采用貝葉斯方法可以得到更好的估計結(jié)果［3］。貝葉斯層次模型是在貝葉斯理論的基礎(chǔ)上通過構(gòu)造多層先驗分布將模型以多層次的形式表現(xiàn)出來，即當參數(shù)的先驗分布含有超參數(shù)時，對超參數(shù)再給出一個先驗，稱為超先驗。由先驗和超先驗共同構(gòu)建的層次先驗作為模型的新先驗，并采用貝葉斯統(tǒng)計推斷方法對模型進行參數(shù)估計。這種方式構(gòu)造的先驗性質(zhì)穩(wěn)健，是社會科學(xué)研究領(lǐng)域廣泛使用的建模方法之一［4－5］。嵌套結(jié)構(gòu)是利用分層的思想將研究對象逐步劃分成若干層，每一次劃分都是建立在上一次分層的基礎(chǔ)上，形成一層套一層類似蜂巢的結(jié)構(gòu)。本文則是將這種思想應(yīng)用于分層隨機抽樣調(diào)查數(shù)據(jù)的建模分析，利用多層貝葉斯統(tǒng)計方法，構(gòu)建適應(yīng)嵌套結(jié)構(gòu)數(shù)據(jù)的層次貝葉斯模型，并探討模型在該類型數(shù)據(jù)中的應(yīng)用。

貝葉斯層次模型構(gòu)建

1.抽樣調(diào)查研究數(shù)據(jù)的嵌套結(jié)構(gòu)特點

在抽樣調(diào)查研究中，經(jīng)常會選擇分層隨機抽樣方法，即將調(diào)查總體（N）按其屬性特征劃分成若干層，然后在每個層中進行隨機抽樣。有時為了樣本的代表性和實際需要，可進行多次分層，即在上一級分層基礎(chǔ)上進行再分層，形成一層套一層的嵌套結(jié)構(gòu)［6］。

2.構(gòu)建貝葉斯層次模型

根據(jù)抽樣調(diào)查研究數(shù)據(jù)的嵌套結(jié)構(gòu)特點，對嵌套結(jié)構(gòu)的每一層構(gòu)建模型。下面以三層嵌套結(jié)構(gòu)為例，闡述構(gòu)建估計總體陽性率的貝葉斯層次模型的思路及具體方法。

模型一：基于廣義線性模型的各嵌套層效應(yīng)分解加成模型，簡稱效應(yīng)分解加成模型。該模型以分層效應(yīng)作為各嵌套層效應(yīng)差別的主要原因，各嵌套層某分層單元的總效應(yīng)等于其上一層總效應(yīng)與該層分層效應(yīng)之和。三層嵌套結(jié)構(gòu)下各層所對應(yīng)的效應(yīng)模型分別為：

（1）～（4）式中，θ0表示總體陽性率，μ0為總體陽性率對應(yīng)的平均效應(yīng)；αi表示第一層第i分層單元的效應(yīng)，且表示嵌套于第一層第i分層單元下第二層第j分層單元的效應(yīng)，且γi（j（k））表示嵌套于第一層第i分層單元第二層第j分層單元下第三層第k分層單元的效應(yīng)，且＝0。θi、θi（j）、θi（j（k））分別為相應(yīng)嵌套層的總體陽性率。實際中，只需要擬合（4）式的模型即可估計出各層的參數(shù)。

模型參數(shù)的先驗信息設(shè)定：μ0作為總體陽性率的平均效應(yīng)，本文假設(shè)其服從均值為μ的正態(tài)分布，即μ0～norm（u，δ2）；α、β、γ作為各嵌套層的效應(yīng)，可假設(shè)其服從均值為0的正態(tài)分布。

模型超參數(shù)設(shè)定：μ一般設(shè)定為正態(tài)分布或均勻分布，其均值可根據(jù)經(jīng)驗或文獻資料確定；方差先驗通常設(shè)置為逆伽瑪分布（inverse gamma distribution）［7］，即方差的倒數(shù)（精度參數(shù)τ）服從伽瑪分布τ～G（r，mu）。當沒有方差先驗信息可以利用時，通常需要將方差成分設(shè)置得足夠大，例如106，精度系數(shù)趨于0，此時的（逆）伽瑪分布就等價于一個無信息先驗分布。

模型二：基于廣義線性模型的各嵌套層效應(yīng)逐級傳遞模型，簡稱效應(yīng)逐級傳遞模型。該模型以各嵌套層效應(yīng)僅受其上一層效應(yīng)影響為出發(fā)點，各嵌套層間的嵌套關(guān)系通過層次先驗來反映。前述三層嵌套結(jié)構(gòu)對應(yīng)的效應(yīng)模型為：

（5）～（8）式中，θ0、θi、θi（j）、θi（j（k））、μ0含義與模型一相同；μi為第一層第i分層單元的總效應(yīng)；μi（j）為嵌套于第一層第i分層單元下第二層第j分層單元的總效應(yīng)；μi（j（k））為嵌套于第一層第i分層單元第二層第j單元下第三層第k單元的總效應(yīng)。

模型參數(shù)的先驗信息設(shè)定：根據(jù)該模型構(gòu)建思路，可假設(shè)每一嵌套層的效應(yīng)服從以其上一級的效應(yīng)為均值的正態(tài)分布［8］。即

模型超參數(shù)設(shè)定：與模型一類似，可根據(jù)先驗信息對超參數(shù)μ設(shè)定一個正態(tài)分布或均勻分布，方差均采用無信息先驗分布。

實例分析

1.資料來源

數(shù)據(jù)來源于2010年重慶市0～5歲兒童出生缺陷基線調(diào)查。該研究采用多階段分層隨機抽樣，第一階段根據(jù)重慶市“一圈兩翼”發(fā)展戰(zhàn)略將全部40個區(qū)縣劃分為“一小時經(jīng)濟圈”、“渝東北翼”和“渝東南翼”三個經(jīng)濟區(qū)，第二階段以區(qū)縣為單位，按比例從三個經(jīng)濟區(qū)總共隨機抽取16個區(qū)縣作為調(diào)查點。

2.建立兩層嵌套結(jié)構(gòu)數(shù)據(jù)的貝葉斯層次模型

重慶市0～5歲兒童出生缺陷基線調(diào)查數(shù)據(jù)可視為兩層嵌套結(jié)構(gòu)的數(shù)據(jù)，即將重慶市（θ0）按照“一圈兩翼”劃分的三個經(jīng)濟區(qū)（θi）作為第一層嵌套，每個經(jīng)濟區(qū)按比例抽取出的區(qū)縣（θi（j））作為第二層嵌套。設(shè)每個調(diào)查區(qū)縣0～5歲兒童調(diào)查人數(shù)為ni（j），出生缺陷人數(shù)為xi（j），則可假設(shè)其服從二項分布，即

根據(jù)前面介紹的貝葉斯層次模型構(gòu)建方法，對應(yīng)的兩類模型分別為：

模型一（效應(yīng)分解加成模型）：

模型二（效應(yīng)逐級傳遞模型）：

通過文獻查閱，我國每年新增的先天殘疾和智力缺陷兒童總數(shù)高達120萬，約占每年出生人口的4%～6%［9－10］，鑒于重慶市歷年來通過醫(yī)院監(jiān)測到的圍產(chǎn)兒出生缺陷發(fā)生率均在1%左右［11－12］，而以醫(yī)院為基礎(chǔ)的監(jiān)測方法由于研究對象受限以及篩查工作的不完善，統(tǒng)計出的出生缺陷發(fā)生率通常要比實際情況低。因此，本文認為總體出生缺陷發(fā)生率很可能介于1%～6%之間，其對應(yīng)的u值范圍為－4.595～－2.752（根據(jù)logit（θ）＝u推算），故可假設(shè)μ服從－4.595～－2.752之間的均勻分布，即μ～uniform（－4.595，－2.752）。方差先驗采用逆伽瑪分布，比較常用的設(shè)置有r＝mu＝0.001［13－14］或r＝mu＝1、r＝mu＝0.1、r＝mu＝0.01［15］或r＝0.5，mu＝0.005［16］，本文對以上設(shè)置進行了模擬研究，發(fā)現(xiàn)模型一和模型二在r＝mu＝0.01時相對較優(yōu)。因此，本研究的方差先驗采用τ～G（0.01，0.01）。

3.敏感性分析

敏感性分析是為了評估不同先驗分布對模型后驗的影響，即考察模型的穩(wěn)健性（robustness）。當有可利用的先驗信息時，敏感性分析主要考察先驗分布的設(shè)置對模型后驗的影響，當無信息先驗被利用時，則主要考察不同無信息先驗的選擇對模型后驗的影響［17］。鑒于本文所構(gòu)建的模型中參數(shù)較多，且兩種先驗均有利用，因此本文選取對模型影響較大的超參數(shù)μ以及采用無信息先驗的方差進行敏感性分析。

超參數(shù)μ在模型一和模型二中均表示總體效應(yīng)μ0的先驗均值，本文通過考查總體出生缺陷率θ0的先驗信息將超參數(shù)μ設(shè)置如下五組先驗范圍：（1）通過經(jīng)驗和歷史數(shù)據(jù)估計重慶市θ0的可能范圍為1%～6%；（2）將（1）中設(shè)置的范圍擴大到0.1%～20%，即理論上大多數(shù)地區(qū)出生缺陷發(fā)生率均在此范圍；（3）設(shè)置一個與（1）無交叉，取值較小且包含于（2）中的范圍：0.1%～1%；（4）設(shè)置一個與（1）無交叉，取值較大且包含于（2）中的范圍：10%～20%；（5）設(shè)置一個無信息先驗，即0%～100%。上述范圍對應(yīng)的超參數(shù)μ的先驗分別為（1）μ～unif（－4.595，－2.752）（θ0介于1%～6%的均勻分布）；（2）μ～unif（－6.907，－1.386）（θ0介于0.1%～20%的均勻分布）；（3）μ～ unif（－6.907，－4.595）（θ0介于0.1%～1%的均勻分布）；（4）μ～unif（－2.197，－1.386）（θ0介于10%～20%的均勻分布）；（5）μ～unif（－10，10）（θ0介于0%～100%的均勻分布，即無信息均勻分布）。

將各種先驗設(shè)置下參數(shù)的估計值與無信息先驗下參數(shù)的估計值進行比較，以考察模型對超參數(shù)μ先驗設(shè)置的穩(wěn)健性。

4.模型擬合及診斷

所有模型均在OpenBUGS 3.2.2中運行，參數(shù)估計采用MCMC。平行模擬初始值不同的兩條鏈，以克服不同初始值對模型的影響；采用Gelman-Rubin統(tǒng)計量目測各參數(shù)的迭代時序圖，以評判模型是否收斂；模型穩(wěn)定后迭代20000次用于參數(shù)估計和DIC的計算［17］。

5.結(jié)果及比較

（1）出生缺陷率的估計

以后驗中位數(shù)作為出生缺陷率的點估計值，模型一和模型二的估計值分別為39.75‰和39.00‰，估計結(jié)果非常接近。模型一與模型二的DIC值分別為101.8和101.6，基本相等，說明兩個模型擬合效果相當（表1）。此外，按照傳統(tǒng)統(tǒng)計方法計算的出生缺陷率為37.92‰，模型一與模型二對全市出生缺陷率的估計值略高于傳統(tǒng)統(tǒng)計方法的估計結(jié)果。

表1 出生缺陷發(fā)生率估計值（‰）及DIC值比較結(jié)果

（2）超參數(shù)μ不同先驗的敏感性分析

由于模型一和模型二中參數(shù)較多，為了便于對比，本文僅選取了主要參數(shù)θ0的估計值（包括后驗均值和后驗中位數(shù)）進行呈現(xiàn)（表2）。與無信息先驗相比，當先驗設(shè)置為0.1%～1%時，模型一和模型二θ0估計值均較小，后驗均數(shù)及后驗中位數(shù)的差值（10－4）分別為－0.6和－0.5（模型一）、－25.2和－7.1（模型二）；當先驗設(shè)置為10%～20%時，模型一和模型二θ0估計值均較大，差值（10－4）分別為2.4和2.1（模型一）、43.2和13.5（模型二）。從標準差這一角度看，模型一各先驗設(shè)置對應(yīng)參數(shù)θ0的后驗均值和后驗中位數(shù)標準差（10－4）分別為1.191和1.038，差別不大；模型二分別為27.546和7.617，差別較大，但后驗中位數(shù)較后驗均值更穩(wěn)定?？梢钥闯觯瑑蓚€模型后驗估計值均在一定程度上受到先驗信息的影響，但模型一在各種先驗下的后驗均值與后驗中位數(shù)基本相等，且變異較小，模型表現(xiàn)穩(wěn)健。模型二則受到不同先驗設(shè)置的影響，參數(shù)估計值變異較大，欠穩(wěn)定，尤其是后驗均值受先驗信息影響更明顯。

表2 超參數(shù)μ不同先驗的敏感性分析結(jié)果

（3）方差先驗的敏感性分析

對模型一和模型二采用前述5種方差先驗進行擬合度的敏感性分析。結(jié)果顯示（表3）五種方差先驗信息下，模型一或模型二的DIC值相差均不明顯。因此，模型一和模型二對方差先驗表現(xiàn)穩(wěn)健。考慮到先驗IG（0.01，0.01）在兩個模型中的DIC值均較小，故選擇該先驗作為模型擬合時的方差先驗。

表3 五種方差先驗下的模型DIC值

討 論

在抽樣調(diào)查研究中，采用分層隨機抽樣方法獲取的數(shù)據(jù)具有一層套一層的嵌套結(jié)構(gòu)特點。本研究通過分析數(shù)據(jù)的嵌套結(jié)構(gòu)特點，構(gòu)建了兩個基于嵌套結(jié)構(gòu)的貝葉斯層次模型：模型一（效應(yīng)分解加成模型）和模型二（效應(yīng)逐級傳遞模型），并利用重慶市出生缺陷調(diào)查數(shù)據(jù)，分析比較兩個模型的優(yōu)缺點。研究結(jié)果顯示，模型一和模型二在擬合優(yōu)度方面相差不大，模型二的DIC值略小于模型一，提示兩個模型均能較好地擬合數(shù)據(jù)；穩(wěn)健性方面，模型一具有較好的穩(wěn)鍵性，模型二穩(wěn)健性較差，尤其是后驗均值極易受μ先驗分布影響。另外，從模型的潛力看，模型一將各嵌套層的分層效應(yīng)從該層總效應(yīng)中分解出來，作為各嵌套層不同分層單元陽性率差異的主要原因，有利于探討模型的內(nèi)在結(jié)構(gòu)及各層的效應(yīng)貢獻。在實例分析中，分層效應(yīng)是由于區(qū)域間的地理位置差異造成，若是將地理位置細化（空間效應(yīng)），或是考慮更多相關(guān)影響因素（如時間效應(yīng)、經(jīng)濟效應(yīng)等），將會提高模型預(yù)測的準確性［17］。模型二則是在各層的總效應(yīng)水平上，通過構(gòu)造各參數(shù)先驗分布的嵌套關(guān)系進行模型構(gòu)建，下一嵌套層的效應(yīng)僅依賴于其上一層的效應(yīng)，且通過逐級隨機傳遞的方式進行影響，其獨特的層次先驗結(jié)構(gòu)使模型得以簡化的同時，失去了繼續(xù)優(yōu)化的可能。

關(guān)于超參數(shù)μ先驗分布的選擇，如果可以通過經(jīng)驗和文獻資料收集到關(guān)于μ0充分的先驗信息，則可利用這些信息確定μ值或其先驗分布；若缺乏μ0的相關(guān)信息，則可根據(jù)研究具體情況選擇一個合理范圍的均勻分布或無信息先驗。本文不建議在μ0無先驗信息可利用的情況下對μ直接賦值0，即μ0～norm（0，δ2）這相當于將陽性率θ0默認設(shè)置為50%。由于模型一和模型二均受μ0先驗的影響，尤其是模型二，這將會導(dǎo)致模型中主要參數(shù)估計過高。

綜上所述，當有合理的先驗信息（如合理的μ值范圍）可利用時，模型一與模型二均能應(yīng)用于嵌套結(jié)構(gòu)數(shù)據(jù)建模分析及預(yù)測；但當先驗信息不明確或欠缺時，考慮到模型二對先驗信息過于敏感，選擇模型一會更恰當，且模型一可通過細化分層效應(yīng)以達到提高模型預(yù)測準確性的目的。

1.Smith AFM，Roberts GO.Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods.Journal of the Royal Statistical Society（Series B），1993，55（1）：3-23.

2.張堯庭，陳漢峰.貝葉斯統(tǒng)計推斷.科學(xué)出版社，1991.

3.Tu XM，Kowalski J，Jia G.Bayesian analysis of prevalence with covariates using simulation based techniques：applications to HIV screening. Statistics in Medicine，1999，18（22）：3059-3073.

4.王哲.多層貝葉斯方法在消費者行為中的應(yīng)用研究.南京：南京航空航天大學(xué)，2012.

5.Lindley DV，Smith AFM.Bayes estimates for the linear model.Journal of the Royal Statistical Society（Series B），1972，34（1）：1-41.

6.王建華.流行病學(xué).北京：人民衛(wèi)生出版社，2008.

7.Daniels MJ.A prior for the variance in hierarchical models.Canadian Journal of Statistics，1999，27（3）：567-578.

8.王顯紅.日本血吸蟲病貝葉斯時空模型的建立.北京：中國疾病預(yù)防控制中心，2007.

9.毛萌，朱軍.出生缺陷檢測研究現(xiàn)狀.實驗兒科臨床雜志，2009，24（11）：801-803.

10.李?；?，田宏，陳艷玲，等.遼寧省2011年度出生缺陷監(jiān)測數(shù)據(jù).中國衛(wèi)生統(tǒng)計，2012，29（3）：410-411.

11.王繼林，樊欣.重慶市1996－2001年出生缺陷監(jiān)測結(jié)果分析.中國兒童保健雜志，2003，11（3）：207-209.

12.張高東，周文正，周曉軍，等.重慶市2004－2010年出生缺陷發(fā)生情況分析.華中科技大學(xué)學(xué)報，2012，41（6）：759-762.

13.Su Z，Peterman RM，Haeseker SL.Spatial hierarchical Bayesian models for stock-recruitment analysis of pink salmon（Oncorhynchusgorbuscha）.Canadian Journal of Fisheries and Aquatic Sciences，2004，61（12）：2471-2486.

14.Kazembe LN，Namangale JJ.A Bayesian multinomial model to analyse spatial patterns of childhood co-morbidity in Malawi.European journal of Epidemiology，2007，22（8）：545-556.

15.Su Z，Adkison MD，Van Alen BW.A hierarchical Bayesian model for estimating historical salmon escapement and escapement timing.Canadian Journal of Fisheries and Aquatic Sciences，2001，58（8）：1648-1662.

16.Kelsall JE，Diggle PJ.Spatial variation in risk of disease：a nonparametric binary regression approach.Journal of the Royal Statistical Society（Series C），1998，47（4）：559-573.

17.Ntzoufras I.Bayesian modeling using WinBUGS.Wiley.com，2011.

18.張俊輝，馮子健，楊超，等.基于層次貝葉斯時空模型的空間多尺度聯(lián)合分析模型的構(gòu)建及應(yīng)用研究.中國衛(wèi)生統(tǒng)計，2013，30（2）：199-202.

（責(zé)任編輯：劉壯）

Application of Hierarchical Bayesian Model for Nested Structural Epidemiological Data

Wen Wen，Wen Xiaoyan，Hu Shan，et al.（Department of Health Statistics，School of Public Health and Management，Chongqing Medical University（400016），Chongqing）

ObjectiveTo develop two hierarchical Bayesian models for the epidemiological data with focusing on its nested structure；as well as to explore the pros and cons of them.MethodsRelationships among nested layers of nested structural data are taken into account when developing the two hierarchical Bayesian models.The first model focuses on the stratification effect of each nested layer for differentiation between the layers.The second model focuses on the transmission effect between the father layer and its son layers.Open BUGS software and a birth defects survey data were used to fit and evaluate the two hierarchical Bayesian models；and the deviance information criterion（DIC）was used for measuring the goodness-of-fit of them.A sensitivity analysis was conducted with different sets of prior information on hyper parameter of the population rateμ.ResultsThe DIC of the two models are 101.8 and 101.6，respectively，which shows almost the same goodness-of-fit of them.The sensitivity analysis shows that the standard deviation of the two models for the posterior mean of estimated population rate are（10－4）1.191and 27.546，respectively，for the posterior median of them are（10－4）1.038 and 7.617，respectively.Both results of posterior mean and posterior median say that the first model has smaller standard deviation under different prior information scenario.ConclusionBoth models can be used to model nested structural epidemiological data.However，the first model is affected by prior information much less than the second model does.Thus，the first model is more stable and is better to model nested structural survey data when little prior information is available.

Nested structural data；Stratified sampling；Hierarchical Bayesian model；OpenBUGS

國家自然科學(xué)基金（81373103）、重慶市科委基礎(chǔ)與前沿研究計劃項目（cstc2013jcyjA10009）

△通信作者：彭斌，Email：pengbin＠cqmu.edu.cn