張晟周潔何書(shū)汪徐林沈毅

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

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

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

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

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

廣義貝葉斯隨機(jī)效應(yīng)模型的meta分析

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

(1)

θi～N(θ,τ2)

(2)

(3)

或者等價(jià)于：

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

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

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

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

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

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

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

generate live0=total0-death0

(4)

generate live1=total0-death1

(5)

metan death0 liver0 death1 liver1,randomi or

(6)

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

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

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

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

fvset base none study

(7)

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

set seed 14

(8)

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

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

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

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

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

block({sig2},gibbs)///

block(0qee4ec,gibbs)///

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

(9)

通過(guò)該命令可獲得后驗(yàn)均數(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)。計(jì)算d的反對(duì)數(shù)值(OR)為0.67，與原文中傳統(tǒng)meta分析的結(jié)果0.69十分相近。CrIs與可信區(qū)間(confidence intervals,CIs)作用相似，但其對(duì)結(jié)果解釋性更直接和相關(guān)。

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

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

bayesstats ess o0cqmiu{sig2}

(10)

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

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

bayesgraph diagnosticsmawkukm{sig2}

(11)

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

圖1　貝葉斯收斂診斷圖

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

上述貝葉斯meta分析方法是基于預(yù)處理后的數(shù)據(jù)進(jìn)行的。此外，stata14.0還提供了一套更為簡(jiǎn)便的方法，可直接對(duì)陽(yáng)性率建模而無(wú)需轉(zhuǎn)化為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(q02mmos,{sig2}))///

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

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

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

block({sig2},gibbs)///

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

(15)

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

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

metan Total1 Mean1 SD1 Total2 Mean2 SD2,randomi nostandard

(16)

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

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

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

討　　論

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

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

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

貝葉斯meta分析也有一些自身的局限：不同的先驗(yàn)分布有可能得到不同結(jié)果，如何選擇合適的先驗(yàn)一直是貝葉斯方法使用的絆腳石；復(fù)雜的計(jì)算過(guò)程可能消耗更多的迭代時(shí)間。

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