董圣杰，冷衛(wèi)東，田家祥，曾憲濤

貝葉斯Meta分析（Bayesian Meta-Analysis）是近年來基于貝葉斯統(tǒng)計(jì)發(fā)展起來的一種新型的Meta分析方法，主要采用“馬爾科夫鏈—蒙特卡羅”（Markov chain Monte Carlo，MCMC）方法、使用WinBUGS軟件[1]進(jìn)行。經(jīng)典統(tǒng)計(jì)學(xué)派的統(tǒng)計(jì)量，往往不易找到其精確的有限樣本分布，因此多數(shù)情況下是基于大樣本漸近分布做出統(tǒng)計(jì)推斷，而貝葉斯學(xué)派則可直接計(jì)算精確的有限樣本分布，并不依賴于漸近理論，且充分考慮了模型的不確定性，故認(rèn)為Meta分析貝葉斯估計(jì)更可靠、更合理，特別在有序數(shù)據(jù)及網(wǎng)狀Meta分析[2]中有傳統(tǒng)Meta分析無法企及的優(yōu)點(diǎn)。當(dāng)前，貝葉斯Meta分析已得到愈發(fā)廣泛的應(yīng)用，本文將簡要介紹貝葉斯Meta分析與WinBUGS軟件。

1 起源與發(fā)展

英國數(shù)學(xué)家Bayes T[3]于1763年在《論有關(guān)機(jī)遇問題的求解》中提出了貝葉斯公式和一種歸納推理的理論（但可能因其認(rèn)為該理論尚存在不完善的地方，在其生前并未發(fā)表），后被Laplace PC等一些統(tǒng)計(jì)學(xué)者發(fā)展為一種系統(tǒng)的統(tǒng)計(jì)推斷方法，稱為貝葉斯方法。

在20世紀(jì)80年代之前貝葉斯統(tǒng)計(jì)因高維函數(shù)的積分，無法給出恰當(dāng)?shù)慕馕鼋?，求解這些積分成為其發(fā)展的障礙，因此一直停留在理論階段。20世紀(jì)90年代起MCMC方法廣泛應(yīng)用于貝葉斯統(tǒng)計(jì)，成功地解決了限制貝葉斯統(tǒng)計(jì)發(fā)展的高維積分運(yùn)算問題，為貝葉斯統(tǒng)計(jì)帶來了革命性的突破，進(jìn)而使其在更多的領(lǐng)域得到應(yīng)用。在醫(yī)學(xué)領(lǐng)域，貝葉斯方法廣泛應(yīng)用于不同的數(shù)據(jù)類型統(tǒng)計(jì)分析中，如遺傳數(shù)據(jù)、縱向數(shù)據(jù)、生存數(shù)據(jù)及缺少數(shù)據(jù)等；同時也應(yīng)用于不同的研究方法，如臨床試驗(yàn)實(shí)施[4]循證醫(yī)學(xué)[5]等。

2 貝葉斯統(tǒng)計(jì)基礎(chǔ)

2.1 基本思想 貝葉斯統(tǒng)計(jì)學(xué)派認(rèn)為在觀察到樣本之前，對于任一未知的參數(shù)θ一般有一定的了解，即已經(jīng)積累一些關(guān)于參數(shù)θ的信息——“先驗(yàn)信息”，在對未知參數(shù)進(jìn)行統(tǒng)計(jì)推斷時應(yīng)綜合先驗(yàn)信息，也應(yīng)考慮樣本信息。用統(tǒng)計(jì)學(xué)語言可描述為：θ作為一個隨機(jī)變量，有一定的先驗(yàn)分布，其分布密度為θ～π(θ)。在獲得樣本之后（給定的樣本信息），θ的后驗(yàn)分布π(θ/x)應(yīng)包含θ的綜合信息，關(guān)于參數(shù)θ的統(tǒng)計(jì)推斷均基于θ的后驗(yàn)分布進(jìn)行。因此，貝葉斯統(tǒng)計(jì)方法的關(guān)鍵在于所作出的任何推斷都只須根據(jù)后驗(yàn)分布π(θ/x)，而不再涉及樣本x的分布。

