蘭州大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)研究所(730000)

目前國內(nèi)研究較多的是沈可[1]等人實現(xiàn)OpenBUGS軟件在網(wǎng)狀Meta分析中的應(yīng)用，鄭曉鴛[7]實現(xiàn)OpenBUGS軟件在跳擴散Shibor模型中的應(yīng)用。本文利用非條件Logisitic回歸模型，研究少年兒童肥胖與膽固醇、甘油三酯、年齡和性別之間的關(guān)系，來實現(xiàn)OpenBUGS在非條件Logisitic回歸中的應(yīng)用。數(shù)據(jù)[8]見表2所示。

·計算機應(yīng)用·

OpenBUGS軟件介紹及應(yīng)用*

蘭州大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)研究所(730000)

張繼巍 高文龍 秦天燕 劉建正 李學(xué)朝 拉扎提木拉提 李娟生△

OpenBUGS軟件是在WinBUGS軟件基礎(chǔ)上研制的一款實現(xiàn)貝葉斯統(tǒng)計推斷的工具軟件，它是以MCMC(Markov Chain Monte Carlo)方法為基礎(chǔ)，將所有未知或不確定的參數(shù)都視為隨機變量，并對此種類型的概率模型進行求解[1]。它廣泛地應(yīng)用于醫(yī)學(xué)、經(jīng)濟學(xué)、生命科學(xué)、心理學(xué)、社會科學(xué)等多個領(lǐng)域[1-4]。目前已經(jīng)更新到OpenBUGS3.2.3版本，可以從其官網(wǎng)(http://www.openbugs.net/w/Downloads)上免費下載使用[5]。但是目前國內(nèi)關(guān)于這款軟件的應(yīng)用報道比較少，對軟件的基本情況也沒有一個詳細的介紹。因此，本文主要對該軟件的功能和具體操作做一個簡單的匯總和介紹，希望有興趣的研究者能較快掌握這個軟件，并起到拋磚引玉的作用。

OpenBUGS軟件及特點

1.基本原理： 通過Gibbs 抽樣和Metropolis算法，從完全條件概率分布中抽樣，從而生成馬爾科夫鏈，通過迭代，最終估計出模型的參數(shù)。

2.作用及過程：OpenBUGS 軟件是通過Gibbs抽樣和MCMC(Markov Chain Monte Carlo)方法從任意復(fù)雜模型的后驗分布中產(chǎn)生樣本，提供了一個有效的方法估計貝葉斯統(tǒng)計模型，極大地推進了貝葉斯的應(yīng)用。因為它將所有的未知參數(shù)都視為隨機變量，故可以方便的對許多缺失數(shù)據(jù)、異常數(shù)值、小樣本資料和非參數(shù)模型及分布進行Gibbs抽樣，實現(xiàn)Gibbs抽樣的建模，并對其進行擬合和分析，最終得到模型參數(shù)的估計，并給出參數(shù)的Gibbs抽樣動態(tài)圖、核密度圖、均數(shù)及其置信區(qū)間的變化圖等，使抽樣結(jié)果更直觀、可靠。Gibbs抽樣收斂后，還可以方便地得到參數(shù)后驗分布的均數(shù)、標準差、95%置信區(qū)間和中位數(shù)等信息[6]。

3.分布類型：OpenBUGS中提供了29種常用的分布(用于指定似然函數(shù)和先驗分布)，主要有離散型單變量分布(discrete univariate distribution)、連續(xù)型單變量分布(continuous univariate distribution)、離散型多變量分布(discrete multivariate distribution)、連續(xù)型多變量分布(continuous multivariate distribution)等四種。其中最常用的有17種，如表1所示。

表1 OpenBUGS中的常用分布

4.基本操作

OpenBUGS的基本操作主要包括Doodle、Model、Attributes、Update、Inference、Info、Stats、DIC和Log等命令。

(1)Doodle：用于構(gòu)建Doodle模型圖，以節(jié)點、箭頭和平板的方式出現(xiàn)。

(2)Model：通過Specificaion用于對模型進行校驗，包括模型的檢查、數(shù)據(jù)的加載、編譯和初始值的生成。

(3)Attributes：用于改變所選文本的字體、格式、字體的大小和顏色等。

(4)Update：用于對Doodle模型進行迭代，可以設(shè)置迭代的次數(shù)和步長。

(5)Inference：用于指定想要考察的變量和迭代計算結(jié)果的輸出。

(6)Info：用于獲得節(jié)點和數(shù)據(jù)的信息。

(7)Stats：用于輸出參數(shù)后驗分布的統(tǒng)計量，比如均數(shù)、標準差、95%置信區(qū)間和中位數(shù)等。

(8)DIC：用于對多個待定模型的擬合優(yōu)度進行比較選擇，DIC值越小，模型的擬合度越好。

(9)Log：用于查看以上所有操作的具體信息。

應(yīng)用實例

目前國內(nèi)研究較多的是沈可[1]等人實現(xiàn)OpenBUGS軟件在網(wǎng)狀Meta分析中的應(yīng)用，鄭曉鴛[7]實現(xiàn)OpenBUGS軟件在跳擴散Shibor模型中的應(yīng)用。本文利用非條件Logisitic回歸模型，研究少年兒童肥胖與膽固醇、甘油三酯、年齡和性別之間的關(guān)系，來實現(xiàn)OpenBUGS在非條件Logisitic回歸中的應(yīng)用。數(shù)據(jù)[8]見表2所示。

表2 少年兒童肥胖危險因素調(diào)查資料表

引入自變量X1[i]，X2[i]，X3[i]，X4[i]，其中

用ri表示第i組少年兒童肥胖的陽性數(shù)，ni表示第i組觀察總例數(shù)，pi表示第i組兒童肥胖的陽性率，則有ri～Binomial(pi,ni)。

建立模型 logit(pi)=β0+β1X1[i]+β2X2[i]+β3X3[i]+β4X4[i]

其中beta用β表示；β0為常數(shù)項，又稱截距；βj(j=1,2,3,4)為自變量Xj的偏回歸系數(shù)；模型的參數(shù)β0，β1，β2，β3和β4已經(jīng)給定了獨立的“non-informative”先驗分布。利用OpenBUGS軟件建模和抽樣，步驟如圖1所示：

圖1 OpenBUGS操作流程圖

1.模型的建立(new model)：通過 OpenBUGS軟件進行Doodle建模，用來指定各參數(shù)的分布類型和邏輯關(guān)系，如圖2所示：

圖2 Doodle模型結(jié)構(gòu)示意圖

圖2中每一個橢圓表示一個節(jié)點，矩形表示常數(shù)節(jié)點，單線箭頭表示從父節(jié)點到隨機型子節(jié)點，雙線空心箭頭表示從父節(jié)點到邏輯型子節(jié)點，邊緣較粗的矩形為平板，其左下角的 “for(i IN 1:N)”表示for循環(huán)，實質(zhì)是計算N個組的似然函數(shù)，從而獲得整個數(shù)據(jù)的似然函數(shù)。

2.模型的檢查(check model)：模型設(shè)定完成后，點擊Model下拉菜單Specification中的check model 按鈕對模塊進行檢查，若檢查無誤，則軟件窗口左下方會提示“Model is syntactically correct”，同時按鈕load data和compile被激活，提示可以加載數(shù)據(jù)并進行編譯。

3.數(shù)據(jù)的加載(load data)和編譯(compile model)：本例數(shù)據(jù)為

list(x1=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2),

x2=c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4),

