传统 Meta 分析的统计基础是频率学理论，这种 Meta 分析通常采用 Q 统计量衡量不同研究之间的异质性大小，然后根据异质性检验结果来确定采用何种模型（固定效应模型或随机效应模型）进行 Meta 分析。经典频率学随机效应模型可通过 Q 检验获得研究间方差的矩估计，但是无法获得研究间协方差的 95% 可信区间（confidence interval，CI），而且基于频率学的 Meta 分析检验效能较低。在传统的 Meta 分析模型中，效应量指标的分布一般假设服从正态分布。例如，对于二分类资料的 Meta 分析而言，最常用的效应量为比值比（odds ratio，OR）^[1]。经典 Meta 分析方法是将 OR 对数化，然后将 log（OR）设为效应值，并假设其服从正态分布^[2]。这种 Meta 分析模型是以正态分布假设为前提的，因此在存在小样本资料时，如果不符合近似正态的条件，经典方法合并效应值存在一定误差^[3]。此外，经典 Meta 分析模型很难胜任复杂的统计模型，并且无法识别和处理极端值^[4]。

贝叶斯 Meta 分析（Bayesian meta-analysis，BMA）是近年来以贝叶斯统计理论为基础发展起来的一种新型的 Meta 分析方法，主要采用“马尔科夫链-蒙特卡罗”（Markov chain Monte Carlo，MCMC）方法进行 Meta 分析^[5]。贝叶斯方法与频率论（经典）方法的原理完全不同，频率论认为概率反映的是某事件在重复试验中出现的次数与试验次数之比，也就是事件发生的频率，具有随机性和稳定性特点，这是频率论的统计基础；而贝叶斯方法认为概率是数值的可信程度。这意味着，频率论将模型参数设置为固定，而将数据设置为随机；贝叶斯则与之相反^[6]。贝叶斯统计学与传统频率学在统计推断的理念和方法上存在本质上的区别，包括对信息的利用、对未知参数 θ 的解释及统计推断理念上的差异等^{[4, 7]}。

Meta 分析同样可采用贝叶斯分析实现。贝叶斯统计将参数 θ 作为一个随机变量，有一定的先验分布，在获得样本之后（给定的样本信息），θ 的后验分布 π（θ|x）应包含 θ 的综合信息，可从 θ 的后验分布获得参数 θ 的贝叶斯估计。贝叶斯分析能很好地解决经典 Meta 分析存在的缺陷，近年来，采用贝叶斯方法的 Meta 分析引起了广泛重视，基于贝叶斯的网络 Meta 分析是 Meta 分析领域的研究热点。贝叶斯 Meta 分析最常用的统计软件是 WinBUGS、R 软件，但实际上 SAS 提供了两种过程来实现贝叶斯 Meta 分析，分别是 PROC GENMOD 与 PROC MCMC^[8]。本研究团队先前已介绍了利用 PROC MCMC 实现诊断性试验的贝叶斯 Meta 分析方法^[6]。然而，目前尚无学者介绍利用 SAS 软件实现二分类数据贝叶斯 Meta 分析的方法。本文通过实例数据介绍利用 SAS 软件中的 PROC MCMC 实现二分类数据贝叶斯 Meta 分析，同时将结果与基于频率学 Meta 分析方法的结果进行对比，相关的 SAS 代码在文后附录中给出。

1   贝叶斯 Meta 分析理论简介

1.1   贝叶斯统计理论简介

贝叶斯分析是基于贝叶斯定理发展起来用于系统阐述和解决统计问题的方法^[3]。其理论简述如下：根据既往经验，假设某事件 θ 发生的概率为 P（θ），我们这个概率称为先验概率。然后我们现有一组新的数据 y，则 y 在 y/θ 的前提下发生的条件概率记为 P（y/θ），此条件概率称为似然。根据先验概率和似然可计算出概率 P（θ/y），表示 θ 在 y 存在的前提下发生事件的可能性大小，称为后验概率。后验概率与先验概率与似然的乘积成正比，即 P（θ/y）=α×P（y/θ）×P（θ）。在贝叶斯框架下，参数 θ 是一个随机变量，服从一定的先验分布，在给定的样本信息（抽样）条件下，通过抽样获取样本，则 θ 的后验分布 应包含 θ 的综合信息，因此可通过 θ 的后验分布来获得参数 θ 的贝叶斯估计。因此，贝叶斯统计是由模型、参数和似然 3 个要素构成的统计方法。

1.2   二分类数据贝叶斯 Meta 分析的建模方法

二分类数据 Meta 分析最常用的效应指标是 OR，经过对数转换后的 OR 值服从正态分布。与频数法相对应，本文贝叶斯 Meta 分析的建模方法同样基于 OR 的对数转换值。具体模型构建如下：

假设某 Meta 分析共纳入 n 个研究，nt、rt、nc、rc 分别为治疗组和对照组总人数及事件发生例数。 为第 i 个研究的效应量，本文为 OR 的对数转换值， 为第 i 个研究的效应量的研究内方差，theta 为研究的真实效应值。上述参数建模如下：

固定效应模型假设所有纳入 Meta 分析的研究来自同一个总体，因此效应值固定。对于固定效应模型，，其中 是固定效应模型背景下的真实效应值，是待估计的固定效应模型参数。因此，贝叶斯固定效应模型可建模如下：

，需要获得的数据集包括：；。模型待估计的参数为 。似然函数为正态分布密度 。后验分布为 。

但由于纳入 Meta 分析的原始研究通常存在异质性，导致使用固定效应模型的假设条件往往不成立。这时采用随机效应模型，即假设纳入 Meta 分析的研究可能来自不同的研究总体，每一个研究都拥有不同的效应值，这些效应值服从同一分布的总体。对于随机效应模型，，而，方差 反映了研究间变异的度量，代表随机效应模型的随机部分。贝叶斯随机效应模型可建模如下：

，。需要获得的数据集包括：；。模型待估计的参数包括 ， 和 。其中 表示每个研究的效应， 为真实效应值； 为研究的方差。似然函数与固定效应模型相似。后验分布为 。

2   资料与方法

2.1   示例数据介绍

该示例数据来源于陈艳丽等^[9]发表的 Meta 分析，该 Meta 分析旨在评价腹腔镜下根治性全子宫切除对比开腹手术对早期宫颈癌的疗效和预后的影响。本研究提取该文中并发症的发生率数据（表 1），采用贝叶斯 Meta 分析比较两种手术方式对术后并发症发生率的影响。

