如何用stata軟件實現(xiàn)貝葉斯meta分析

張晟周潔何書汪徐林沈毅

南通大學公共衛(wèi)生學院流行病與衛(wèi)生統(tǒng)計學系(226019)

系統(tǒng)評價與meta分析是整合研究結果并進行定量分析的方法，是循證醫(yī)學不可或缺的部分[1]。貝葉斯統(tǒng)計以參數(shù)θ為有先驗分布的隨機變量，在獲得樣本之后，依據(jù)貝葉斯定理，得到θ的后驗分布π(θ|x) ，并根據(jù)其后驗分布獲得參數(shù)θ[2]。20世紀90年代起馬爾可夫鏈蒙特卡羅(Markov chain Monte Carlo,MCMC)算法成功解決了限制其發(fā)展的高維積分運算問題，為貝葉斯統(tǒng)計帶來了革命性突破[3]。

對meta分析而言，貝葉斯方法在處理復雜的隨機效應、分層結構或稀疏數(shù)據(jù)時比經(jīng)典的頻率學派更具吸引力[4-6]。特別是納入文獻的樣本量較小時，正態(tài)性檢驗將無法進行；但不滿足正態(tài)性，經(jīng)典的頻率學派模型合并估計效應時就有可能存在一定誤差[7]。另外，如果納入研究的數(shù)據(jù)中存在較多極端值，經(jīng)典方法很難識別其隨機效應。貝葉斯統(tǒng)計則可通過指定先驗分布和超先驗分布，很好地解決經(jīng)典meta分析存在的缺陷，不僅提高模型估計的精度，而且建模方式更為靈活。Rebecca等人指出，如果納入的文獻量較少，或者文獻中數(shù)據(jù)較為稀疏時，貝葉斯方法可對組間方差選擇合適的信息先驗(informative prior)，獲得更為精確的合并效應值[6]。

近年來，盡管模擬算法不斷改進，但貝葉斯方法的計算復雜性仍使眾多研究者望而卻步。此外，數(shù)據(jù)處理軟件的匱乏也成為貝葉斯統(tǒng)計發(fā)展的桎梏。最常用于貝葉斯分析的BUGS類軟件是一種為Windows系統(tǒng)編寫的免費軟件，內(nèi)含一系列抽樣方法，可對給定問題自動挑選最佳解決方案。但該軟件編程復雜，且無法直接提供可視化meta分析結果，如森林圖、漏斗圖、模型診斷圖等，而這一點卻是meta分析所必不可少的。

2015年更新的stata軟件14.0版本新增貝葉斯meta分析方法，并與該軟件所特有的可視化圖形展示相結合，從而構造了一個易學易用的貝葉斯meta分析平臺[8-9]。本文主要介紹如何使用stata 14.0版本實現(xiàn)貝葉斯meta分析。

廣義貝葉斯隨機效應模型的meta分析

廣義線性模型(generalized linear models,GLM)作為一般線性模型的推廣，將諸多不同的分布函數(shù)統(tǒng)一到指數(shù)族框架內(nèi)。無論因變量y是連續(xù)性還是離散性隨機變量，均可表示為指數(shù)族的概率分布形式。在廣義貝葉斯隨機效應模型中，首先，假定每個研究效應量的估計值yi(i=1,…,k)服從均數(shù)為θi，方差為si2的正態(tài)分布：

(1)

θi～N(θ,τ2)

(2)

(3)

或者等價于：

式中ui為模型的隨機部分，而τ2代表研究間的變異程度。在貝葉斯meta中，對τ2可以選擇log-normal，inverse-gamma或gamma分布。對于ui往往選擇無信息先驗，以使其結果更為客觀。

對于參數(shù)后驗分布的估計有三種方法：①直接利用MCMC獲得對θ和τ2聯(lián)合后驗分布的參數(shù)估計；②通過數(shù)值積分法計算θ和τ2聯(lián)合后驗分布的矩估計和百分位數(shù)；③運用重要性抽樣法推導θ和τ2的聯(lián)合后驗分布[6]。

而stata 14.0是采用第一種方法，即MCMC完成貝葉斯分析，其主要命令為bayesmh。該命令使用調(diào)整后的Metropolis-Hastings(MH)和倒伽馬(inverse-gamma) 方法來計算一系列的貝葉斯回歸模型。同時，stata 14.0也支持Gibbs抽樣以選擇似然函數(shù)及合并先驗值。另外，它還提供了一系列檢查MCMC收斂性和效率，總結參數(shù)及參數(shù)函數(shù)的后驗信息，假設檢驗以及模型比較的命令。

stata實現(xiàn)貝葉斯meta分析

1.二分類數(shù)據(jù)

二分類數(shù)據(jù)來源于一篇關于人工肝支持系統(tǒng)(ALSS)與標準藥物治療對慢加急性肝衰竭死亡率影響的meta分析，該研究共納入了9篇文獻，包含6個隨機對照試驗和3個隊列研究[10]。數(shù)據(jù)整理格式如表1。