2.2 貝葉斯公式 貝葉斯統(tǒng)計(jì)學(xué)的基礎(chǔ)是貝葉斯公式和貝葉斯定理。貝葉斯公式是基于條件概率的定義及全概率公式推導(dǎo)而得的，因此是貝葉斯公式的事件形式，如下：

設(shè)試驗(yàn)E的樣本空間為S，A為E的事件B1,B2,...,Bn為樣本空間S的一個劃分，且P(A)＞0，P(Bi)＞0 （i=1,2,...,n），則由條件概率的定義及全概率公式可得：

貝葉斯公式的密度函數(shù)形式如下：

設(shè)x=(x1,x2,..,xn)是來自某總體的一個樣本，該總體的概率密度函數(shù)為p(x/θ)，當(dāng)給定一組觀察值x=(x1,x2,..,xn) ，θ的條件概率分布為：

即在樣本x=(x1,x2,..,xn)下θ的后驗(yàn)分布。其中π(θ)為參數(shù)θ的先驗(yàn)分布；

為樣本 x=(x1,x2,..,xn)的聯(lián)合條件密度函數(shù)，也即似然函數(shù)；

為x的邊緣密度函數(shù)，是一個與θ無關(guān)的量。

2.3 先驗(yàn)分布 貝葉斯統(tǒng)計(jì)分析中最重要且最受經(jīng)典統(tǒng)計(jì)學(xué)派批判的一點(diǎn)即先驗(yàn)分布的選取。統(tǒng)計(jì)學(xué)家提出多種方法，但至今理論上仍無統(tǒng)一的、完整的、不失一般性的確定先驗(yàn)分布的方法，但是在實(shí)用的范圍內(nèi)，常見的問題所采用的先驗(yàn)分布，已經(jīng)得到正確的驗(yàn)證[6]，這里僅對眾多方法做一簡述，有興趣的讀者可參考相關(guān)書籍[7]。

常用的選取先驗(yàn)分布的方法包括三種：

①根據(jù)貝葉斯假設(shè)選擇先驗(yàn)分布：貝葉斯假設(shè)是指參數(shù)θ的無信息先驗(yàn)分布π(θ)應(yīng)在θ的取值范圍內(nèi)是“均勻”分布的，在貝葉斯假設(shè)下，似然函數(shù)l(θ/x) 即為后驗(yàn)密度的核。

②共軛先驗(yàn)分布：設(shè)θ是總體分布中的參數(shù)，π(θ)是θ的先驗(yàn)密度函數(shù)，假如由抽樣信息算得的后驗(yàn)密度函數(shù)與π(θ)有相同的函數(shù)形式，則稱為θ的共軛先驗(yàn)分布，也即π(θ)與π(θ/x)屬于同一類分布族。共軛先驗(yàn)分布有兩個優(yōu)點(diǎn)：一是計(jì)算簡便；二是后驗(yàn)分布中的一些參數(shù)可以得到很好的解釋。雖然共軛先驗(yàn)分布計(jì)算簡便，但應(yīng)以合理性為首要原則。常用的共軛先驗(yàn)分布見表1。

表1 常用共軛先驗(yàn)分布

③采用Jeffreys原則確定無信息先驗(yàn)分布：Jeffreys對先驗(yàn)分布的確定做出重大貢獻(xiàn)，利用Fisher信息矩陣給出了確定無信息先驗(yàn)分布的一般方法。Jeffreys原則包括兩部分：一是對先驗(yàn)分布應(yīng)有一個合理的要求；二是給出一個具體的方法去求的合乎要求的先驗(yàn)分布。Jeffreys利用了Fisher信息矩陣的一個不變性質(zhì)，發(fā)現(xiàn)θ的無信息先驗(yàn)分布應(yīng)以信息陣I(θ)的行列式的平方根為核。

3 貝葉斯與傳統(tǒng)Meta分析