将表 1 的数据录入 SAS 软件，并存储在名为 Dichotomous_data 的数据集中，数据集录入的代码见附录。为拟合本文的 BMA 模型，需要的数据集变量包括 author、year、rt、Nt、rc、Nc；分别代表作者、发表年份、研究组事件数（在本示例中，事件代表术后发生并发症）、研究组样本量、对照组事件数和对照组样本量。

2.2   SAS 代码介绍

利用附录代码 1 完成原始数据集的录入之后，再利用附录代码 2 计算效应值 di 及标准误 sigmai2。需要注意的是，本文示例数据不存在零事件试验，对于零事件试验，需要进行连续性校正，校正方法是 2×2 四格表每个单元格增加 0.5；然后使用 PROC MCMC 程序，分别拟合固定效应模型和随机效应模型的贝叶斯 Meta 分析。

在固定效应模型 BMA 中，模型将迭代 50 000 次，前 5 000 次退火用于消除初始值的影响。固定效应值 θ 的先验分布被设置为均数为 0，方差为 1×10⁶的正态分布。不同研究观测的效应值 di 被建模为模型均数为 θ，方差为 1×10⁶的正态分布。与固定效应模型不同的是，随机效应模型引入一个反映研究间变异的度量 tau2，并且 tau2 的先验分布被假设服从形状参数（shape parameter）为 0.001、逆尺度参数为 0.001 的伽玛分布（Gamma distribution）。在随机效应模型背景下，通过“random deltai ～ normal（mean=theta，var=tau2）subject=study monitor=（deltai）；”语句设置随机效应。

为保证基于 SAS 软件的贝叶斯 Meta 分析具有实用性，本文给出了 SAS 软件背景下的森林图绘制代码^[10]。森林图绘制代码将笔者先前发表的 SAS 代码包装成两个宏程序，分别为&forestplot 和&forestdata。前者用于绘制森林图，其宏参数是需要绘制森林图的 SAS 数据集 Forest；后者用于构建森林图数据集，并调用宏程序 forestplot 实现森林图绘制；其宏参数是 PROC MCMC 过程步完成数据分析后输出的后验摘要数据集 PostSummaries_data 和指定的 Meta 分析模型 type（type=1 表示固定效应模型，type=2 表示随机效应模型）。

3   结果

3.1   固定效应模型贝叶斯 Meta 分析结果

固定效应模型贝叶斯 Meta 分析结果见表 2，表中的 theta 和 OR 是对数转换关系，即 OR=exp（theta）。表 3 是 PROC MCMC 提供的参数值的 95% 可信区间。图 1 是 MCMC 采样分析图，分别包括了 Gibbs 抽样动态踪迹图、自相关图和后验分布核密度估计图。从图 1 中的 Gibbs 动态图可看出，参数 OR 经过 50 000 次迭代后基本收敛。在自相关图中，OR 的验后自相关非常小，模型可靠性高。核密度图显示，OR 的后验分布类似于正态分布，与研究预期相似。

3.2   随机效应模型贝叶斯 Meta 分析结果

表 4 给出了随机效应模型贝叶斯 Meta 分析结果。除了给出 θ 的估计值外，随机效应模型还给出了各个研究的随机效应值。表 4 中，除了 OR 值外，其他需要经过拟对数转换后应用于分析。随机效应模型的 MCMC 采样分析图见图 2，除了图 2，随机效应模型同时还会对其他参数进行 MCMC 采样分析，读者可自行运行本文的 SAS 代码查看。

3.3   不同 Meta 分析方法的结果比较

表 5 给出了频率学方法和贝叶斯方法的 Meta 分析结果，两种方法给出的模型参数估计值极其相似。

3.4   贝叶斯 Meta 分析的森林图结果

森林图可非常直观地观察到 Meta 分析的结果。但 PROC MCMC 过程步无法直接给出森林图，为便于读者利用 SAS 软件开展研究，本文开发了 SAS 软件绘制森林图的代码（见附录），运行代码后，可获得固定效应模型和随机效应模型的森林图。图 3 显示了固定效应模型的森林图。

4   讨论

贝叶斯统计和经典统计学派（频率学派）是目前两大主要的统计学派。贝叶斯统计比频率学统计方法需要更强大的计算能力来完成推断，因此贝叶斯统计是以计算机运算能力发展为前提的。随着计算机技术的发展，贝叶斯统计学在医学领域内的应用越来越广泛，其中也包括 Meta 分析。贝叶斯统计是综合未知参数的先验信息和样本信息，根据贝叶斯定理，求出后验分布，然后根据后验分布对未知参数进行统计推断，最后获得参数估计值的统计方法。在方法学上，频率统计学方法以样本为基础，通过假设检验和统计计算的方法进行统计推断，最后获得结论；而贝叶斯方法则是基于贝叶斯原理，根据先验概率推断后验概率^[11]。虽然传统 Meta 分析通常采用频率学方法实现，但越来越多的学者在探索利用贝叶斯方法实现 Meta 分析的方法。目前已有不少统计软件可实现贝叶斯 Meta 分析，包括 R、Stata、WinBugs 等众多软件^{[4, 8, 12-16]}。Winbugs 是用于 Gibbs 抽样的专用软件包，是免费软件，目前已广泛应用于实施贝叶斯方法^[4]。既往学者研究表明，与频率学 Meta 分析相比，贝叶斯 Meta 分析具有以下优点：① 建模方式更灵活，在多种贝叶斯分析软件中可操作性强；② 贝叶斯 Meta 分析充分考虑了模型参数（如方差等）的不确定性，例如可估计研究间方差的不精确性，最终获得的合并效应量点估计及 95% 可信区间来源于参数的后验分布，结果可靠；③ 贝叶斯 Meta 分析提供的统计检验，反映在某种程度上可改变人们的观念^[17]。此外，贝叶斯分析方法最大的优势是不受样本量限制，因为贝叶斯的基本原理就是通过先验概率推断后验概率；而频率学派反映的是抽样概率，因此需要样本量的支持。而实践中往往纳入 Meta 分析的原始研究数量极为有限，这为贝叶斯统计提供了更广阔的应用范围。贝叶斯 Meta 分析不足之处有以下两点：① 频率学派认为一个事件的概率可通过大量重复试验下事件出现的频率来解释，这种概率不依赖于主观性，被认为是客观概率，而贝叶斯统计需要指定先验分布，而这种分布带有主观性，因此被称为主观概率；因此，其参数的先验分布受统计分析师建模思想的影响；② 尽管贝叶斯统计并不接受频率观点，但贝叶斯统计仍然需要依赖样本分布，而样本则来源于抽样。在实际应用中，往往很难完全区分频率学派和贝叶斯统计孰优孰劣。在 Meta 分析领域，初学者往往难以把握统计建模，应以频率学统计结果为准；而具有统计建模基础的系统评价员，则应更多尝试贝叶斯统计方法。

