貝葉斯統(tǒng)計分析的有力工具
——OpenBUGS軟件*

東南大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)系（210009） 楊維維 劉 沛 巢健茜 楊 靚

貝葉斯統(tǒng)計分析的有力工具
——OpenBUGS軟件*

東南大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)系（210009） 楊維維 劉 沛△巢健茜 楊 靚

貝葉斯統(tǒng)計學(xué)與頻率統(tǒng)計學(xué)是當(dāng)今世界主要的兩大統(tǒng)計學(xué)派，二者在統(tǒng)計推斷的理論和方法上存在較大差異。隨著現(xiàn)代計算機(jī)的高速發(fā)展，貝葉斯統(tǒng)計的研究不再只停滯于理論階段，馬爾可夫鏈蒙特卡羅（Markov Chain Monte Carlo，MCMC）方法的應(yīng)用，解決了后驗(yàn)分布復(fù)雜高維計算的瓶頸問題，使得貝葉斯統(tǒng)計在理論和方法上均取得了快速發(fā)展［1］。特別是基于貝葉斯統(tǒng)計理論的BUGS軟件、W inBUGS軟件和OpenBUGS軟件的成功開發(fā)，為貝葉斯統(tǒng)計分析提供了強(qiáng)有力工具。

BUGS是Bayesian inference using gibbs sampling的縮寫，最初由英國醫(yī)學(xué)研究協(xié)會（medical research council，MRC）于1989年研制，適用于 DOS操作系統(tǒng)［2－3］。W inBUGS軟件是 BUGS軟件在 W indows操作系統(tǒng)下的版本，最早出現(xiàn)于1989年，最新版本為1.4.3，于 2007年停止更新［3－4］。2004年，Andrew Thomas在赫爾辛基大學(xué)開始研制OpenBUGS軟件，經(jīng)過十余年的開發(fā)，目前這一軟件已越來越顯示出其獨(dú)特的優(yōu)勢。雖然該軟件基于的理論方法與W in-BUGS軟件基本相同［3］，但用戶能夠根據(jù)自己的需求進(jìn)行二次開發(fā)，也無需考慮軟件操作系統(tǒng)問題。本文將對目前國際應(yīng)用較廣而國內(nèi)介紹較少的OpenBUGS軟件從基本情況、操作步驟和實(shí)例應(yīng)用方面進(jìn)行介紹，并將其與目前國內(nèi)應(yīng)用較多的W inBUGS軟件進(jìn)行比較。

OpenBUGS軟件概述

OpenBUGS軟件是采用MCMC方法對復(fù)雜統(tǒng)計模型進(jìn)行貝葉斯推斷的專業(yè)工具，由Pascal語言編寫源代碼并且開放源代碼。OpenBUGS軟件可以在Windows、Unix、Linux（包括通過 wine模擬器軟件）操作系統(tǒng)下使用，也可以通過R軟件的程序包（例如R2OpenBUGS）調(diào)用 OpenBUGS軟件實(shí)現(xiàn)統(tǒng)計分析［5］。OpenBUGS軟 件 可 從 http://www.openbugs.net/w/Downloads免費(fèi)下載，最新版本為2014年3月15日發(fā)布的3.2.3，不需要注冊即可使用軟件的全部功能。

OpenBUGS軟件基本操作步驟

在采用OpenBUGS軟件構(gòu)建貝葉斯模型時，其操作流程與 W inBUGS軟件基本一致［6－8］，主要分為 6個步驟:

步驟1:模型的構(gòu)建與數(shù)據(jù)的輸入

在文件窗口中編寫模型程序，完成數(shù)據(jù)的輸入并賦予初始值，F(xiàn)ile菜單中選擇“Save AS”將編寫好的程序保存到文檔中（后綴名為.odc）。

步驟2:模型的定義

選擇菜單Model＼Specification，光標(biāo)移動到模型框架內(nèi)或者選中“model”，單擊對話框中的“check model”按鈕，若無語法錯誤，窗口底部將顯示“model is syntactically correct”，然后依次進(jìn)行加載數(shù)據(jù)（loda data）、給定模擬鏈數(shù)、編譯程序（compile）、加載初始值（lodainits）或者由系統(tǒng)自動產(chǎn)生初始值（gen inits），窗口底部將顯示“model is initialized”。