傳統(tǒng)的Meta分析是基于經(jīng)典的頻率學(xué)派統(tǒng)計(jì)理論，在固定效應(yīng)模型或隨機(jī)效應(yīng)模型的選擇是基于統(tǒng)計(jì)量Q進(jìn)行。經(jīng)典的隨機(jī)效應(yīng)模型中可以通過Q檢驗(yàn)獲得研究間方差的矩估計(jì)，但是較難獲得其95%可信區(qū)間（CI），且檢驗(yàn)功效較低。再者，傳統(tǒng)的Meta分析很難應(yīng)對復(fù)雜的模型，如對于二分類資料的Meta分析選擇效應(yīng)量為OR，經(jīng)典方法是將OR對數(shù)化，然后假設(shè)log（OR），并服從正態(tài)分布，相應(yīng)的計(jì)算都是基于正態(tài)分布假設(shè)的前提下，因此在存在小樣本資料時，如果不符合正態(tài)近似的條件，經(jīng)典方法不能處理。此外當(dāng)納入的研究存在較多的極端值時，經(jīng)典方法很難識別隨機(jī)效應(yīng)[8]。

貝葉斯統(tǒng)計(jì)將參數(shù)θ作為一個隨機(jī)變量，有一定的先驗(yàn)分布，在獲得樣本之后（給定的樣本信息），θ的后驗(yàn)分布π(θ/x）應(yīng)包含θ的綜合信息，可從θ的后驗(yàn)分布獲得參數(shù)θ的貝葉斯估計(jì)。因此，可很好的解決經(jīng)典Meta分析存在的缺陷。Carlin[9]研究了2×2四格表的貝葉斯Meta分析，采用可交換先驗(yàn)分布，用隨機(jī)模擬的方法得到參數(shù)的后驗(yàn)分布。Warn和Thompson等給出了二分類變量絕對風(fēng)險(xiǎn)差（ARD）及相對危險(xiǎn)度（RR）隨機(jī)效應(yīng)模型貝葉斯估計(jì)方法[8]。

總的來說，兩者只是在行Meta分析統(tǒng)計(jì)方法上的不同，在Meta分析的其他步驟及報(bào)告規(guī)范上仍是相同的。

4 模型的構(gòu)建及WinBUGS代碼

本文以二分類資料為例簡單介紹貝葉斯統(tǒng)計(jì)在Meta分析中應(yīng)用及實(shí)現(xiàn)，數(shù)據(jù)源于Colditz GA等《Efficacy of BCG vaccine in the prevention of tuberculosis》[10]一文。

對于二分類資料模型構(gòu)建可以基于率的logit轉(zhuǎn)換，實(shí)驗(yàn)組及對照組的率的logit轉(zhuǎn)化的差值服從正態(tài)分布。設(shè)共納入n個研究，nt、rt、nc、rc分別為治療組和對照組總?cè)藬?shù)及事件發(fā)生例數(shù)；設(shè)治療組事件數(shù)rt，對照組事件發(fā)生數(shù)為rc，則rt、rc均服從二項(xiàng)分布，即：

其中pti、pci分別為治療組和對照組事件發(fā)生率。率的logit轉(zhuǎn)換為：

故該方法可以按照率的logit轉(zhuǎn)換的差值來進(jìn)行估計(jì)，相應(yīng)固定效應(yīng)模型及隨機(jī)效應(yīng)模型如下：

?

5 WinBUGS軟件

5.1 簡介 BUGS（Bayesian inference using gibbs sampling）軟件最初于1989年由位于英國劍橋的生物統(tǒng)計(jì)學(xué)研究所（Biostatistics the Medical Research Council, Cambridge,United Kingdom）研制的，現(xiàn)在由這個研究所和位于倫敦的S.t Mary＇s皇家學(xué)院醫(yī)學(xué)分院（the Imperial College School of Medicine）共同開發(fā)。BUGS的運(yùn)行以MCMC方法為基礎(chǔ)，它將所有未知參數(shù)都看做隨機(jī)變量，然后對此種類型的概率模型進(jìn)行求解。它所使用的編程語言非常容易理解，允許使用者直接對研究的概率模型作出說明。WinBUGS是在BUGS基礎(chǔ)上開發(fā)面向?qū)ο蠼换ナ降腤indows軟件版本，提供了圖形界面，允許用鼠標(biāo)的點(diǎn)擊直接建立研究模型。WinBUGS使用MCMC技術(shù)從復(fù)雜模型的后驗(yàn)分布中產(chǎn)生樣本，提供了一個有效的方法估計(jì)貝葉斯模型。