虽然 SAS 软件在统计学领域具有广泛的应用基础，但基于 SAS 软件的贝叶斯 Meta 分析却起步较晚。目前介绍 SAS 软件实现贝叶斯 Meta 分析的文献并不多。在 SAS 9.2 版本之前，SAS 无法实现贝叶斯模型分析，但在 SAS 9.2 之后，SAS 软件提供了多种过程步来实现贝叶斯 Meta 分析，最常用的过程步包括 PROC GENMOD 和 PROC MCMC。PROC GENMOD 可通过“BAYES”语句实现贝叶斯 Meta 分析；但这种贝叶斯分析与 PROC MCMC 的贝叶斯分析基础不同，前者类似于频率学方法，后者是基于 MCMC 的统计方法。因此，在 SAS 中，贝叶斯 Meta 分析首选 PROC MCMC^[18]。国内外已有不少学者介绍了基于 PROC MCMC 实现不同数据的贝叶斯 Meta 分析^{[6, 8, 19]}；但目前尚无介绍 PROC MCMC 实现二分类数据 Meta 分析的文献。鉴于二分类数据在医学领域极为常见，有必要介绍二分类数据的贝叶斯 Meta 分析。曾平等^[20]详细介绍了二分类数据的贝叶斯 Meta 分析的建模方法，本文的侧重介绍该模型在 SAS 软件中的实现过程。

本研究详细介绍了利用 PROC MCMC 实现二分类数据贝叶斯 Meta 分析的 SAS 程序。通过示例数据分析显示，贝叶斯 Meta 分析几乎提供了与频率学 Meta 分析完全一致的结果。但需要注意几点：① PROC MCMC 是以马尔科夫链-蒙特卡罗方法为基础的。MCMC 是在贝叶斯理论框架下，通过计算机进行模拟的蒙特卡洛方法（Monte Carlo）。该方法将马尔科夫（Markov）过程引入到 Monte Carlo 模拟中，实现抽样分布随模拟的进行而改变的动态模拟，弥补了传统的蒙特卡罗积分只能静态模拟的缺陷。② 在多数情况下，Metropolis-Hastings 抽样（M-H sampling）是 MCMC 常用的算法；但对于高维的情形，M-H 算法效率太低。Meta 分析正是一种特殊的高维形式，因此 Gibbs 抽样应用更广泛。③ PROC MCMC 提供了两种不同的后验可信区间（posterior intervals），分别为等尾可信区间（equal-tail intervals，ET）和最高后验密度可信区间（highest posterior density intervals，HPD）。读者在使用可信区间时，应该首选 ET 可信区间，因为 ET 可信区间更窄，而 HPD 可信区间更宽。④ 贝叶斯 Meta 分析不提供统计学概率值（P值），因此不能利用P值来判断效应值是否具有统计学意义，一般来说，贝叶斯通过后验分布获得的可信区间来反映效应值的统计学意义。以 OR 指标为例，根据既往统计学理论，OR 效应值是否具有统计学意义以是否跨越“1”为标准，当可信区间包含了“1”，就认为效应值无统计学意义。因此，结合本文的示例数据，随机效应模型贝叶斯 Meta 分析获得的 OR 值为 0.53，95% 可信区间为（0.30，0.85），该可信区间并未包含 1，因此具有统计学意义。其临床意义为：与开腹手术相比，腹腔镜下宫颈癌根治术的术后并发症发生率明显降低，差异具有统计学意义。最后，与目前多数贝叶斯统计软件相比，SAS 软件编程功能强大。基于 SAS 的编程基础，本研究开发了一段 SAS 宏程序，可非常方便地实现森林图绘制，这有效解决了目前多数贝叶斯统计软件无法获得森林图的问题。

总之，基于 SAS 软件强大的编程能力，PROC MCMC 可非常方便地实现二分类数据的贝叶斯 Meta 分析，随着贝叶斯统计理论的快速发展，贝叶斯 Meta 分析在 Meta 分析领域将发挥越来越重要的作用。

附录代码

/*附录代码 1：数据集建立及资料录入*/

Data Dichotomous_data；

　length author $20.；

　study=_N_；

　input author year rtNtrcNc；

　cards；

Lee2002430930

……

；

Run；

/*Proc print data=Dichotomous_data noobs；run；*/

/*附录代码 2：计算效应值 di 及标准误 sigmai2*/

Data Bayesian_Dichotomous_data；

　set Dichotomous_data；

　Di=log（（rt/（Nt-rt））/（rc/（Nc-rc）））；

　sigmai2=1/rt+1/（Nt-rt）+1/rc+1/（Nc-rc）；

run；

proc print data=Bayesian_Dichotomous_data noobs；run；

/*附录代码 3：固定效应模型贝叶斯 Meta 分析*/

ods output PostSummaries =mcmc_Summaries_fixed；

ods graphics on；

Proc mcmc data=Bayesian_Dichotomous_data nbi=5000 nmc=50000 seed=1234

monitor=（theta OR）stats（percentage=（2.5 97.5））=all dic；

parms theta 0；

prior theta ～ normal（mean = 0，var = 1.0e6）；

model di ～ normal（theta，var = sigmai2）；

OR=exp（theta）；

run；

ods graphics off；

proc print data=mcmc_Summaries_fixed；run；

/*附录代码 4：森林图绘制代码*/

%macro forestplot（forestdata）；

data forest；

　　set & forestdata.；

　　format OddsRatio 6.2；

　　format Lowercl 6.2；format Uppercl 6.2；

　　format lcl2 6.2；format ucl2 6.2；