x3=c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2),

x4=c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2),

r=c(8,1,0,0,12,8,1,12,5,9,2,4,22,4,0,1,7,0,2,0,12,7,2,5,2,7,2,1,4,3,0,0),

n=c(98,16,2,2,83,18,22,39,75,23,13,19,232,55,5,2,93,19,4,2,66,22,18,32,75,28,5,10,185,82,3,4),N=32)

將鼠標移到數(shù)據(jù)指示語“l(fā)ist”并選中，點擊load data(加載數(shù)據(jù))按鈕，若左下方提示“data loaded”，則點擊compile按鈕進行編譯，編譯成功后，窗口底部左下方會提示“model compiled”，隨即load-inits和gen inits按鈕被激活 。

4.初始值的設(shè)定(initial values)：本例初始值為

list(beta0=0,beta1=0,beta2=0,beta3=0,beta4=0)

與加載數(shù)據(jù)類似，選中指示語 “l(fā)ist”，點擊load inits按鈕加載(如果只產(chǎn)生一條鏈，此時窗口左下方出現(xiàn)“model is initialized”，gen inits按鈕變灰，提示初始值設(shè)置完成；如果產(chǎn)生2條或2條以上的鏈，此時窗口左下方提示“initial values loaded but chain contains uninitialized variables”，繼續(xù)點擊gen load按鈕，提示“initial values generated,model is initialized ”)，初始值設(shè)定完成。

5.模型的退火(burn in)：點擊Inference下拉菜單中的Samples按鈕，在彈出的 “Sample Monitor Tool” 對話框中的beg處輸入相應(yīng)的數(shù)字，如輸入1000則表示拋去前1000次抽樣以消除初始值對抽樣的影響，實現(xiàn)對模型的退火，從1001次開始抽樣。

6.變量的監(jiān)控(monitor nodes)：通過點擊Inference下拉菜單中的samples，在彈出的對話框node處，添加自己想要觀察的變量名(本例主要觀察的變量為β0，β1，β2，β3，β4)，逐一按set鍵確定，同時可以在percentiles下方的列表中選擇給出的置信區(qū)間，默認為95%(val2.5pc，val97.5pc)置信區(qū)間。

7.模型的迭代(updates )：通過點擊Model下拉菜單中的updates按鈕對模型進行迭代運算，迭代的次數(shù)和步長可以自行設(shè)置。