WinBUGS可以從http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/WinBUGS14.exe免費(fèi)下載，目前最新的版本為WinBUGS 1.4.3，下載完成后，雙擊WinBUGS14.exe，按提示進(jìn)行安裝。安裝完成后，需要進(jìn)行注冊，否則無法使用軟件的全部功能，注冊碼可從http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/WinBUGS14_immortality_key.txt下載。

5.2 分析流程 WinBUGS分析的基本流程如圖1所示，主要分為5大步。

圖1 WinBUGS操作一般流程圖

第一步，構(gòu)建模型及數(shù)據(jù)輸入（以固定效應(yīng)模型為例），如圖2；

圖2 模型的建立及數(shù)據(jù)輸入

第二步，模型的確定：選擇M o d e l下拉菜單中Specification選項(xiàng)，會跳出Specification Tool對話框（圖3）；

圖3 選擇Model對話框

模型的確定分以下4步：

①模型的檢查（Check model）：將光標(biāo)移到描述統(tǒng)計(jì)模型的語句model處，選中model字樣，再點(diǎn)擊Specification Tool對話框的check model處，若對模型描述的語法正確的話，則窗口底部左下角會提示 Model is syntactically correct（圖4）。

圖4 檢查模型

②加載數(shù)據(jù)（Load data）：將光標(biāo)移到數(shù)據(jù)的語句前面的list處，選中l(wèi)ist字樣，再點(diǎn)擊Specification Tool對話框的load data，若對模型描述的語法正確的話，則窗口底部左下角會提示data loaded（圖5）。

圖5 加載數(shù)據(jù)

③編譯（Compile）：點(diǎn)擊Specification Tool對話框的compile，編譯成功后，窗口底部左下角會提示model compiled，同時也激活了初始值的按鈕（圖6）

圖6 模型的編譯

④加載初始值（Load inits）：與加載數(shù)據(jù)類似，選定初始值的list，點(diǎn)擊Specification Tool對話框的load inits，加載成功后，窗口底部左下角會提示model is initialized（圖7）。

圖7 加載初始值，成功建立模型

在定義模型的過程中，可根據(jù)實(shí)際定義Markov鏈數(shù)，在定義一條以上鏈時，當(dāng)各條鏈都收斂到差不多的后驗(yàn)量時，可以認(rèn)為迭代收斂，結(jié)果是正確的。有時在定義多條鏈時或者有遺漏的初始值，可以點(diǎn)擊gen inits來生成初始值。

第三步，指定要考察的參數(shù)：從Inference下拉菜單中選中Samples選項(xiàng)，出現(xiàn)Sample Monitor Tool對話框（圖8）；

圖8 參數(shù)設(shè)定

于node后的輸入框中依次輸出delta，而后點(diǎn)擊set按鈕。beg處可以輸入相應(yīng)數(shù)字，如輸入1000則表示前1000次退火用以消除初始值的影響，從1001次后開始抽樣。本例題只是演示故未設(shè)置。其他參數(shù)采用默認(rèn)值即可，也可根據(jù)實(shí)際需要設(shè)定（在次不作詳述）。

第四步，迭代運(yùn)算：判斷迭代運(yùn)算的收斂性，可從總運(yùn)算結(jié)果中判斷，因此再次直接正式迭代運(yùn)算。選擇 Model下拉菜單中的Updatd選項(xiàng)，會彈出Update Tool對話框，在該對話框的updates處寫入所需迭代的次數(shù)，默認(rèn)為1000次，本例指定迭代10000次，然后點(diǎn)擊update按鈕（圖9）。

圖9 模型迭代設(shè)置

第五步，輸出迭代計(jì)算結(jié)果：從Inference下拉菜單中選中Samples選項(xiàng)，出現(xiàn)Sample Monitor Tool對話框，然后于node后的輸入框中輸入“*”（“*”代表指定的所有參量）。可以獲得相應(yīng)的后驗(yàn)分布的相關(guān)統(tǒng)計(jì)量以及迭代是否收斂，如點(diǎn)擊trace給出Gibbs動態(tài)抽樣圖，點(diǎn)擊stats會給出參數(shù)的計(jì)算結(jié)果（圖10）。

圖10 抽樣結(jié)果

