貝葉斯統計分析的有力工具
——OpenBUGS軟件*

東南大學公共衛生學院流行病與衛生統計學系（210009） 楊維維 劉 沛 巢健茜 楊 靚

貝葉斯統計分析的有力工具
——OpenBUGS軟件*

東南大學公共衛生學院流行病與衛生統計學系（210009） 楊維維 劉 沛△巢健茜 楊 靚

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

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

OpenBUGS軟件概述

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

OpenBUGS軟件基本操作步驟

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

步驟1:模型的構建與數據的輸入

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

步驟2:模型的定義

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

步驟3:考察參數的選定

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

步驟4:迭代運算

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

步驟5:收斂性診斷

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

步驟6:后驗分析

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

OpenBUGS軟件與W inBUGS軟件區別

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

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

圖1 OpenBUGS軟件中Updater options界面

圖2 WinBUGS軟件中Updater options界面

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

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

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

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

應用舉例

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

式中yij為反應變量，β0為基線時研究對象yij的均數，β1為協變量年齡參數，β2為協變量性別參數，b0i為隨機效應參數，εij為隨機誤差。

將上述貝葉斯混合模型寫成如下的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）＃膽固醇攝入量服從正態分布

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）＃設方差倒數為無信息先驗 Ga（0.001，0.001）

sigma.i＜－1/sqrt（tau.i）＃通過方差的倒數求解方差和標準差

sigma.i2＜－1/tau.i

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

sigma.b0～dunif（0，100）＃設隨機效應參數為無信息先驗 U（0，100）

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

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

｝

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

1.收斂性的判別

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

圖3 對各參數后驗分布的估計

2.貝葉斯統計推斷

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

表1 對各參數的后驗估計值

討 論

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

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

［1］劉沛，陳啟光.貝葉斯統計及其在診斷和篩檢試驗評價中的應用.中國衛生統計，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］茆詩松，湯銀才主編.貝葉斯統計.北京:中國統計出版社，2012，275-280.

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

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

［7］孟海英，劉桂芬，羅天娥，等.W inBUGS軟件應用.中國衛生統計，2006，（4）:375-377.

［8］饒克勤主編.衛生統計方法與應用進展（第2卷）.北京:人民衛生出版社，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.

國家自然科學基金（81273189）；“中央高校基本科研業務費專項資金”和“江蘇省普通高校研究生科研創新計劃資助”項目（SJLX＿0108）

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

（責任編輯:郭海強）