　　format SE 6.2；format weight percent5.；

　　　format Q1 Q3 4.2；

　if studyname="" then do；grp=2；study2="Overall"；end；

　else do；grp=1；end；

　　ObsId=_N_；

　　if mod（obsid，2）= 0 then studyref=studyname；

　　sum_weight=& sum_weight.；

　　weight=weight0/sum_weight；

　　lcl2=lowercl；

　　ucl2=uppercl；

　if grp=1 then do；

　　weightQ=weight0/sum_weight*0.5；

　　Q1=OddsRatio-OddsRatio*weightQ；

　　Q3=OddsRatio+OddsRatio*weightQ；

　　study2=""；

　　end；

　　OR='OR'；LCL='可信区间低限'；UCL='可信区间上限'；WT='权重'；

run；

data _null_；

　pct=0.75/nobs；height=3；dpi=100；

　call symputx（"pct"，pct）；

　call symputx（"pct2"，2*pct）；

　call symputx（"dpi"，dpi）；

　call symputx（"height"，height）；

　call symputx（"heightin"，height || "in"）；

　call symputx（"thickness"，floor（height*dpi*pct））；

　set forest nobs=nobs；

run；

proc template；

　define style styles.ForestColor93；

　　parent = Styles.htmlBlue；

　　　style GraphFonts from GraphFonts /

　　　　'GraphDataFont' =（"<sans-serif>，<MTsans-serif>"，7pt）

　　　　'GraphValueFont' =（"<sans-serif>，<MTsans-serif>"，7pt）；

　　　style GraphData1 from GraphData1 /

　　　　contrastcolor = GraphColors（'gcdata2'）

　　　　color = GraphColors（'gdata2'）；

　　　style GraphData2 from GraphData2 /

　　　　contrastcolor = GraphColors（'gcdata1'）

　　　　color = GraphColors（'gdata1'）；

　　end；

run；

/*设置森林图样式，DPI，图像参数和标题等*/

ods graphics / reset width=6in height=& heightin imagename="ForestPlotColor" IMAGEFMT=EMF；

ods listing sge=off style=Styles.ForestColor93 image_dpi=& dpi；

title h=12pt "二分类数据的贝叶斯 Meta 分析的森林图"；

title2 h=8pt '比值比 OR 及其 95% 可信区间'；

proc sgplot data=forest noautolegend nocycleattrs；

　refline studyref / lineattrs=（thickness=& thickness）transparency=0.85；

　scatter y=study2 x=oddsratio / markerattrs=（symbol=diamondfilled size=10）；

　highlow y=studyname low=lcl2 high=ucl2 / type=line；

　highlow y=studyname low=q1 high=q3 / type=bar；

　scatter y=studyname x=or / markerchar=oddsratio x2axis；

　scatter y=studyname x=lcl / markerchar=lowercl x2axis；

　scatter y=studyname x=ucl / markerchar=uppercl x2axis；

　scatter y=studyname x=wt / markerchar=weight x2axis；

　scatter y=study2 x=or / markerchar=oddsratio x2axis；

　scatter y=study2 x=lcl / markerchar=lowercl x2axis；

　scatter y=study2 x=ucl / markerchar=uppercl x2axis；

　scatter y=study2 x=wt / markerchar=weight x2axis；

　refline 1 100 / axis=x；

　refline 0.01 0.1 10 / axis=x lineattrs=（pattern=shortdash）transparency=0.5；

　inset '开腹手术组' / position=bottom；

　inset '腹腔镜下手术组'/ position=bottomleft；

　xaxis type=log offsetmin=0 offsetmax=0.35 min=0.01 max=100 minor display=（nolabel）；

　x2axis offsetmin=0.7 display=（noticks nolabel）；

　yaxis display=（noticks nolabel）offsetmin=& pct offsetmax=& pct2 reverse；

run；

Title；

%mend；

%macro forestdata（PostSummaries_data，type）；

%if & type=1 %then %do；

　data temp0；

　　set & PostSummaries_data.；

　　if _N_= 2；

　run；

　%let tau=0；

%end；

%if & type=2 %then %do；

　data temp0；

　　set & PostSummaries_data.；

　　if _N_= 3；

　run；

　data temp_tau；

　　set & PostSummaries_data.；

　　if _N_= 2；

　　keep mean；

　run；

　proc sql noprint；

　　select mean into：tau0 from temp_tau；

　quit；

　%put & tau0.；

　%let tau=& tau0.；

%end；

data tem1；

　　set temp0；

　　keep mean stddev P2_5 P97_5；

　　rename mean=OddsRatio；rename StdDev=SE；

　　rename P2_5=Lowercl；rename P97_5=Uppercl；

run；

data Forest_plot_data；

　set Bayesian_Dichotomous_data；

/*计算效应值 OR 及标准误*/

　OddsRatio=exp（di）；

　Lowercl=round（exp（di-1.96*sqrt（sigmai2）），.01）；

　Uppercl=round（exp（di+1.96*sqrt（sigmai2）），.01）；

　SE=（Uppercl-Lowercl）/（2*1.96）；

　studyname=author；

　weight0=1/（sigmai2+& tau.）；

　keep studyname OddsRatio SE Lowercl Uppercl weight0；

run；

proc sql noprint；

　select sum（weight0）into：sum_weight from Forest_plot_data；

quit；

　%put & sum_weight.；

data forest；

　set Forest_plot_data tem1；

run；

%forestplot（forest）；

%mend；

%forestdata（mcmc_Summaries_fixed，1）；

/*附录代码 5：随机效应模型贝叶斯 Meta 分析*/

ods output PostSummaries =mcmc_Summaries_random；

ods graphics on；

Proc mcmc data=Bayesian_Dichotomous_data nbi=5000 nmc=50000 seed=246810

monitor=（theta tau2 OR）stats（percentage=（2.5 97.5））=all dic diagnostics=all；

parms theta=0 tau2 1；

prior theta ～ normal（mean = 0，var = 1e6）；

prior tau2 ～ igamma（0.001，scale =0.001）；

random deltai ～ normal（mean = theta，var = tau2）subject=study monitor=（deltai）；

model di ～ normal（deltai，var = sigmai2）；

OR=exp（theta）；

run；

ods graphics off；

proc print data=mcmc_Summaries_random；run；

%forestdata（mcmc_Summaries_random，2）；