步驟3:考察參數(shù)的選定

選擇菜單 Inference/Samples，在 Sample Monitor Tool對話框中的node處輸入需要考察的參數(shù)，每輸入一個參數(shù)名均單擊set，在node中輸入“*”即指定所有需考察的未知參數(shù)，trace、history等按鈕點(diǎn)擊后不能打開相應(yīng)的窗口，需在迭代更新后打開。

步驟4:迭代運(yùn)算

選擇菜單 Model/Update，在 updates中輸入 MCMC預(yù)迭代次數(shù)，單擊update按鈕開始模擬運(yùn)算，若要中途停止更新，再次單擊update按鈕即可。

步驟5:收斂性診斷

點(diǎn)擊Sample Monitor Tool對話框中history按鈕觀察迭代歷史圖，trace按鈕給出Gibbs動態(tài)抽樣蹤跡圖，如果迭代歷史圖和蹤跡圖趨于穩(wěn)定，說明收斂性較好，不收斂則重復(fù)步驟4，增加迭代次數(shù)，若迭代很多次后仍不收斂，則需考慮對模型進(jìn)行相應(yīng)的修改。

步驟6:后驗(yàn)分析

在beg中輸入丟棄初始迭代結(jié)果的次數(shù)，減少初始值的影響，單擊Sample Monitor Tool對話框中stats按鈕輸出后驗(yàn)參數(shù)的描述性統(tǒng)計量，包括均數(shù)、中位數(shù)、標(biāo)準(zhǔn)差、MC誤差等。coda按鈕可將模擬結(jié)果保存到外部文件，供R等軟件進(jìn)一步分析和做圖使用。

OpenBUGS軟件與W inBUGS軟件區(qū)別

OpenBUGS軟件的開發(fā)主要以操作系統(tǒng)運(yùn)行環(huán)境的擴(kuò)充、其他軟件的調(diào)用和開放源代碼為主要目的，同時，OpenBUGS軟件的開發(fā)應(yīng)用并不會對已安裝W in-BUGS軟件的用戶產(chǎn)生任何危害［9］。W inBUGS1.4.3軟件可以免費(fèi)獲得，但用戶必須安裝W inBUGS軟件無約束適用的“key”，否則只能運(yùn)行100個節(jié)點(diǎn)以下的模型，而OpenBUGS軟件則不需要，用戶可以免費(fèi)獲得其全部的功能。除了在授權(quán)和操作系統(tǒng)范圍方面不一樣外，兩者之間其他方面的區(qū)別如下［10－12］:

1.二者間最主要的不同在于模型更新算法方面，OpenBUGS軟件的 Updater Options包含“updaters”和“parameters”兩個標(biāo)簽（圖1），updaters標(biāo)簽下的列表框給出大量可用的更新算法系統(tǒng)名稱，可以改變當(dāng)前node的更新方法（在模型編譯后），并且可通過菜單Info/Updaters（by name）來查看各 node算法更新的相關(guān)信息；parameters標(biāo)簽下對話框給出為解決當(dāng)前問題所采用的更新方法的系統(tǒng)名稱，可改變算法的默認(rèn)參數(shù)值（如iterations等），而WinBUGS軟件的更新選項(xiàng)界面如圖2所示，其只有“parameter”對話框的作用，可見OpenBUGS軟件擁有更好的靈活性和可擴(kuò)展性［5，10－11］。

圖1 OpenBUGS軟件中Updater options界面

圖2 WinBUGS軟件中Updater options界面

2.OpenBUGS軟件的Model菜單除了增加“Updater Options”外，還增加了 Externalize和 Internalize兩個選項(xiàng)，Externalize可以將編譯好模型的所有內(nèi)部數(shù)據(jù)結(jié)構(gòu)信息以bug文件形式保存在OpenBUGS軟件安裝目錄下的“Restar”文件夾中，Internalize選項(xiàng)可以從“Restar”文件夾讀入編譯好的模型數(shù)據(jù)結(jié)構(gòu)信息，這兩個選項(xiàng)為保存和調(diào)用模型提供了便捷，不需要再對模型從頭進(jìn)行模型驗(yàn)證、數(shù)據(jù)加載、編譯等步驟，直接可以進(jìn)行迭代和觀測后驗(yàn)參數(shù)的描述性統(tǒng)計量，當(dāng)然，Externalize和Internalize選項(xiàng)也會產(chǎn)生很大的文件。Model菜單中的Latex選項(xiàng)將給出Latex代碼來描述當(dāng)前的模型。