8.結(jié)果的輸出(output)：通過Inference下拉菜單中Samples選項，在打開的“Sample Monitor Tool”對話框中的node處輸入前面指定的觀察參數(shù)，也可以直接輸入“*”(“*”代表所有指定的參數(shù))，就可以獲得各參數(shù)后驗分布統(tǒng)計量，如圖3～圖8所示。

圖3 各參數(shù)后驗分布統(tǒng)計量

圖4 核密度圖

圖5 迭代軌跡圖

圖6 分位數(shù)圖

圖7 自相關(guān)函數(shù)圖

其中：MC＿error表示蒙特卡羅模擬的誤差，用來度量由模擬引起的參數(shù)均值的方差；val2.5pc和 val97.5pc分別表示中位數(shù)的95%置信區(qū)間的下限和上限；median表示中位數(shù)，通常比mean更穩(wěn)定；start表示Gibbs抽樣得起始點，本例拋去前1000次用于消除初始值對抽樣的影響，從1001次后開始抽樣；sample表示總共抽取的樣本數(shù)。

圖8 迭代歷史圖

9.收斂性的判別(convergence)：模型的收斂性可以通過觀察遍歷均值進行判斷，在得到的鏈中每隔一段距離計算所觀察參數(shù)的遍歷均值，當這樣算的均值穩(wěn)定后，可認為模型收斂；也可以通過迭代軌跡圖和迭代歷史圖來進行判斷，當?shù)壽E和歷史迭代趨于穩(wěn)定，可以認為模型收斂[6]。 由圖5，圖8可知模型基本收斂，得到的回歸方程為：Logit(pi)=-2.14-0.52X1-0.23X2+0.26X3+0.79X4。

討 論

本文僅僅對OpenBUGS軟件做了一個簡單的介紹，通過實現(xiàn)OpenBUGS軟件在非條件Logistic回歸分析中的應(yīng)用，使讀者對該軟件有一個初步的了解，并能掌握OpenBUGS軟件的基本操作，變量的分布類型和作用以及熟悉OpenBUGS軟件的操作流程。在運用該軟件進行貝葉斯統(tǒng)計分析時，最重要的是對模型似然函數(shù)和參數(shù)先驗分布的指定[9]，這就要求使用者熟練掌握貝葉斯統(tǒng)計推斷的理論知識，結(jié)合研究的具體內(nèi)容，反復(fù)摸索。

OpenBUGS軟件在克服WinBUGS缺點的基礎(chǔ)上，對模型的代碼和操作菜單做了相關(guān)優(yōu)化和調(diào)整，使運行更加平穩(wěn)，操作更為簡單[1]。與其他貝葉斯分析軟件(如BACC/BMA)相比，OpenBUGS軟件的亮點就在于其具有很強的靈活性，能夠使貝葉斯分析中復(fù)雜的數(shù)值計算簡單化。但其在網(wǎng)狀關(guān)系圖、森林圖等圖形的繪制方面仍存在著不足[1]，并且在代碼建模過程中容易出現(xiàn)錯誤而不易被發(fā)現(xiàn)。目前，OpenBUGS還在更新和維護當中，其強大的數(shù)據(jù)分析能力，在不久的將來，會成為貝葉斯統(tǒng)計分析的主流軟件。

[1]沈可,王芬,張超,等.應(yīng)用OpenBUGS軟件實現(xiàn)網(wǎng)狀Meta分析.湖北醫(yī)藥學(xué)報,2013,32(6):476-479.

[2]Carroll R,Lawson AB,Faes C,et al.Comparing INLA and OpenBUGS for hierarchical Poisson modeling in disease mapping.Spatial and Spatio-temporal Epidemiology,2015,14(15):5-54.

[3]Lyle W,Konigsberg,Frankenberg.Bayes in Biological Anthropology.American Journal of Physical Anthropology,2013,152(57):153-184.

[4]Eitzel M,Battles J,York R,et al.Estimating tree growth from complex forest monitoring data.Ecological Applications,2013,23(6):1288-1296.

[5]OpenBUGS 3.2.3 user manual.

[6]孟海英,劉桂芬,羅天娥.Winbugs軟件應(yīng)用.中國衛(wèi)生統(tǒng)計,2006,23(4):375-377.

[7]鄭曉鴛.MCMC方法在跳擴散Shibor模型參數(shù)模擬中的應(yīng)用.時代金融,2014,562(8):31-33.

[8]孫振球,徐勇勇,主編.醫(yī)學(xué)統(tǒng)計學(xué).第4版.北京:人民衛(wèi)生出版社,2014,244-261.

[9]曾平,王婷,何鵬.非標準分布貝葉斯分析的WinBUGS軟件實現(xiàn).中國衛(wèi)生統(tǒng)計,2012,29(4):614-617.

(責任編輯：郭海強)

教育部人文社會科學(xué)研究西部和邊疆地區(qū)項目(項目批準號：15XJC910001)

△通信作者:李娟生，E-mail:ljsh@lzu.edu.cn