传统 Meta 分析的统计基础是频率学理论，这种 Meta 分析通常采用 Q 统计量衡量不同研究之间的异质性大小，然后根据异质性检验结果来确定采用何种模型（固定效应模型或随机效应模型）进行 Meta 分析。经典频率学随机效应模型可通过 Q 检验获得研究间方差的矩估计，但是无法获得研究间协方差的 95% 可信区间（confidence interval，CI），而且基于频率学的 Meta 分析检验效能较低。在传统的 Meta 分析模型中，效应量指标的分布一般假设服从正态分布。例如，对于二分类资料的 Meta 分析而言，最常用的效应量为比值比（odds ratio，OR）^[1]。经典 Meta 分析方法是将 OR 对数化，然后将 log（OR）设为效应值，并假设其服从正态分布^[2]。这种 Meta 分析模型是以正态分布假设为前提的，因此在存在小样本资料时，如果不符合近似正态的条件，经典方法合并效应值存在一定误差^[3]。此外，经典 Meta 分析模型很难胜任复杂的统计模型，并且无法识别和处理极端值^[4]。

贝叶斯 Meta 分析（Bayesian meta-analysis，BMA）是近年来以贝叶斯统计理论为基础发展起来的一种新型的 Meta 分析方法，主要采用“马尔科夫链-蒙特卡罗”（Markov chain Monte Carlo，MCMC）方法进行 Meta 分析^[5]。贝叶斯方法与频率论（经典）方法的原理完全不同，频率论认为概率反映的是某事件在重复试验中出现的次数与试验次数之比，也就是事件发生的频率，具有随机性和稳定性特点，这是频率论的统计基础；而贝叶斯方法认为概率是数值的可信程度。这意味着，频率论将模型参数设置为固定，而将数据设置为随机；贝叶斯则与之相反^[6]。贝叶斯统计学与传统频率学在统计推断的理念和方法上存在本质上的区别，包括对信息的利用、对未知参数 θ 的解释及统计推断理念上的差异等^{[4, 7]}。

Meta 分析同样可采用贝叶斯分析实现。贝叶斯统计将参数 θ 作为一个随机变量，有一定的先验分布，在获得样本之后（给定的样本信息），θ 的后验分布 π（θ|x）应包含 θ 的综合信息，可从 θ 的后验分布获得参数 θ 的贝叶斯估计。贝叶斯分析能很好地解决经典 Meta 分析存在的缺陷，近年来，采用贝叶斯方法的 Meta 分析引起了广泛重视，基于贝叶斯的网络 Meta 分析是 Meta 分析领域的研究热点。贝叶斯 Meta 分析最常用的统计软件是 WinBUGS、R 软件，但实际上 SAS 提供了两种过程来实现贝叶斯 Meta 分析，分别是 PROC GENMOD 与 PROC MCMC^[8]。本研究团队先前已介绍了利用 PROC MCMC 实现诊断性试验的贝叶斯 Meta 分析方法^[6]。然而，目前尚无学者介绍利用 SAS 软件实现二分类数据贝叶斯 Meta 分析的方法。本文通过实例数据介绍利用 SAS 软件中的 PROC MCMC 实现二分类数据贝叶斯 Meta 分析，同时将结果与基于频率学 Meta 分析方法的结果进行对比，相关的 SAS 代码在文后附录中给出。

1   贝叶斯 Meta 分析理论简介

1.1   贝叶斯统计理论简介

贝叶斯分析是基于贝叶斯定理发展起来用于系统阐述和解决统计问题的方法^[3]。其理论简述如下：根据既往经验，假设某事件 θ 发生的概率为 P（θ），我们这个概率称为先验概率。然后我们现有一组新的数据 y，则 y 在 y/θ 的前提下发生的条件概率记为 P（y/θ），此条件概率称为似然。根据先验概率和似然可计算出概率 P（θ/y），表示 θ 在 y 存在的前提下发生事件的可能性大小，称为后验概率。后验概率与先验概率与似然的乘积成正比，即 P（θ/y）=α×P（y/θ）×P（θ）。在贝叶斯框架下，参数 θ 是一个随机变量，服从一定的先验分布，在给定的样本信息（抽样）条件下，通过抽样获取样本，则 θ 的后验分布 应包含 θ 的综合信息，因此可通过 θ 的后验分布来获得参数 θ 的贝叶斯估计。因此，贝叶斯统计是由模型、参数和似然 3 个要素构成的统计方法。

1.2   二分类数据贝叶斯 Meta 分析的建模方法

二分类数据 Meta 分析最常用的效应指标是 OR，经过对数转换后的 OR 值服从正态分布。与频数法相对应，本文贝叶斯 Meta 分析的建模方法同样基于 OR 的对数转换值。具体模型构建如下：

假设某 Meta 分析共纳入 n 个研究，nt、rt、nc、rc 分别为治疗组和对照组总人数及事件发生例数。 为第 i 个研究的效应量，本文为 OR 的对数转换值， 为第 i 个研究的效应量的研究内方差，theta 为研究的真实效应值。上述参数建模如下：

固定效应模型假设所有纳入 Meta 分析的研究来自同一个总体，因此效应值固定。对于固定效应模型，，其中 是固定效应模型背景下的真实效应值，是待估计的固定效应模型参数。因此，贝叶斯固定效应模型可建模如下：

，需要获得的数据集包括：；。模型待估计的参数为 。似然函数为正态分布密度 。后验分布为 。

但由于纳入 Meta 分析的原始研究通常存在异质性，导致使用固定效应模型的假设条件往往不成立。这时采用随机效应模型，即假设纳入 Meta 分析的研究可能来自不同的研究总体，每一个研究都拥有不同的效应值，这些效应值服从同一分布的总体。对于随机效应模型，，而，方差 反映了研究间变异的度量，代表随机效应模型的随机部分。贝叶斯随机效应模型可建模如下：

，。需要获得的数据集包括：；。模型待估计的参数包括 ， 和 。其中 表示每个研究的效应， 为真实效应值； 为研究的方差。似然函数与固定效应模型相似。后验分布为 。

2   资料与方法

2.1   示例数据介绍

该示例数据来源于陈艳丽等^[9]发表的 Meta 分析，该 Meta 分析旨在评价腹腔镜下根治性全子宫切除对比开腹手术对早期宫颈癌的疗效和预后的影响。本研究提取该文中并发症的发生率数据（表 1），采用贝叶斯 Meta 分析比较两种手术方式对术后并发症发生率的影响。