3.OpenBUGS軟件中Script命令的語法與W in-BUGS軟件的 Script有很大差別，與其在R中調(diào)用OpenBUGS所產(chǎn)生的腳本命令相似，通過Script命令可以自動進(jìn)行一些常規(guī)分析和建立可重復(fù)的分析，通過R軟件可以將W inBUGS軟件中的腳本轉(zhuǎn)化為OpenBUGS軟件的腳本。為了更好地展現(xiàn)不同nodes迭代的信息，OpenBUGS軟件的Info菜單中增加了很多選項(xiàng)，例如:Updaters by name和Updaters by depth，將nodes分別按名稱和更新深度列出。選擇Info菜單中的Show data選項(xiàng)，將打開一個窗口顯示模型中的數(shù)據(jù)部分。

4.OpenBUGS軟件新增了一些函數(shù)功能和分布，例 如 eigen.vals， integral， gammap， ode， prod，p.valueM，solution等等，同時也對一些分布進(jìn)行重命名使其符合R中的相應(yīng)分布規(guī)范命名，例如把狄里雷特分布從 ddirch改為 ddirich等。除此之外，Open-BUGS軟件在處理刪失C（lower，upper）和截尾分布T（lower，upper）采用的函數(shù)不一樣，而 W inBUGS軟件則均通過I（lower，upper）來處理刪失和截尾分布，通常不贊成這種做法。

5.OpenBUGS軟件Help菜單中增加了索引選項(xiàng)，可以輸入關(guān)鍵詞在指南和案例中快速檢索，將函數(shù)功能和分布單獨(dú)作為選項(xiàng)，并且增加了Examples菜單，提供更多的使用例子，這些改變都更方便用戶的使用和學(xué)習(xí)。

應(yīng)用舉例

本研究以貝葉斯混合線性模型為例介紹Open-BUGS軟件的應(yīng)用及實(shí)現(xiàn)，例子來源于DenHond E，Lesaffre E等《The Inter-regional Belgian Bank Employee Nutrition Study（IBBENS）》［13－14］一文，目的是比較年齡與性別對8家銀行子公司成員膽固醇攝入量的影響，其中7家位于 Flanders，1家位于 Wallonie。將年齡與性別作為協(xié)變量，構(gòu)建貝葉斯混合線性模型如下:

式中yij為反應(yīng)變量，β0為基線時研究對象yij的均數(shù)，β1為協(xié)變量年齡參數(shù)，β2為協(xié)變量性別參數(shù)，b0i為隨機(jī)效應(yīng)參數(shù)，εij為隨機(jī)誤差。

將上述貝葉斯混合模型寫成如下的OpenBUGS程序代碼:

model＃OpenBUGS代碼以 model開始，｛…｝為模型程序部分

｛

for（ic in 1:C）｛

for（j in 1:nc［ic］）｛

chol［cumnc［ic］＋j］～dnorm（mu［cumnc［ic］＋j］，tau.i）＃膽固醇攝入量服從正態(tài)分布

mu［cumnc［ic］＋j］＜－beta［1］＋beta［2］*（age［cumnc［ic］＋j］-mean（age［］））＋beta［3］*gender［cumnc［ic］＋j］＋b0［ic］｝＃混合線性模型

b0［ic］～dnorm（0.0，tau.b02）

｝

tau.i～dgamma（0.001，0.001）＃設(shè)方差倒數(shù)為無信息先驗(yàn) Ga（0.001，0.001）

sigma.i＜－1/sqrt（tau.i）＃通過方差的倒數(shù)求解方差和標(biāo)準(zhǔn)差

sigma.i2＜－1/tau.i

for（i in 1:3）｛beta［i］～dnorm（0.0，1.0E－6）｝＃設(shè)定3個回歸系數(shù)為無信息先驗(yàn)N（0，106）

