R 軟件 bmeta 程序包是一款通過調用 JAGS 軟件來實現貝葉斯 Meta 分析和 Meta 回歸的程序包，該程序基于“馬爾可夫鏈-蒙特卡羅”（MCMC）算法來合并不同類型資料（二分類、連續和計數）的各種效應量（OR、MD 和 IRR）。該程序包具有命令函數參數少、提供模型豐富、繪圖功能強大、易于理解和掌握等優點。本文將結合實例介紹展示 bmeta 程序包實現貝葉斯 Meta 分析與 Meta 回歸的完整操作流程。
引用本文： 石豐豪, 孟蕊, 芮明軍, 馬愛霞. 應用 R 軟件 bmeta 程序包實現貝葉斯 Meta 分析與 Meta 回歸. 中國循證醫學雜志, 2020, 20(12): 1482-1488. doi: 10.7507/1672-2531.202004178 復制
Meta 分析作為一種整合單個研究效應量進行證據合并的常用統計方法，在循證醫學中占有重要地位[1]。貝葉斯 Meta 分析是基于貝葉斯統計發展起來的一種的 Meta 分析方法，主要采用“馬爾科夫鏈—蒙特卡羅”（Markov chain Monte Carlo，MCMC）方法，因其在處理復雜隨機效應、分層結構或是稀疏數據時比頻率學 Meta 分析方法更有優勢，目前越來越受歡迎。董圣杰等[2]介紹了貝葉斯 Meta 分析的建模過程及如何在 WinBUGS 軟件實現貝葉斯 Meta 分析。但是 BUGS 編程復雜，操作過程繁瑣并且無法繪制森林圖、漏斗圖等 Meta 分析必不可少的圖形。R 軟件作為一款開源且繪圖功能強大的軟件能夠實現幾乎所有的 Meta 分析類型[3, 4]。然而目前 R 軟件基于貝葉斯方法計算的程序包更多是用來進行網狀 Meta 分析[5-7]，用來進行傳統頭對頭 Meta 分析的程序包的理論原理均基于經典的頻率學派[8, 9]。在這樣的背景需求下，bmeta 程序包被研發出來，該程序包采用簡潔的代碼輸入的方式，并通過調用 JAGS 來實行貝葉斯 Meta 分析和 Meta 回歸。本文以《Efficacy of BCG vaccine in the prevention of tuberculosis meta-analysis of the published literature》[10]一文的數據為例來展示該程序包的具體操作。
1 軟件安裝及程序包加載
1.1 R 軟件安裝
截至完稿時，R 軟件最新版本為 R-4.0.0，用戶可在下面的網址根據自己計算機的操作系統下載最新版本：https://www.r-project.org/。
1.2 JAGS 軟件的下載
JAGS（Just another Gibbs sampler）軟件是一款基于貝葉斯理論運算的編程軟件，但尚不完善的數據輸出與圖形繪制功能極大程度上限制了該程序的使用與推廣[11]。為完善這一不足，研發者們借助 R 軟件研發了 R2jags、rjags 及 runjags 等程序包調用 JAGS 軟件來同時實現數據運算與圖形繪制。bmeta 程序包內部數據運算是基于貝葉斯理論且依賴 JAGS 軟件實現的，因此在安裝與加載 bmeta 程序包的同時還需安裝與加載 JAGS 軟件和 R2jags 程序包。目前 JAGS 軟件的最新版本為 JAGS-4.3.0.exe，可從下面的網址下載安裝：https://sourceforge.net/projects/mcmc-jags/files/。
1.3 程序包的加載
當軟件下載完成后，可用下面的命令下載和加載相應的程序包（在加載 bmeta 程序包時會自動加載 R2jags 程序包，因此無需自己加載 R2jags 程序包）：
install.packages("R2jags")
install.packages("bmeta")
library("bmeta")
2 數據預處理與加載
該程序包數據預處理與其他程序包類似，但是要注意如果協變量是分類資料時，需要設置一個基線變量，將其設置為 0，處理后的結果見表 1，然后運行下述代碼加載表 1 的數據：
表1
納入的 13 個預防肺結核病的臨床試驗數據
data = read.csv("bin.csv")
與其他 Meta 分析程序包要求的數據格式不同，bmeta 程序包的數據格式為列表格式，通過下面代碼可以完成從數據框格式到列表格式的轉換：
data.list <- list(y0=data$y0,y1=data$y1,n0=data$n0,n1=data$n1)
3 數據分析
當數據加載完畢時便可進行數據分析，bmeta 程序包進行數據分析的主要函數為“bmeta（）”，該函數可以完成所有類型數據的 Meta 分析和 Meta 回歸，其具體參數如下：
bmeta(data, outcome = c("bin", "ctns", "count"), model = c("std.norm", "std.dt", "reg.norm", "reg.dt", "std.ta", "std.mv", "reg.ta", "reg.mv", "std", "std.unif", "std.hc", "reg", "reg.unif", "reg.hc"), type = c("fix", "ran"), n.iter = 10000, n.burnin = 5000, n.samples = 1000, n.chains = 2, model.file = "model.txt")
在上述代碼中 data 為指定數據集；outcome 為資料類型，包括二分類資料（bin）、連續型資料（ctns）、計數型資料（count）；model 為模型，括號內是適配不同資料的 14 個模型，其中二分類資料包括正態分布和 t 分布的 Meta 分析（std.norm,、std.dt）及 Meta 回歸（reg.norm,、reg.dt）；type 為模型的類型，分為固定效應模型（fix）和隨機效應模型（ran）；n.iter 為迭代的次數，默認值為 10 000；n.burnin 為退火次數，默認為 5 000；n.chains 為運算鏈數，默認值為 2；model.file 為該程序包運行時所產生的 JAGS 代碼文件。
3.1 Meta 分析
應用 bmeta 程序包完成隨機效應模型和固定效應模型比較簡單，當研究沒有合適的先驗分布時，程序包會自動選取無信息的先驗分布進行分析，研究者可在當前工作的文件夾查看其產生 JAGS 代碼文件，本次研究使用 50 000 次的迭代次數，其他參數為默認值，運行的代碼如下：
ran<- bmeta(data=data.list,outcome="bin",model="std.norm",type="ran"，n.iter = 50000)
fix<- bmeta(data=data.list,outcome="bin",model="std.norm",type="fix"，n.iter = 50000)
該程序包運行時的迭代進程見圖 1，運行完畢的結果見圖 2（由于隨機效應模型的結果過長，本文只展示其核心研究結果）、圖 3，其中 alpha[]為固定效應模型中每個研究的 OR 值，gamma[]為隨機效應模型中每個研究的 OR 值；rho 為合并后的總體 OR 值。偏差信息量準則（deviance information criterion，DIC）是等級模型化的赤池信息量準則（Akaike information criterion，AIC），被廣泛應用于由 MCMC 模擬出的后驗分布的貝葉斯模型選擇問題，和 AIC 一樣，DIC 是隨樣本容量增加的漸近近似，只應用于后驗分布呈多元正態分布的情況，一般情況下該值越小證明模型擬合越好。
圖1
bmeta 程序包運行貝葉斯 Meta 分析時的迭代進程圖
圖2
隨機效應模型結果圖
圖3
固定效應模型結果圖
本次 JAGS 運行的隨機效應模型代碼如下：
model {
for (s in 1:S){
y0[s]~dbin(pi0[s],n0[s])
y1[s]~dbin(pi1[s],n1[s])
logit(pi0[s])<-alpha[s]+Z0[s,]%*%beta0
logit(pi1[s])<-alpha[s]+Z0[s,]%*%beta0+delta[s]
alpha[s]~dnorm(0,0.0001)
delta[s]~dnorm(mu,tau)
gamma[s]<-exp(delta[s])
}
###先驗分布###
beta0[1:J]~dmnorm(m.beta0[],tau.beta0[,])
mu~dnorm(0,0.0001)
sigma~dunif(0,5)
tau<-pow(sigma,-2)
rho<-exp(mu)
}
由圖 2 和圖 3 可知本次研究的 Meta 分析結果，隨機效應模型：OR=0.474，95%CI（0.299，0.710）；固定效應模型：OR=0.633，95%CI（0.560，0.839）。固定效應模型和隨機效應模型結果相差比較大，其中固定效應模型的 DIC 值 35 379 908.4，模型擬合效果非常差，主要是由于該數據并不符合固定效應模型假設前提（即所納入的研究均來自同一總體，研究間差異反應的是隨機誤差結果）。
3.2 Meta 回歸
當研究有協變量時可在原有的數據里通過添加 X=cbind(data$) 來進行 Meta 回歸，本研究數據中包括了緯度差異及患者分配方法，將兩個協變量及其組合分別納入進行 Meta 回歸，代碼如下：
d1 <- list(y0=data$y0,y1=data$y1,n0=data$n0,n1=data$n1,X=cbind(data$latitude))
d2 <- list(y0=data$y0,y1=data$y1,n0=data$n0,n1=data$n1,X=cbind(data$allocation))
d3 <- list(y0=data$y0,y1=data$y1,n0=data$n0,n1=data$n1,X=cbind(data$latitude, data$
allocation))
m1 <- bmeta(d1, outcome="bin", model="reg.norm", type="ran",n.iter = 50000)
m2 <- bmeta(d2, outcome="bin", model="reg.norm", type="ran",n.iter = 50000)
m3 <- bmeta(d3, outcome="bin", model="reg.norm", type="ran",n.iter = 50000)
與 Metafor 基于頻率學進行 Meta 回歸不一樣，貝葉斯 Meta 回歸中沒有 I2 和 Q 等統計量，也無法對每個變量做檢驗判斷其是否統計學意義。Bmeta 程序包主要通過 DIC 值來判斷協變量的加入對模型的影響，所以其結果跟上述隨機效應模型的 Meta 分析結果類似，僅在參數上有細微差異，三個模型的核心結果見表 2，由表可知 m2 的 DIC 值和不加入協變量的隨機效應模型相同（DIC=190.5），m1 和 m3 的模型擬合相似且 DIC 值降低，說明分配方式斜協變量的加入對原始模型沒影響，緯度差異的協變量加入模型擬合情況更佳，可以選擇 m1 模型為最終模型。
表2
貝葉斯 Meta 回歸結果
4 圖形繪制
bmeta 包自身不僅可繪制展示結果的森林圖、后驗分布圖等，還可繪制用于檢查模型擬合狀況的軌跡圖和診斷圖等，接下來用上面的隨機效應模型“ran”模型來展示這些圖形的繪制。
4.1 后驗分布圖
由于 bmeta 程序包使用完整的貝葉斯方法進行效應合并，所以可通過繪制后驗分布圖（圖 4）判斷 OR 值的后驗分布情況，代碼如下：
圖4
OR 值的后驗分布圖
posterior.plot(ran,xlim=c(-0.5,2.5),xlab="odds ratio",main="Posterior distribution of pooled odds ratio", scale="exp")
除了通過 OR 值得到后驗分布圖外，bmeta 程序包還可使用研究之間標準差（standard deviation，SD）后驗分布圖來檢驗研究之間的異質性（圖 5）。SD 值越小說明研究之間的異質性越小，當研究之間的 SD 接近 0 時可以使用固定效應模型，實現代碼如下：
圖5
研究之間標準差的后驗分布圖
posterior.plot(ran,heterogeneity=TRUE,xlim=c(0,3),xlab="between-study standard deviation")
4.2 森林圖
在 bmeta 包里森林圖（圖 6）的繪制代碼如下：
圖6
森林圖
forest.plot(ran,boxsize=0.2,study.label=c(paste0(data$Study,",",data$Year),"Summary estimate"), title="forestplot for BCG Vaccine")
在上面代碼中加入 add.null=TURE 后可以繪制隨機效應模型和沒有合并效應量（no-pooling effects）的模型對比森林圖（圖 7），完整代碼如下：
圖7
對比森林圖
forest.plot(ran,add.null=TRUE,boxsize=0.2,study.label=c(paste0(data$Study,",",data$Year),"Summary estimate"), title="Two-line forestplot for comparison")
4.3 漏斗圖
漏斗圖（圖 8）的繪制代碼為：funnel.plot(ran,title="funnel plot")。
圖8
漏斗圖
4.4 軌跡圖
軌跡圖主要用來評估模型參數的 MCMC 收斂情況，其橫軸顯示迭代的次數，縱軸顯示參數的后驗分布值，當 MCMC 達到穩態時，參數的模擬值會在均值附近以相同的幅度上下波動（圖 9）。在 bmeta 包里繪制軌跡圖的代碼為：traceplot.bmeta(ran,"mu",lab="mu")，其中“mu”為想要評估的模型參數，lab 后面為縱軸的名稱。
圖9
軌跡圖
4.5 診斷圖
診斷圖主要用來評估模型參數的收斂情況，bmeta 程序包采用 Gelman-Rubin 準則來判斷參數的收斂情況。在收斂的情況下，Ghat 值應該隨著交互次數的增加逐漸減小到零，至少應低于 1.05，如果收斂情況不好，則要增加迭代次數。該圖的實現代碼為 diag.plot(ran)，通過圖 10 可以看出本次模型參數收斂情況較好。
圖10
收斂診斷圖
4.6 自相關圖
自相關性圖橫軸表示滯后時間，縱軸仍為后驗分布值。該圖顯示了 MCMC 樣本自相關性從滯后 0 開始的一系列滯后范圍，在滯后 0 這一點上的值等于 MCMC 的樣本方差（圖 11）。在 bmeta 包里繪制該圖的代碼為：acf.plot(ran,"gamma[1]")，其中 gamma[1]為想要評估的參數。
圖11
自相關圖
5 討論
貝葉斯分析是用于統計建模、結果解釋和數據預測的強大分析工具。相比于經典的頻率學派，貝葉斯 Meta 分析有以下優點：① 可充分結合先驗信息、樣本信息及總體信息，獲得后驗分布[12]；② 當納入研究數量較少時，貝葉斯方法在計算研究之間的異質性和合并效應值方面更有優勢[13]；③ 頻率學派大多數情況下是基于大樣本漸近分布做出統計推斷，僅能計算效應量的置信區間，而貝葉斯學派則可直接計算精確的有限樣本分布，且充分考慮了模型的不確定性[14]。
目前，R 軟件里用于頭對頭 Meta 分析的 metafor 和 meta 程序包都是基于頻率學派[8, 9]，用于貝葉斯 Meta 分析的程序包有 bmeta 和 bayesmeta 程序包，國內學者張天嵩介紹了如何使用 bayesmeta 程序包實現單臂試驗連續型數據的 Meta 分析[15]。本研究通過一個經典的例子展示了如何使用 bmeta 程序包實現二分類資料的貝葉斯 Meta 分析與 Meta 回歸。通過本文演示可看出當前的 bmeta 程序包相對 BUGS 軟件具有以下優點：① 操作簡單，相比 bayesmeta 程序包，該程序包僅用一個函數便可實現所有類型數據的模型分析，并省去了用戶在 BUGS 軟件建模的過程，大大簡化了操作流程；② R 軟件里強大的繪圖功能在 bmeta 程序包中體現，用戶可根據自己的研究需求和喜好來定制個性化圖形。另外除了二分類資料外，bmeta 程序包也可實現連續性資料和計數資料及其相應模型的 Meta 分析，由于篇幅所限本文并未展示，讀者可參考其幫助文檔。另外該程序包還存在一些缺點，首先，該程序包采用貝葉斯方法，每次運行隨機抽樣使結果都有細微變化，讀者可使用 seed（）函數設立種子數或增大迭代次數來使結果穩定；其次，該程序包繪制的圖形雖然豐富，但其功能仍需進一步完善，才能滿足使用者的需求。
綜上，R 軟件 bmeta 程序包充分結合了 JAGS 軟件的計算功能與 R 軟件特有的數據整合與強大的圖形繪制等功能，從而實現了結果便利輸出與配套優質結果圖的繪制。對那些希望實現貝葉斯 Meta 分析但不擅長使用 BUGS 的研究者，bmeta 程序包無疑是更好的選擇。
Meta 分析作為一種整合單個研究效應量進行證據合并的常用統計方法，在循證醫學中占有重要地位[1]。貝葉斯 Meta 分析是基于貝葉斯統計發展起來的一種的 Meta 分析方法，主要采用“馬爾科夫鏈—蒙特卡羅”（Markov chain Monte Carlo，MCMC）方法，因其在處理復雜隨機效應、分層結構或是稀疏數據時比頻率學 Meta 分析方法更有優勢，目前越來越受歡迎。董圣杰等[2]介紹了貝葉斯 Meta 分析的建模過程及如何在 WinBUGS 軟件實現貝葉斯 Meta 分析。但是 BUGS 編程復雜，操作過程繁瑣并且無法繪制森林圖、漏斗圖等 Meta 分析必不可少的圖形。R 軟件作為一款開源且繪圖功能強大的軟件能夠實現幾乎所有的 Meta 分析類型[3, 4]。然而目前 R 軟件基于貝葉斯方法計算的程序包更多是用來進行網狀 Meta 分析[5-7]，用來進行傳統頭對頭 Meta 分析的程序包的理論原理均基于經典的頻率學派[8, 9]。在這樣的背景需求下，bmeta 程序包被研發出來，該程序包采用簡潔的代碼輸入的方式，并通過調用 JAGS 來實行貝葉斯 Meta 分析和 Meta 回歸。本文以《Efficacy of BCG vaccine in the prevention of tuberculosis meta-analysis of the published literature》[10]一文的數據為例來展示該程序包的具體操作。
1 軟件安裝及程序包加載
1.1 R 軟件安裝
截至完稿時，R 軟件最新版本為 R-4.0.0，用戶可在下面的網址根據自己計算機的操作系統下載最新版本：https://www.r-project.org/。
1.2 JAGS 軟件的下載
JAGS（Just another Gibbs sampler）軟件是一款基于貝葉斯理論運算的編程軟件，但尚不完善的數據輸出與圖形繪制功能極大程度上限制了該程序的使用與推廣[11]。為完善這一不足，研發者們借助 R 軟件研發了 R2jags、rjags 及 runjags 等程序包調用 JAGS 軟件來同時實現數據運算與圖形繪制。bmeta 程序包內部數據運算是基于貝葉斯理論且依賴 JAGS 軟件實現的，因此在安裝與加載 bmeta 程序包的同時還需安裝與加載 JAGS 軟件和 R2jags 程序包。目前 JAGS 軟件的最新版本為 JAGS-4.3.0.exe，可從下面的網址下載安裝：https://sourceforge.net/projects/mcmc-jags/files/。
1.3 程序包的加載
當軟件下載完成后，可用下面的命令下載和加載相應的程序包（在加載 bmeta 程序包時會自動加載 R2jags 程序包，因此無需自己加載 R2jags 程序包）：
install.packages("R2jags")
install.packages("bmeta")
library("bmeta")
2 數據預處理與加載
該程序包數據預處理與其他程序包類似，但是要注意如果協變量是分類資料時，需要設置一個基線變量，將其設置為 0，處理后的結果見表 1，然后運行下述代碼加載表 1 的數據：
表1
納入的 13 個預防肺結核病的臨床試驗數據
data = read.csv("bin.csv")
與其他 Meta 分析程序包要求的數據格式不同，bmeta 程序包的數據格式為列表格式，通過下面代碼可以完成從數據框格式到列表格式的轉換：
data.list <- list(y0=data$y0,y1=data$y1,n0=data$n0,n1=data$n1)
3 數據分析
當數據加載完畢時便可進行數據分析，bmeta 程序包進行數據分析的主要函數為“bmeta（）”，該函數可以完成所有類型數據的 Meta 分析和 Meta 回歸，其具體參數如下：
bmeta(data, outcome = c("bin", "ctns", "count"), model = c("std.norm", "std.dt", "reg.norm", "reg.dt", "std.ta", "std.mv", "reg.ta", "reg.mv", "std", "std.unif", "std.hc", "reg", "reg.unif", "reg.hc"), type = c("fix", "ran"), n.iter = 10000, n.burnin = 5000, n.samples = 1000, n.chains = 2, model.file = "model.txt")
在上述代碼中 data 為指定數據集；outcome 為資料類型，包括二分類資料（bin）、連續型資料（ctns）、計數型資料（count）；model 為模型，括號內是適配不同資料的 14 個模型，其中二分類資料包括正態分布和 t 分布的 Meta 分析（std.norm,、std.dt）及 Meta 回歸（reg.norm,、reg.dt）；type 為模型的類型，分為固定效應模型（fix）和隨機效應模型（ran）；n.iter 為迭代的次數，默認值為 10 000；n.burnin 為退火次數，默認為 5 000；n.chains 為運算鏈數，默認值為 2；model.file 為該程序包運行時所產生的 JAGS 代碼文件。
3.1 Meta 分析
應用 bmeta 程序包完成隨機效應模型和固定效應模型比較簡單，當研究沒有合適的先驗分布時，程序包會自動選取無信息的先驗分布進行分析，研究者可在當前工作的文件夾查看其產生 JAGS 代碼文件，本次研究使用 50 000 次的迭代次數，其他參數為默認值，運行的代碼如下：
ran<- bmeta(data=data.list,outcome="bin",model="std.norm",type="ran"，n.iter = 50000)
fix<- bmeta(data=data.list,outcome="bin",model="std.norm",type="fix"，n.iter = 50000)
該程序包運行時的迭代進程見圖 1，運行完畢的結果見圖 2（由于隨機效應模型的結果過長，本文只展示其核心研究結果）、圖 3，其中 alpha[]為固定效應模型中每個研究的 OR 值，gamma[]為隨機效應模型中每個研究的 OR 值；rho 為合并后的總體 OR 值。偏差信息量準則（deviance information criterion，DIC）是等級模型化的赤池信息量準則（Akaike information criterion，AIC），被廣泛應用于由 MCMC 模擬出的后驗分布的貝葉斯模型選擇問題，和 AIC 一樣，DIC 是隨樣本容量增加的漸近近似，只應用于后驗分布呈多元正態分布的情況，一般情況下該值越小證明模型擬合越好。
圖1
bmeta 程序包運行貝葉斯 Meta 分析時的迭代進程圖
圖2
隨機效應模型結果圖
圖3
固定效應模型結果圖
本次 JAGS 運行的隨機效應模型代碼如下：
model {
for (s in 1:S){
y0[s]~dbin(pi0[s],n0[s])
y1[s]~dbin(pi1[s],n1[s])
logit(pi0[s])<-alpha[s]+Z0[s,]%*%beta0
logit(pi1[s])<-alpha[s]+Z0[s,]%*%beta0+delta[s]
alpha[s]~dnorm(0,0.0001)
delta[s]~dnorm(mu,tau)
gamma[s]<-exp(delta[s])
}
###先驗分布###
beta0[1:J]~dmnorm(m.beta0[],tau.beta0[,])
mu~dnorm(0,0.0001)
sigma~dunif(0,5)
tau<-pow(sigma,-2)
rho<-exp(mu)
}
由圖 2 和圖 3 可知本次研究的 Meta 分析結果，隨機效應模型：OR=0.474，95%CI（0.299，0.710）；固定效應模型：OR=0.633，95%CI（0.560，0.839）。固定效應模型和隨機效應模型結果相差比較大，其中固定效應模型的 DIC 值 35 379 908.4，模型擬合效果非常差，主要是由于該數據并不符合固定效應模型假設前提（即所納入的研究均來自同一總體，研究間差異反應的是隨機誤差結果）。
3.2 Meta 回歸
當研究有協變量時可在原有的數據里通過添加 X=cbind(data$) 來進行 Meta 回歸，本研究數據中包括了緯度差異及患者分配方法，將兩個協變量及其組合分別納入進行 Meta 回歸，代碼如下：
d1 <- list(y0=data$y0,y1=data$y1,n0=data$n0,n1=data$n1,X=cbind(data$latitude))
d2 <- list(y0=data$y0,y1=data$y1,n0=data$n0,n1=data$n1,X=cbind(data$allocation))
d3 <- list(y0=data$y0,y1=data$y1,n0=data$n0,n1=data$n1,X=cbind(data$latitude, data$
allocation))
m1 <- bmeta(d1, outcome="bin", model="reg.norm", type="ran",n.iter = 50000)
m2 <- bmeta(d2, outcome="bin", model="reg.norm", type="ran",n.iter = 50000)
m3 <- bmeta(d3, outcome="bin", model="reg.norm", type="ran",n.iter = 50000)
與 Metafor 基于頻率學進行 Meta 回歸不一樣，貝葉斯 Meta 回歸中沒有 I2 和 Q 等統計量，也無法對每個變量做檢驗判斷其是否統計學意義。Bmeta 程序包主要通過 DIC 值來判斷協變量的加入對模型的影響，所以其結果跟上述隨機效應模型的 Meta 分析結果類似，僅在參數上有細微差異，三個模型的核心結果見表 2，由表可知 m2 的 DIC 值和不加入協變量的隨機效應模型相同（DIC=190.5），m1 和 m3 的模型擬合相似且 DIC 值降低，說明分配方式斜協變量的加入對原始模型沒影響，緯度差異的協變量加入模型擬合情況更佳，可以選擇 m1 模型為最終模型。
表2
貝葉斯 Meta 回歸結果
4 圖形繪制
bmeta 包自身不僅可繪制展示結果的森林圖、后驗分布圖等，還可繪制用于檢查模型擬合狀況的軌跡圖和診斷圖等，接下來用上面的隨機效應模型“ran”模型來展示這些圖形的繪制。
4.1 后驗分布圖
由于 bmeta 程序包使用完整的貝葉斯方法進行效應合并，所以可通過繪制后驗分布圖（圖 4）判斷 OR 值的后驗分布情況，代碼如下：
圖4
OR 值的后驗分布圖
posterior.plot(ran,xlim=c(-0.5,2.5),xlab="odds ratio",main="Posterior distribution of pooled odds ratio", scale="exp")
除了通過 OR 值得到后驗分布圖外，bmeta 程序包還可使用研究之間標準差（standard deviation，SD）后驗分布圖來檢驗研究之間的異質性（圖 5）。SD 值越小說明研究之間的異質性越小，當研究之間的 SD 接近 0 時可以使用固定效應模型，實現代碼如下：
圖5
研究之間標準差的后驗分布圖
posterior.plot(ran,heterogeneity=TRUE,xlim=c(0,3),xlab="between-study standard deviation")
4.2 森林圖
在 bmeta 包里森林圖（圖 6）的繪制代碼如下：
圖6
森林圖
forest.plot(ran,boxsize=0.2,study.label=c(paste0(data$Study,",",data$Year),"Summary estimate"), title="forestplot for BCG Vaccine")
在上面代碼中加入 add.null=TURE 后可以繪制隨機效應模型和沒有合并效應量（no-pooling effects）的模型對比森林圖（圖 7），完整代碼如下：
圖7
對比森林圖
forest.plot(ran,add.null=TRUE,boxsize=0.2,study.label=c(paste0(data$Study,",",data$Year),"Summary estimate"), title="Two-line forestplot for comparison")
4.3 漏斗圖
漏斗圖（圖 8）的繪制代碼為：funnel.plot(ran,title="funnel plot")。
圖8
漏斗圖
4.4 軌跡圖
軌跡圖主要用來評估模型參數的 MCMC 收斂情況，其橫軸顯示迭代的次數，縱軸顯示參數的后驗分布值，當 MCMC 達到穩態時，參數的模擬值會在均值附近以相同的幅度上下波動（圖 9）。在 bmeta 包里繪制軌跡圖的代碼為：traceplot.bmeta(ran,"mu",lab="mu")，其中“mu”為想要評估的模型參數，lab 后面為縱軸的名稱。
圖9
軌跡圖
4.5 診斷圖
診斷圖主要用來評估模型參數的收斂情況，bmeta 程序包采用 Gelman-Rubin 準則來判斷參數的收斂情況。在收斂的情況下，Ghat 值應該隨著交互次數的增加逐漸減小到零，至少應低于 1.05，如果收斂情況不好，則要增加迭代次數。該圖的實現代碼為 diag.plot(ran)，通過圖 10 可以看出本次模型參數收斂情況較好。
圖10
收斂診斷圖
4.6 自相關圖
自相關性圖橫軸表示滯后時間，縱軸仍為后驗分布值。該圖顯示了 MCMC 樣本自相關性從滯后 0 開始的一系列滯后范圍，在滯后 0 這一點上的值等于 MCMC 的樣本方差（圖 11）。在 bmeta 包里繪制該圖的代碼為：acf.plot(ran,"gamma[1]")，其中 gamma[1]為想要評估的參數。
圖11
自相關圖
5 討論
貝葉斯分析是用于統計建模、結果解釋和數據預測的強大分析工具。相比于經典的頻率學派，貝葉斯 Meta 分析有以下優點：① 可充分結合先驗信息、樣本信息及總體信息，獲得后驗分布[12]；② 當納入研究數量較少時，貝葉斯方法在計算研究之間的異質性和合并效應值方面更有優勢[13]；③ 頻率學派大多數情況下是基于大樣本漸近分布做出統計推斷，僅能計算效應量的置信區間，而貝葉斯學派則可直接計算精確的有限樣本分布，且充分考慮了模型的不確定性[14]。
目前，R 軟件里用于頭對頭 Meta 分析的 metafor 和 meta 程序包都是基于頻率學派[8, 9]，用于貝葉斯 Meta 分析的程序包有 bmeta 和 bayesmeta 程序包，國內學者張天嵩介紹了如何使用 bayesmeta 程序包實現單臂試驗連續型數據的 Meta 分析[15]。本研究通過一個經典的例子展示了如何使用 bmeta 程序包實現二分類資料的貝葉斯 Meta 分析與 Meta 回歸。通過本文演示可看出當前的 bmeta 程序包相對 BUGS 軟件具有以下優點：① 操作簡單，相比 bayesmeta 程序包，該程序包僅用一個函數便可實現所有類型數據的模型分析，并省去了用戶在 BUGS 軟件建模的過程，大大簡化了操作流程；② R 軟件里強大的繪圖功能在 bmeta 程序包中體現，用戶可根據自己的研究需求和喜好來定制個性化圖形。另外除了二分類資料外，bmeta 程序包也可實現連續性資料和計數資料及其相應模型的 Meta 分析，由于篇幅所限本文并未展示，讀者可參考其幫助文檔。另外該程序包還存在一些缺點，首先，該程序包采用貝葉斯方法，每次運行隨機抽樣使結果都有細微變化，讀者可使用 seed（）函數設立種子數或增大迭代次數來使結果穩定；其次，該程序包繪制的圖形雖然豐富，但其功能仍需進一步完善，才能滿足使用者的需求。
綜上，R 軟件 bmeta 程序包充分結合了 JAGS 軟件的計算功能與 R 軟件特有的數據整合與強大的圖形繪制等功能，從而實現了結果便利輸出與配套優質結果圖的繪制。對那些希望實現貝葉斯 Meta 分析但不擅長使用 BUGS 的研究者，bmeta 程序包無疑是更好的選擇。