将表 1 的数据录入 SAS 软件，并存储在名为 Dichotomous_data 的数据集中，数据集录入的代码见附录。为拟合本文的 BMA 模型，需要的数据集变量包括 author、year、rt、Nt、rc、Nc；分别代表作者、发表年份、研究组事件数（在本示例中，事件代表术后发生并发症）、研究组样本量、对照组事件数和对照组样本量。

2.2   SAS 代码介绍

利用附录代码 1 完成原始数据集的录入之后，再利用附录代码 2 计算效应值 di 及标准误 sigmai2。需要注意的是，本文示例数据不存在零事件试验，对于零事件试验，需要进行连续性校正，校正方法是 2×2 四格表每个单元格增加 0.5；然后使用 PROC MCMC 程序，分别拟合固定效应模型和随机效应模型的贝叶斯 Meta 分析。

在固定效应模型 BMA 中，模型将迭代 50 000 次，前 5 000 次退火用于消除初始值的影响。固定效应值 θ 的先验分布被设置为均数为 0，方差为 1×10⁶的正态分布。不同研究观测的效应值 di 被建模为模型均数为 θ，方差为 1×10⁶的正态分布。与固定效应模型不同的是，随机效应模型引入一个反映研究间变异的度量 tau2，并且 tau2 的先验分布被假设服从形状参数（shape parameter）为 0.001、逆尺度参数为 0.001 的伽玛分布（Gamma distribution）。在随机效应模型背景下，通过“random deltai ～ normal（mean=theta，var=tau2）subject=study monitor=（deltai）；”语句设置随机效应。

为保证基于 SAS 软件的贝叶斯 Meta 分析具有实用性，本文给出了 SAS 软件背景下的森林图绘制代码^[10]。森林图绘制代码将笔者先前发表的 SAS 代码包装成两个宏程序，分别为&forestplot 和&forestdata。前者用于绘制森林图，其宏参数是需要绘制森林图的 SAS 数据集 Forest；后者用于构建森林图数据集，并调用宏程序 forestplot 实现森林图绘制；其宏参数是 PROC MCMC 过程步完成数据分析后输出的后验摘要数据集 PostSummaries_data 和指定的 Meta 分析模型 type（type=1 表示固定效应模型，type=2 表示随机效应模型）。

3   结果

3.1   固定效应模型贝叶斯 Meta 分析结果

固定效应模型贝叶斯 Meta 分析结果见表 2，表中的 theta 和 OR 是对数转换关系，即 OR=exp（theta）。表 3 是 PROC MCMC 提供的参数值的 95% 可信区间。图 1 是 MCMC 采样分析图，分别包括了 Gibbs 抽样动态踪迹图、自相关图和后验分布核密度估计图。从图 1 中的 Gibbs 动态图可看出，参数 OR 经过 50 000 次迭代后基本收敛。在自相关图中，OR 的验后自相关非常小，模型可靠性高。核密度图显示，OR 的后验分布类似于正态分布，与研究预期相似。

3.2   随机效应模型贝叶斯 Meta 分析结果

表 4 给出了随机效应模型贝叶斯 Meta 分析结果。除了给出 θ 的估计值外，随机效应模型还给出了各个研究的随机效应值。表 4 中，除了 OR 值外，其他需要经过拟对数转换后应用于分析。随机效应模型的 MCMC 采样分析图见图 2，除了图 2，随机效应模型同时还会对其他参数进行 MCMC 采样分析，读者可自行运行本文的 SAS 代码查看。

3.3   不同 Meta 分析方法的结果比较

表 5 给出了频率学方法和贝叶斯方法的 Meta 分析结果，两种方法给出的模型参数估计值极其相似。

3.4   贝叶斯 Meta 分析的森林图结果

森林图可非常直观地观察到 Meta 分析的结果。但 PROC MCMC 过程步无法直接给出森林图，为便于读者利用 SAS 软件开展研究，本文开发了 SAS 软件绘制森林图的代码（见附录），运行代码后，可获得固定效应模型和随机效应模型的森林图。图 3 显示了固定效应模型的森林图。

4   讨论

贝叶斯统计和经典统计学派（频率学派）是目前两大主要的统计学派。贝叶斯统计比频率学统计方法需要更强大的计算能力来完成推断，因此贝叶斯统计是以计算机运算能力发展为前提的。随着计算机技术的发展，贝叶斯统计学在医学领域内的应用越来越广泛，其中也包括 Meta 分析。贝叶斯统计是综合未知参数的先验信息和样本信息，根据贝叶斯定理，求出后验分布，然后根据后验分布对未知参数进行统计推断，最后获得参数估计值的统计方法。在方法学上，频率统计学方法以样本为基础，通过假设检验和统计计算的方法进行统计推断，最后获得结论；而贝叶斯方法则是基于贝叶斯原理，根据先验概率推断后验概率^[11]。虽然传统 Meta 分析通常采用频率学方法实现，但越来越多的学者在探索利用贝叶斯方法实现 Meta 分析的方法。目前已有不少统计软件可实现贝叶斯 Meta 分析，包括 R、Stata、WinBugs 等众多软件^{[4, 8, 12-16]}。Winbugs 是用于 Gibbs 抽样的专用软件包，是免费软件，目前已广泛应用于实施贝叶斯方法^[4]。既往学者研究表明，与频率学 Meta 分析相比，贝叶斯 Meta 分析具有以下优点：① 建模方式更灵活，在多种贝叶斯分析软件中可操作性强；② 贝叶斯 Meta 分析充分考虑了模型参数（如方差等）的不确定性，例如可估计研究间方差的不精确性，最终获得的合并效应量点估计及 95% 可信区间来源于参数的后验分布，结果可靠；③ 贝叶斯 Meta 分析提供的统计检验，反映在某种程度上可改变人们的观念^[17]。此外，贝叶斯分析方法最大的优势是不受样本量限制，因为贝叶斯的基本原理就是通过先验概率推断后验概率；而频率学派反映的是抽样概率，因此需要样本量的支持。而实践中往往纳入 Meta 分析的原始研究数量极为有限，这为贝叶斯统计提供了更广阔的应用范围。贝叶斯 Meta 分析不足之处有以下两点：① 频率学派认为一个事件的概率可通过大量重复试验下事件出现的频率来解释，这种概率不依赖于主观性，被认为是客观概率，而贝叶斯统计需要指定先验分布，而这种分布带有主观性，因此被称为主观概率；因此，其参数的先验分布受统计分析师建模思想的影响；② 尽管贝叶斯统计并不接受频率观点，但贝叶斯统计仍然需要依赖样本分布，而样本则来源于抽样。在实际应用中，往往很难完全区分频率学派和贝叶斯统计孰优孰劣。在 Meta 分析领域，初学者往往难以把握统计建模，应以频率学统计结果为准；而具有统计建模基础的系统评价员，则应更多尝试贝叶斯统计方法。