sigma.b0～dunif（0，100）＃設(shè)隨機(jī)效應(yīng)參數(shù)為無信息先驗(yàn) U（0，100）

sigma.b02＜-pow（sigma.b0，2）

tau.b02＜-pow（sigma.b02，-1）

｝

注:C為子公司數(shù)，Nc為子公司的人數(shù)，Cumnc為子公司累計人數(shù)，Chol為每個人的膽固醇攝入量。

1.收斂性的判別

所有基于MCMC結(jié)果的貝葉斯推斷是在假設(shè)鏈已經(jīng)收斂的前提下進(jìn)行的，可見馬爾可夫鏈的收斂性判別在貝葉斯統(tǒng)計推斷中非常關(guān)鍵［8］。本文采用的例子中對各參數(shù)均取以無信息先驗(yàn)，同時模擬3條鏈并給出3組初始值，設(shè)定預(yù)迭代模擬5000次，迭代完成后得到一系列診斷收斂性的圖，包括蹤跡圖、Gelman-Rubin診斷圖和自相關(guān)圖。各參數(shù)的蹤跡圖顯示，三條鏈很好地混合在一起，且圖形呈水平趨勢，波動幅度大致一樣，趨于穩(wěn)定，各參數(shù)的Gelman-Rubin診斷圖顯示在迭代5000次后各參數(shù)的三條鏈均處于平穩(wěn)狀態(tài)，統(tǒng)計量R所對應(yīng)的線基本保持在1處，并且呈水平線。同時，各參數(shù)的自相關(guān)函數(shù)圖很快接近于0，這些特征均說明鏈已經(jīng)趨于收斂。

圖3 對各參數(shù)后驗(yàn)分布的估計

2.貝葉斯統(tǒng)計推斷

模型擬合穩(wěn)定后正式迭代10000次用于各參數(shù)的估計，5000次預(yù)迭代將用以退火，得到估計結(jié)果如圖3和表1，beta［］和sigma.i的后驗(yàn)核密度函數(shù)圖形呈正態(tài)分布，后驗(yàn)中位數(shù)估計與后驗(yàn)均值估計重合，sigma.b0呈偏度分布，取中位數(shù)作為貝葉斯估計效果更好。在舍去預(yù)迭代5000次后，beta［2］和 beta［3］的后驗(yàn)均值（標(biāo)準(zhǔn)差）分別為－0.6919（0.5759）和 －62.64（10.71），95%可信區(qū)間分別為（－1.82，0.4315）和（－83.67，－41.47），說明在考慮到公司效應(yīng)的情況下，協(xié)變量年齡對膽固醇攝入量沒有影響，而不同性別的膽固醇攝入量之間存在統(tǒng)計學(xué)差異，女性膽固醇攝入量低于男性。

表1 對各參數(shù)的后驗(yàn)估計值

討 論

雖然，OpenBUGS軟件與W inBUGS軟件擁有相似的界面和操作過程，但是作為BUGS軟件工程現(xiàn)在和未來開發(fā)的重點(diǎn)，OpenBUGS軟件已經(jīng)顯示出其優(yōu)越性、成熟性和穩(wěn)定性。從開放性講，OpenBUGS是一款開放源代碼的軟件，用戶可根據(jù)自己的需求進(jìn)行二次開發(fā)，而W inBUGS軟件只是免費(fèi)卻不開放其源代碼，并且已經(jīng)停止更新，其可擴(kuò)展性遠(yuǎn)不如OpenBUGS軟件。從靈活性說，相比W inBUGS軟件，OpenBUGS軟件提供大量的解決當(dāng)前問題所需的更新方法，同時可以將已經(jīng)編譯好的模型信息保存，在下次使用時避免操作的重復(fù)性。除此之外，OpenBUGS在操作系統(tǒng)應(yīng)用范圍和運(yùn)行時間效率方面都優(yōu)于W inBUGS軟件，在 OpenBUGS軟件網(wǎng)站上已給出 OpenBUGS 3.2.3與W inBUGS 1.4.3關(guān)于59個軟件案例運(yùn)行時間的比較結(jié)果，其中50個案例在OpenBUGS軟件中運(yùn)行時間短于W inBUGS軟件，有的時間甚至達(dá)到6倍左右。

