如何用stata軟件實現貝葉斯meta分析

張晟周潔何書汪徐林沈毅

南通大學公共衛生學院流行病與衛生統計學系(226019)

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

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

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

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

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

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

(1)

θi～N(θ,τ2)

(2)

(3)

或者等價于：

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

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

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

stata實現貝葉斯meta分析

1.二分類數據

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

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

generate live0=total0-death0

(4)

generate live1=total0-death1

(5)

metan death0 liver0 death1 liver1,randomi or

(6)

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

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

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

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

fvset base none study

(7)

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

set seed 14

(8)

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

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

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

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

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

block({sig2},gibbs)///

block(g0gggggg,gibbs)///

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

(9)

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

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

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

bayesstats ess g0gggggg{sig2}

(10)

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

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

bayesgraph diagnosticsg0gggggg{sig2}

(11)

表4　貝葉斯meta運行效率參數表

圖1　貝葉斯收斂診斷圖

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

上述貝葉斯meta分析方法是基于預處理后的數據進行的。此外，stata14.0還提供了一套更為簡便的方法，可直接對陽性率建模而無需轉化為OR值，但其數據格式(表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(g0gggggg,{sig2}))///

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

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

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

block({sig2},gibbs)///

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

(15)

連續性數據

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

metan Total1 Mean1 SD1 Total2 Mean2 SD2,randomi nostandard

(16)

表5　連續性變量數據數據分布格式

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

表6　連續性數據貝葉斯meta分析均數差與方差及其95%置信區間

討　　論

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

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

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

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

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