5.3 計(jì)算結(jié)果 于WinBUGS軟件中迭代10000次，前1000次退火用于初始值的影響，得到參數(shù)的后驗(yàn)分布估計(jì)如（表2）。由表2可知：固定效應(yīng)模型計(jì)算結(jié)果OR=0.6198，95%CI為（0.5714～0.6721）；隨機(jī)效應(yīng)模型的OR=0.4776，95%CI為（0.3044～0.6987）；由上述結(jié)果可知，固定效應(yīng)模型與隨機(jī)效應(yīng)模型結(jié)果相差較多，究其原因是固定效應(yīng)模型假設(shè)所納入的研究均來自同一總體，也即，此條件并不符合數(shù)據(jù)的真實(shí)情況；而隨機(jī)效應(yīng)模型假設(shè)所納入的研究來自不同的總體，即，因此隨機(jī)效應(yīng)模型結(jié)論更為可靠。雖然貝葉斯Meta統(tǒng)計(jì)并無I2等經(jīng)典Meta分析模型選擇的統(tǒng)計(jì)量，但基本思想與經(jīng)典統(tǒng)計(jì)學(xué)派的思想是一致的，即模型的選擇依賴于納入的研究是否來自同一總體。參數(shù)的Gibbs動態(tài)圖（圖11）提示迭代基本收斂（圖12）。

表2 計(jì)算結(jié)果

圖11 后驗(yàn)分布核密度估計(jì)圖及Gibbs抽樣動態(tài)蹤跡圖

圖12 迭代歷史圖

6 結(jié)語

傳統(tǒng)的Meta分析通過Q檢驗(yàn)獲得研究間方差的矩估計(jì)，并以此作為固定效應(yīng)模型或隨機(jī)效應(yīng)模型的選擇的依據(jù)，其計(jì)算都是基于正態(tài)分布假設(shè)的前提下，因此在存在小樣本資料時，如果不符合正態(tài)近似的條件，經(jīng)典方法不能處理。貝葉斯統(tǒng)計(jì)方法不受限于經(jīng)典統(tǒng)計(jì)學(xué)派的前提假設(shè)，可以結(jié)合先驗(yàn)信息、樣本信息及總體信息，獲得后驗(yàn)分布并在其基礎(chǔ)上方便地得到效應(yīng)量的合并值、研究間方差等參數(shù)及95%CI。相信貝葉斯統(tǒng)計(jì)方法在循證醫(yī)學(xué)/Meta分析中會有更廣泛的引用。

[1] 曾憲濤,Joey S.W. Kwong,田國祥,等. Meta分析系列之二：Meta分析的軟件[J]. 中國循證心血管醫(yī)學(xué)雜志,2012,4(2):89-91.

[2] 曾憲濤,曹世義,孫鳳,等. Meta分析系列之六：間接比較及網(wǎng)狀Meta分析[J]. 中國循證心血管醫(yī)學(xué)雜志,2012,4(5):399-402

[3] Bayes T. An essay towards solving a problem in the doctrine of chances [J]. Phil Trans,1763,53:370-418.

[4] Spiegelhalter DJ,Freedman LS,Parmar MKB. Bayesian Approaches to Randomized Trials [J]. J R Statist Soc,1994,157(3):357-416.

[5] Smith TC, Spiegelhalter DJ, Thomas A. Bayesian approaches to random-effects meta-analysis: a comparative study [J]. Stat Med,1995,14(24):2685-99.

[6] 張堯庭,陳漢峰. 貝葉斯統(tǒng)計(jì)推斷[M]. 北京:科學(xué)出版社,1991.

[7] 茆詩松. 貝葉斯統(tǒng)計(jì)[M]. 北京:中國統(tǒng)計(jì)出版社,1999.

[8] Warn DE,Thompson SG,Spiegelhalter DJ. Bayesian random effects meta-analysis of trials with binary outcomes: methods for the absolute risk difference and relative risk scales [J]. Stat Med,2002,21(11):1601-23.

[9] Carlin JB. Meta-analysis for 2×2 tables: a Bayesian approach [J].Stat Med,1992,11(2):141-58.

[10] Colditz GA,Brewer TF,Berkey CS,et al. Efficacy of BCG vaccine in the prevention of tuberculosis. Meta-analysis of the published literature [J]. JAMA,1994,271(9):698-702.