虽然 SAS 软件在统计学领域具有广泛的应用基础，但基于 SAS 软件的贝叶斯 Meta 分析却起步较晚。目前介绍 SAS 软件实现贝叶斯 Meta 分析的文献并不多。在 SAS 9.2 版本之前，SAS 无法实现贝叶斯模型分析，但在 SAS 9.2 之后，SAS 软件提供了多种过程步来实现贝叶斯 Meta 分析，最常用的过程步包括 PROC GENMOD 和 PROC MCMC。PROC GENMOD 可通过“BAYES”语句实现贝叶斯 Meta 分析；但这种贝叶斯分析与 PROC MCMC 的贝叶斯分析基础不同，前者类似于频率学方法，后者是基于 MCMC 的统计方法。因此，在 SAS 中，贝叶斯 Meta 分析首选 PROC MCMC^[18]。国内外已有不少学者介绍了基于 PROC MCMC 实现不同数据的贝叶斯 Meta 分析^{[6, 8, 19]}；但目前尚无介绍 PROC MCMC 实现二分类数据 Meta 分析的文献。鉴于二分类数据在医学领域极为常见，有必要介绍二分类数据的贝叶斯 Meta 分析。曾平等^[20]详细介绍了二分类数据的贝叶斯 Meta 分析的建模方法，本文的侧重介绍该模型在 SAS 软件中的实现过程。

本研究详细介绍了利用 PROC MCMC 实现二分类数据贝叶斯 Meta 分析的 SAS 程序。通过示例数据分析显示，贝叶斯 Meta 分析几乎提供了与频率学 Meta 分析完全一致的结果。但需要注意几点：① PROC MCMC 是以马尔科夫链-蒙特卡罗方法为基础的。MCMC 是在贝叶斯理论框架下，通过计算机进行模拟的蒙特卡洛方法（Monte Carlo）。该方法将马尔科夫（Markov）过程引入到 Monte Carlo 模拟中，实现抽样分布随模拟的进行而改变的动态模拟，弥补了传统的蒙特卡罗积分只能静态模拟的缺陷。② 在多数情况下，Metropolis-Hastings 抽样（M-H sampling）是 MCMC 常用的算法；但对于高维的情形，M-H 算法效率太低。Meta 分析正是一种特殊的高维形式，因此 Gibbs 抽样应用更广泛。③ PROC MCMC 提供了两种不同的后验可信区间（posterior intervals），分别为等尾可信区间（equal-tail intervals，ET）和最高后验密度可信区间（highest posterior density intervals，HPD）。读者在使用可信区间时，应该首选 ET 可信区间，因为 ET 可信区间更窄，而 HPD 可信区间更宽。④ 贝叶斯 Meta 分析不提供统计学概率值（P值），因此不能利用P值来判断效应值是否具有统计学意义，一般来说，贝叶斯通过后验分布获得的可信区间来反映效应值的统计学意义。以 OR 指标为例，根据既往统计学理论，OR 效应值是否具有统计学意义以是否跨越“1”为标准，当可信区间包含了“1”，就认为效应值无统计学意义。因此，结合本文的示例数据，随机效应模型贝叶斯 Meta 分析获得的 OR 值为 0.53，95% 可信区间为（0.30，0.85），该可信区间并未包含 1，因此具有统计学意义。其临床意义为：与开腹手术相比，腹腔镜下宫颈癌根治术的术后并发症发生率明显降低，差异具有统计学意义。最后，与目前多数贝叶斯统计软件相比，SAS 软件编程功能强大。基于 SAS 的编程基础，本研究开发了一段 SAS 宏程序，可非常方便地实现森林图绘制，这有效解决了目前多数贝叶斯统计软件无法获得森林图的问题。

总之，基于 SAS 软件强大的编程能力，PROC MCMC 可非常方便地实现二分类数据的贝叶斯 Meta 分析，随着贝叶斯统计理论的快速发展，贝叶斯 Meta 分析在 Meta 分析领域将发挥越来越重要的作用。

附录代码

/*附录代码 1：数据集建立及资料录入*/

Data Dichotomous_data；

　length author $20.；

　study=_N_；

　input author year rtNtrcNc；

　cards；

Lee2002430930

……

；

Run；

/*Proc print data=Dichotomous_data noobs；run；*/

/*附录代码 2：计算效应值 di 及标准误 sigmai2*/

Data Bayesian_Dichotomous_data；

　set Dichotomous_data；

　Di=log（（rt/（Nt-rt））/（rc/（Nc-rc）））；

　sigmai2=1/rt+1/（Nt-rt）+1/rc+1/（Nc-rc）；

run；

proc print data=Bayesian_Dichotomous_data noobs；run；

/*附录代码 3：固定效应模型贝叶斯 Meta 分析*/

ods output PostSummaries =mcmc_Summaries_fixed；

ods graphics on；

Proc mcmc data=Bayesian_Dichotomous_data nbi=5000 nmc=50000 seed=1234

monitor=（theta OR）stats（percentage=（2.5 97.5））=all dic；

parms theta 0；

prior theta ～ normal（mean = 0，var = 1.0e6）；

model di ～ normal（theta，var = sigmai2）；

OR=exp（theta）；

run；

ods graphics off；

proc print data=mcmc_Summaries_fixed；run；

/*附录代码 4：森林图绘制代码*/

%macro forestplot（forestdata）；

data forest；

　　set & forestdata.；

　　format OddsRatio 6.2；

　　format Lowercl 6.2；format Uppercl 6.2；

　　format lcl2 6.2；format ucl2 6.2；