stata在做貝葉斯meta分析之前，首先需通過metan命令對數(shù)據(jù)進行預處理：

generate live0=total0-death0

(4)

generate live1=total0-death1

(5)

metan death0 liver0 death1 liver1,randomi or

(6)

表1　二分類變量數(shù)據(jù)分布格式(1)

表2　二分類變量數(shù)據(jù)分布格式(2)

其中，death0和live0分別指接受ALSS治療后死亡和生存的病人數(shù)，death1 和live1指接受藥物治療后死亡和生存的病人數(shù)。live0和live1可以使用generate命令用total減去death生成(6)，該命令運行后會自動產(chǎn)生比值比(odds Ratio,OR)及其對數(shù)值的標準誤，系統(tǒng)自動將比值比命名為_ES，對數(shù)值的標準誤為_selogES。為方便后續(xù)分析，可對其進行變量變換：D=log(_ES)、var=[(_selogES)]2。

在stata中，i.study通常默認為設置啞變量，但在貝葉斯meta中，i.study則表示納入研究的序號，所以需先通過命令fvset進行設置，以示區(qū)分。

fvset base none study

(7)

貝葉斯計算依賴于隨機MCMC方法，為保證結果的可重復，應先設置隨機種子數(shù)，本例種子數(shù)設為14。接下來對均數(shù)j5i0abt0b和方差{sig2}設置先驗分布，先驗可使用Gibbs抽樣產(chǎn)生。為了提高效率，各參數(shù)可安排在不同模塊(block)下運行。運行前，先設置退火次數(shù)burnin，本例為5000次；迭代次數(shù)mcmcsize則設置為50000次；由于設置的迭代次數(shù)較多，所需時間較長，可通過adaptation選項設置更新次數(shù)和模擬的監(jiān)測過程。本例中，監(jiān)測burnin更新數(shù)設為10(every(10))，表示退火5000/10=500次時，則在結果上展示一個點；迭代更新數(shù)設為2500(every(2500))，每迭代2500次則顯示迭代次數(shù)。

set seed 14

(8)

bayesmh D i.study,likelihood(normal(var))noconstant///

prior({D:i.study},normal(j5i0abt0b,{sig2}))///

prior(j5i0abt0b,normal(0,1000))///

prior({sig2),igamma(0.001,0.001))///

block({D:i,study},split)///

block({sig2},gibbs)///

block(j5i0abt0b,gibbs)///

burn(5000)adaptation(every(10))dots(500,every(2500))mcmcsize(50000)

(9)

通過該命令可獲得后驗均數(shù)d和方差sig2，以及它們的95%置信區(qū)間(credit intervals,CrIs)。如表3所示，本例d為-0.405 (95%CrI:-0.747～-0.135)，sig2為0.064 (95%CrI:0.001～0.400)。計算d的反對數(shù)值(OR)為0.67，與原文中傳統(tǒng)meta分析的結果0.69十分相近。CrIs與可信區(qū)間(confidence intervals,CIs)作用相似，但其對結果解釋性更直接和相關。

表3　二分類數(shù)據(jù)貝葉斯meta分析后驗均數(shù)與方差結果及其95%置信區(qū)間

進行貝葉斯meta分析時，命令運行的效率與預期樣本量(ESS)和自相關性(Corr.time)有關，可用bayesstatsess命令檢查運行效率的相關參數(shù)：

bayesstats ess j5i0abt0b{sig2}

(10)

結果給出Efficiency，ESS和Corr.Time三個評價指標(表4)。ESS的大小是根據(jù)MCMC方法計算所得，其值越接近真實樣本量，模擬結果越好。Corr.Time是Efficiency的倒數(shù)，表示時間序列中某指標過去與將來值的相關性，有時稱為“滯后相關性”或者“序列相關性”，也可以理解為一系列數(shù)字在時間上的相關性，它常用于檢查數(shù)據(jù)的隨機性，其值越低表示效率越高。

貝葉斯的收斂性可通過圖形化展示。命令如下：

bayesgraph diagnosticsj5i0abt0b{sig2}

(11)

表4　貝葉斯meta運行效率參數(shù)表

圖1　貝葉斯收斂診斷圖

其結果包含后驗均數(shù)和方差的軌跡圖，自相關性圖，直方圖以及核密度曲線圖四個部分。圖1軌跡圖(trace)橫軸顯示迭代的次數(shù)，縱軸顯示參數(shù)的后驗分布值，當MCMC達到穩(wěn)態(tài)時，參數(shù)的模擬值會在均值附近以相同的幅度上下波動。自相關性圖(Autocorrelation)橫軸表示滯后時間，縱軸仍為后驗分布值。該圖顯示了MCMC樣本自相關性從滯后0開始的一系列滯后范圍，在滯后0這一點上的值等于MCMC的樣本方差。還有兩個圖形分別為直方圖和核密度曲線圖，主要展示了參數(shù)的邊際后驗分布，是總結MCMC輸出的各種參數(shù)統(tǒng)計量的補充。后驗均數(shù)的單峰直方圖提示j5i0abt0b近似正態(tài)分布，與指定的共軛正態(tài)的邊際后驗正態(tài)分布相一致；{var}的直方圖是單峰右偏態(tài)，與指定先驗的邊際后驗分布為倒伽馬分布相一致。核密度曲線圖是模擬邊際后驗分布的另一種方式，本例j5i0abt0b和{var}的核密度曲線圖與直方圖的結果相一致。

