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 程序包无疑是更好的选择。