同時，OpenBUGS軟件也存在一定的不足。首先，雖然其已經(jīng)做到允許數(shù)據(jù)部分不出現(xiàn)在模型中而保存到數(shù)據(jù)文件里，但是其不能直接調(diào)用外部數(shù)據(jù)文件，在數(shù)據(jù)格式方面有嚴(yán)格的要求，這使得在處理大樣本數(shù)據(jù)時非常不便捷并且易出現(xiàn)錯誤，目前解決這一問題的主要方法是在其他軟件中調(diào)用OpenBUGS軟件。其次，OpenBUGS軟件與W inBUGS軟件在計算偏差信息準(zhǔn)則（DIC）方面能力均較弱，例如在建立混合模型情況下，DIC值是無法計算的，其菜單選項(xiàng)呈現(xiàn)灰色狀態(tài)。這些不足的地方也正是OpenBUGS軟件未來發(fā)展需要解決的?？偟膩碚f，OpenBUGS軟件有著強(qiáng)大的開放性、可擴(kuò)展性、靈活性以及操作系統(tǒng)應(yīng)用范圍廣等優(yōu)勢。目前，其在我國的應(yīng)用還處于起步階段，本文只是對OpenBUGS軟件進(jìn)行初步介紹，相信只有在進(jìn)一步實(shí)際應(yīng)用中，才能不斷加深對其的理解，從而發(fā)揮其在貝葉斯統(tǒng)計分析中的強(qiáng)大功能。

［1］劉沛，陳啟光.貝葉斯統(tǒng)計及其在診斷和篩檢試驗(yàn)評價中的應(yīng)用.中國衛(wèi)生統(tǒng)計，2006，23（4）:361-363.

［2］GilksWR，Thomas A，Spiegel halter DJ，et al.A Language and Program for Complex Bayesian Modelling.Journal of the Royal Statistical Society.Series D（The Statistician），1994，43（1）:169-177.

［3］Serang S，Zhang Z，Helm J，et al.Evaluation of a Bayesian Approach to Estimating Nonlinear Mixed-Effects Mixture Models.Structural E-quation Modeling:A Multidisciplinary Journal，2014，22（2）:202-215.

［4］茆詩松，湯銀才主編.貝葉斯統(tǒng)計.北京:中國統(tǒng)計出版社，2012，275-280.

［5］ OpenBUGS website.Overview.2010.http://www.openbugs.net/w/Overview.

［6］董圣杰，冷衛(wèi)東，田家祥，等.Meta分析系列之五:貝葉斯Meta分析與W inBUGS軟件.中國循證心血管醫(yī)學(xué)雜志，2012，（5）:395-398.

［7］孟海英，劉桂芬，羅天娥，等.W inBUGS軟件應(yīng)用.中國衛(wèi)生統(tǒng)計，2006，（4）:375-377.

［8］饒克勤主編.衛(wèi)生統(tǒng)計方法與應(yīng)用進(jìn)展（第2卷）.北京:人民衛(wèi)生出版社，2008:292-308.

［9］Lunn D，Spiegel halter D，Thomas A，etal.The BUGS project:Evolution，critique and future directions.Statistics in Medicine，2009，28（25）:3049-3067.

［10］Spiegel halter D，Thomas A，Best N，et al.WinBUGS User Manual，Version 1.4，2003.

［11］Spiegel halter D，Thomas A，Best N，et al.OpenBUGS User Manual，Version 3.2.3，2014.

［12］Changes Between WinBUGS and Open BUGS.http://www.openbugs.net/w/OpenVsWin.

［13］Hond ED，Muylaert A，Lesaffre E，et al.The Inter-regional Belgian Bank Employee Nutrition Study（IBBENS）.Eur JClin Nutr，1994，48（2）:106-117.

［14］Lesaffre E.Bayesian Biostatistics.Wiley.com，2012.

國家自然科學(xué)基金（81273189）；“中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金”和“江蘇省普通高校研究生科研創(chuàng)新計劃資助”項(xiàng)目（SJLX＿0108）

△通信作者:劉沛，E-mail:liupeiseu@126.com

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