上述貝葉斯meta分析方法是基于預處理后的數(shù)據(jù)進行的。此外，stata14.0還提供了一套更為簡便的方法，可直接對陽性率建模而無需轉化為OR值，但其數(shù)據(jù)格式(表2)及主要命令與上述方法略有不同：

generate mu=study

(12)

fvset base none mu study

(13)

set seed 14

(14)

bayesmh deaths i.mu 1.treat#i.study,///

likelihood(binlogit(total))noconstant prior({deaths:i.mu},flat)///

prior({deaths:1.treat#i.study},normal(j5i0abt0b,{sig2}))///

prior(j5i0abt0b,normal(0,1000))prior({sig2},igamma(0.001,0.001))///

block({deaths:1.treat#i.study},split)///

block({deaths:i.mu},split)block(j5i0abt0b,gibbs)///

block({sig2},gibbs)///

burnin(5000)adaptation(every(10))dots(500,every(2500))mcmcsize(50000)

(15)

連續(xù)性數(shù)據(jù)

連續(xù)性數(shù)據(jù)資料來源于一篇比較噻托溴銨與三聯(lián)療法對于慢性阻塞性肺病治療效果的meta分析，共納入6篇文獻，皆為隨機對照試驗[11]。原始數(shù)據(jù)如表5。和二分類數(shù)據(jù)一樣，做貝葉斯分析之前，需對數(shù)據(jù)進行預處理，命令如下：

metan Total1 Mean1 SD1 Total2 Mean2 SD2,randomi nostandard

(16)

表5　連續(xù)性變量數(shù)據(jù)數(shù)據(jù)分布格式

該命令產(chǎn)生兩個變量：均數(shù)差(_ES)及均數(shù)差的標準誤(_seES)。用D直接代表_ES，var代表(_seES)2(注意在二分類數(shù)據(jù)中D是_ES的對數(shù)值，這是分類結局與連續(xù)性結局細微的差別)。接著運行(7)-(11)的代碼，就可以完成貝葉斯meta分析。由表6可見均數(shù)差結果為67.90，與原文結果63.68相似。

表6　連續(xù)性數(shù)據(jù)貝葉斯meta分析均數(shù)差與方差及其95%置信區(qū)間

討　　論

本研究通過兩個經(jīng)典的例子展示了如何使用stata實現(xiàn)貝葉斯隨機效應模型meta分析，例1是二分類資料，計算合并OR值，例2是連續(xù)性變量，計算合并均數(shù)差。這兩個例子在原文中均采用傳統(tǒng)meta分析的方法，貝葉斯方法所得結果均與原文獻相似。

貝葉斯分析是用于統(tǒng)計建模，結果解釋和數(shù)據(jù)預測的強大分析工具，它遵循著貝葉斯定理，用一個簡單的公式將先驗信息與來自現(xiàn)有數(shù)據(jù)的證據(jù)相結合，形成模型參數(shù)的后驗分布[12]。貝葉斯meta與頻率派方法的不同之處在于它認為數(shù)據(jù)和模型參數(shù)都是隨機量，似然函數(shù)是基于給定數(shù)據(jù)的合理性而進行設定的[13]。貝葉斯meta分析可提供治療效果的概率解釋，或者是效應大于(或小于)一個特定值的概率，也可以評價療效是否大于最小臨床意義變化值(MCID)[14]。

與頻率學派相比，貝葉斯meta分析優(yōu)勢如下：貝葉斯分析可用于所有參數(shù)模型，適用性較為廣泛；貝葉斯分析對先驗信息的提取兼顧理論合理和原則可行，從而為特定問題獲得更精準結果；先驗信息的納入可降低小樣本量的影響，從而部分解決結構零的問題。

貝葉斯meta分析也有一些自身的局限：不同的先驗分布有可能得到不同結果，如何選擇合適的先驗一直是貝葉斯方法使用的絆腳石；復雜的計算過程可能消耗更多的迭代時間。

stata軟件是一套提供數(shù)據(jù)管理、統(tǒng)計分析和圖表繪制的整合性統(tǒng)計軟件，在社會學、經(jīng)濟學和生物統(tǒng)計學中使用廣泛。過去，stata進行貝葉斯分析時只能通過連接命令借助外部軟件(如BUGS，JAGS或MLwiN)完成[15]。從2015年開始用戶可在stata 14.0版本中直接進行貝葉斯分析，它可分步提供詳細的結果和理想的圖形。對那些希望實現(xiàn)貝葉斯meta分析但對BUGS了解不多的研究者，stata14.0應該是更好的選擇。