　　format SE 6.2；format weight percent5.；

　　　format Q1 Q3 4.2；

　if studyname="" then do；grp=2；study2="Overall"；end；

　else do；grp=1；end；

　　ObsId=_N_；

　　if mod（obsid，2）= 0 then studyref=studyname；

　　sum_weight=& sum_weight.；

　　weight=weight0/sum_weight；

　　lcl2=lowercl；

　　ucl2=uppercl；

　if grp=1 then do；

　　weightQ=weight0/sum_weight*0.5；

　　Q1=OddsRatio-OddsRatio*weightQ；

　　Q3=OddsRatio+OddsRatio*weightQ；

　　study2=""；

　　end；

　　OR='OR'；LCL='可信区间低限'；UCL='可信区间上限'；WT='权重'；

run；

data _null_；

　pct=0.75/nobs；height=3；dpi=100；

　call symputx（"pct"，pct）；

　call symputx（"pct2"，2*pct）；

　call symputx（"dpi"，dpi）；

　call symputx（"height"，height）；

　call symputx（"heightin"，height || "in"）；

　call symputx（"thickness"，floor（height*dpi*pct））；

　set forest nobs=nobs；

run；

proc template；

　define style styles.ForestColor93；

　　parent = Styles.htmlBlue；

　　　style GraphFonts from GraphFonts /

　　　　'GraphDataFont' =（"<sans-serif>，<MTsans-serif>"，7pt）

　　　　'GraphValueFont' =（"<sans-serif>，<MTsans-serif>"，7pt）；

　　　style GraphData1 from GraphData1 /

　　　　contrastcolor = GraphColors（'gcdata2'）

　　　　color = GraphColors（'gdata2'）；

　　　style GraphData2 from GraphData2 /

　　　　contrastcolor = GraphColors（'gcdata1'）

　　　　color = GraphColors（'gdata1'）；

　　end；

run；

/*设置森林图样式，DPI，图像参数和标题等*/

ods graphics / reset width=6in height=& heightin imagename="ForestPlotColor" IMAGEFMT=EMF；

ods listing sge=off style=Styles.ForestColor93 image_dpi=& dpi；

title h=12pt "二分类数据的贝叶斯 Meta 分析的森林图"；

title2 h=8pt '比值比 OR 及其 95% 可信区间'；

proc sgplot data=forest noautolegend nocycleattrs；

　refline studyref / lineattrs=（thickness=& thickness）transparency=0.85；

　scatter y=study2 x=oddsratio / markerattrs=（symbol=diamondfilled size=10）；

　highlow y=studyname low=lcl2 high=ucl2 / type=line；

　highlow y=studyname low=q1 high=q3 / type=bar；

　scatter y=studyname x=or / markerchar=oddsratio x2axis；

　scatter y=studyname x=lcl / markerchar=lowercl x2axis；

　scatter y=studyname x=ucl / markerchar=uppercl x2axis；

　scatter y=studyname x=wt / markerchar=weight x2axis；

　scatter y=study2 x=or / markerchar=oddsratio x2axis；

　scatter y=study2 x=lcl / markerchar=lowercl x2axis；

　scatter y=study2 x=ucl / markerchar=uppercl x2axis；

　scatter y=study2 x=wt / markerchar=weight x2axis；

　refline 1 100 / axis=x；

　refline 0.01 0.1 10 / axis=x lineattrs=（pattern=shortdash）transparency=0.5；

　inset '开腹手术组' / position=bottom；

　inset '腹腔镜下手术组'/ position=bottomleft；

　xaxis type=log offsetmin=0 offsetmax=0.35 min=0.01 max=100 minor display=（nolabel）；

　x2axis offsetmin=0.7 display=（noticks nolabel）；

　yaxis display=（noticks nolabel）offsetmin=& pct offsetmax=& pct2 reverse；

run；

Title；

%mend；

%macro forestdata（PostSummaries_data，type）；

%if & type=1 %then %do；

　data temp0；

　　set & PostSummaries_data.；

　　if _N_= 2；

　run；

　%let tau=0；

%end；

%if & type=2 %then %do；

　data temp0；

　　set & PostSummaries_data.；

　　if _N_= 3；

　run；

　data temp_tau；

　　set & PostSummaries_data.；

　　if _N_= 2；

　　keep mean；

　run；

　proc sql noprint；

　　select mean into：tau0 from temp_tau；

　quit；

　%put & tau0.；

　%let tau=& tau0.；

%end；

data tem1；

　　set temp0；

　　keep mean stddev P2_5 P97_5；

　　rename mean=OddsRatio；rename StdDev=SE；

　　rename P2_5=Lowercl；rename P97_5=Uppercl；

run；

data Forest_plot_data；

　set Bayesian_Dichotomous_data；

/*计算效应值 OR 及标准误*/

　OddsRatio=exp（di）；

　Lowercl=round（exp（di-1.96*sqrt（sigmai2）），.01）；

　Uppercl=round（exp（di+1.96*sqrt（sigmai2）），.01）；

　SE=（Uppercl-Lowercl）/（2*1.96）；

　studyname=author；

　weight0=1/（sigmai2+& tau.）；

　keep studyname OddsRatio SE Lowercl Uppercl weight0；

run；

proc sql noprint；

　select sum（weight0）into：sum_weight from Forest_plot_data；

quit；

　%put & sum_weight.；

data forest；

　set Forest_plot_data tem1；

run；

%forestplot（forest）；

%mend；

%forestdata（mcmc_Summaries_fixed，1）；

/*附录代码 5：随机效应模型贝叶斯 Meta 分析*/

ods output PostSummaries =mcmc_Summaries_random；

ods graphics on；

Proc mcmc data=Bayesian_Dichotomous_data nbi=5000 nmc=50000 seed=246810

monitor=（theta tau2 OR）stats（percentage=（2.5 97.5））=all dic diagnostics=all；

parms theta=0 tau2 1；

prior theta ～ normal（mean = 0，var = 1e6）；

prior tau2 ～ igamma（0.001，scale =0.001）；

random deltai ～ normal（mean = theta，var = tau2）subject=study monitor=（deltai）；

model di ～ normal（deltai，var = sigmai2）；

OR=exp（theta）；

run；

ods graphics off；

proc print data=mcmc_Summaries_random；run；

%forestdata（mcmc_Summaries_random，2）；