网状Meta分析系列文章（二）：贝叶斯框架

作者：石清阳；审稿：武珊珊

统计分析一直以来都分为两个学派：频率和贝叶斯学派。它们最主要的区别就是关于被估参数的认知不同，频率学派认为被估参数是固定的，只不过观测实验有限，所以才会去计算假如进行无限次观测实验后的95%置信区间。而贝叶斯学派认为真正的被估参数无法得到也不需要得到，我们只需要基于现有观测数据计算最接近真实情况的被估参数就可以了，因此被估参数是一个随机变量，我们直接计算它的概率分布。通常情况下这两类方法的选择没有优劣之分，但是在这里的网状Meta分析中，贝叶斯方法有一些独特的优势，比如灵活性高、天然的决策模型等等，因此本文主要介绍贝叶斯网状Meta分析模型。

核心模型框架

首先考虑Pairwise Meta分析，假设有M个试验，均为干预1对比干预2。此时，固定效应模型即假设每个试验i均为固定的效应量d12（即干预2相对于干预1的风险比OR或RR），不同试验中观测到的OR仅因随机误差的原因而不同，简而言之，研究间无效应修饰作用（即无异质性）。然而，当存在明显异质性时，通常选用随机效应模型，即假设每个试验i产生一个效应量δi,12，表示为第i个研究中干预2比干预1的OR。随机模型中的效应量并不要求相等，但假设其满足可交换性，即效应量δ和试验标签i完全独立，也就是不同试验i中的效应量δi之间“可交换”(Bernardo and Smith, 1994)。这个假设也可以表示为效应量δi抽样自一个“公共”分布或“随机效应”分布。通常这个分布选择正态分布，即：

其中d12和σ12表示为正态分布的均值和标准差。当然，其他任何符合情况的分布都可以在这里使用。随机效应中的d12其实就是所有研究最终的合并效应量，也是我们所关注的、需要估计的核心参数。当标准差σ12=0时，模型退化为固定效应，此时δi,12=d12。

在贝叶斯框架中，参数估计可以额外施加一个先验分布，即估计前给定一些限制。然而，通常情况下，或在没有足够把握的情况下，我们希望最终估计值可以完全反应所纳入众多RCT的数据情况，因此会给定一些无信息或弱信息先验分布。无信息的意思是先验分布可以覆盖任何可能的参数取值范围，例如我们想估计最终的合并效应量d12，在连续性分布尺度下（分类变量取对数即为连续性尺度），一个合理的可能是正无穷到负无穷，同时越偏离0点概率越小。因此，一个合适的无信息先验分布即为正态分布：

此外，对于研究间异质性参数σ12的先验分布则需要限制在所有正数之中，通常设定为均匀分布：

参数b为分布上限，需要根据具体效应量尺度决定。

具体实例

下面考虑一个假想的例子，假设一共有m个试验（i=1, 2, 3…, m），每个试验均为干预2对比干预1（或治疗组对比对照组）。在第i个试验中，设rik为事件发生数，nik为总人数，k为干预臂数（k=1为对照组，k=2为治疗组）。对于每个事件r只有发生、不发生两种可能，因此假设其满足二项分布：

则pik为第i个试验第k臂中发生事件的概率（i=1,…,m; k=1, 2）。

此时数据为二分类变量（即只能取1发生；0不发生），因此在线性模型中需要变换为连续分布尺度（即取值范围由0、1变为全体实数）。通常我们选择使用logit函数将其映射至实数范围：

logit函数即为log odds函数，因此参数δi表示第i个试验中第k臂对比第1臂的logOR值（δi,12表示治疗组比对照组的logOR值），同时也是我们所关注的核心参数。δi,11表示为对照组比对照组的logOR值，设为δi,11=0（即效应相同）。因此完整的公式表示为：

对照组：

治疗组：

对于随机效应模型，logOR的参数δi,12=logit(pi2)logit(pi1)取自随机效应分布。对于固定效应模型，δi,1k替换为di,1k且固定为常数（也可表示为随机效应分布方差等于0的情况）。

在模型中可以估计两个参数：对照组的log odds（即μ）、治疗组比对照组的logOR（即d或δ）。然而对于合并效应量，我们仅关注其中的logOR值，也就是治疗组比对照组的效应强度比。对照组单独的log odds并不是我们所关注的，属于冗余参数，通常给定一个无信息先验（即标准差100的正态分布）。

下面通过BUGS的代码展示一下固定效应模型的完整过程：

注意在BUGS的代码中，正态分布的参数指定为均值和精确值，精确值为方差的倒数，因此等于0.0001（即1/1002）。

最终我们关注的是d[2]的后验分布，分布的中位数即是治疗组比对照组的logOR，则exp(d[2])即等于OR值。

随机效应模型类似，但需要将d替换为δ，并且为其指定随机效应分布：

展示BUGS代码最主要的目的是后续说明间接比较和网状分析是如何在此框架下被直接拓展的。

间接比较和网状Meta分析

在前述中，我们假设共有m个试验，每个试验均为两个干预间比较，并且每个试验i中的干预2比干预1的效应量δi,12在均值为d12的分布中满足可交换性。现在我们将其拓展，假设同时存在一些试验为干预1和干预3之间的比较，我们也认为干预3比干预1的效应量δi,13在均值为d13的分布中满足可交换性。此时，假设以干预1为节点的干预比较之间满足传递性（transitivity），则有δi,23=δi,13δi,12，进一步同时假设干预3比干预2的效应量δi,23在均值为d23的分布中满足可交换性，即：

因此，可以推导出以下关系（Lu and Ades, 2009）：

注意一点，在pairwise meta分析中，我们分别假设了δi,12或δi,13满足可交换性，而在网状分析中，我们需要δi,12和δi,13同时满足可交换性，才能有下面的推导。我们也可以将其直接拓展为存在多个干预的情形（k=1,2,3…S），即δ在所有试验中均满足可交换性，则接下来传递性才可以确保间接比较δi,xy满足可交换性，因此，进一步推导出如下一致方程组（consistency equations）（Lu and Ades, 2009）：

上述假设是进行网状Meta分析时所特有的假设，然而当我们给定明确的PICOs时，或给定明确的研究问题和纳排标准时，这并不是额外的假设，而可以看作是原有假设的自然拓展。尽管我们依然需要有一些考虑来论证这些假设成立，或使用一些统计检验来补充验证。

对于所有研究间方差，我们通常认为其相等，即σ212=σ213=σ223=…（Higgins and Whitehead, 1996）。在这个假设下，也意味着不同干预比较之间的相关系数等于0.5（即ρ=0.5）。当我们纳入研究中的人群相似、研究设计相似时，此假设是合理的。或存在一些已知的效应修饰时，我们可以在给定一个修饰水平下认为此假设满足。

在网状比较时，因为不同试验中的干预比较不在相同，因此需要额外的标记来区分不同研究i中不同干预臂k中的不同干预。则试验i中第k臂上的干预比第1臂上的干预δi,1k取自一个公共的随机效应分布：

因为δi,1k总意味着在第i个试验中，第k臂比第1臂的效应量，因此可以省略写为δik。同时dt表示第k臂的干预比第1臂的干预的合并效应量，σ2表示研究间方差，则干预1和干预2对比为d12，干预2和干预3对比为d23。如前述，固定效应模型为方差等于0的特殊例子，即δik=dt。

尽管dt可以表示所有的干预比较，但是只有d1k的估计值是我们所关注的，因此我们为其给定无信息先验：

具体实例

同前述，考虑rik表示试验i中第k臂的事件数，nik表示总人数，同样此过程假定为二项分布，并通过logit函数映射至实数域：

其中μi表示对照组的log odds，δik表示第k臂干预比对照的logOR。同样的，尽管μi在不同的试验中表示不同的对照效应量，但仍不在我们的参数关系范围，属于冗余参数。

在固定效应模型中，公式替换为：

在BUGS中，固定效应模型的过程如下：

上述代码和前述pairwise Meta分析的代码几乎一致，除了我们需要额外应对多臂循环na[i]、多干预比较d[t[i,k]]和多个比较d[k]的先验分布。其中t[i,k]表示第i个试验中第k臂上的干预。

在随机效应模型，我们还需要考虑多臂试验中效应量间的内部相关性，并校正此相关性（具体见附录）。

一些注意事项

“干预1”的选择

在前述框架中，我们需要指定一个“干预1”（或对照组）称为reference treatment，以它为基准进行干预比较，那么这个“干预1”应该如何选择？事实上，在模型中，“干预1”的选择并不会影响数据的联合分布本身，然而有时可能会影响到后验分布，这是因为无信息先验分布并非完全无信息，会施加一些非常弱的限制，因此影响后验分布的估计。尽管这个影响可能比较弱，但为了避免这种情况发生，通常可以选择处在网状节点中心的干预作为“干预1”。这样做有一个明显的好处是可以使MCMC的过程更快收敛，当然，无论如何，巨大的迭代次数下模型结果依然是稳健的。

先验分布的选择

在贝叶斯框架下，先验分布的选择通常是一个额外且特有的问题。尽管前文中已经给出了一些先验分布选择的通用方法，比如logOR的估计使用方差很大的正态分布，对于一些特定范围的效应量估计（如均数差MD），可能需要指定一些不同的先验分布。例如一些实验室指标如血糖，取值的范围是正数，此时如果依然使用方差很大的正态分布就不再合适，一个比较简单的代替是使用均匀分布，比如unif(0,30)对于血糖的估计可能已经足够。

对于方差（标准差）的先验分布通常推荐采用均匀分布，比如σ~unif(0,b)，如何决定上限b是一个比较关键的问题。例如，当中位数OR=2（logOR=0.3），则绝大多数logOR可能就在0.3±b这个区间中，即OR=exp(0.3±b)。当我们期望OR的上限在5左右，则b=1.3，若上限在10左右，则b=2，通常情况下OR超过10的可能性还是比较小的。

此外，尽管对于方差的估计采用无信息先验分布已经成为标准法则，但还是存在一些重要的缺陷。在缺少大样本或纳入试验较少时，方差的后验分布无法很好的被估计，通常会很大或很小。因此，这种情况下，使用外部数据进行先验分布的设定是被推荐的，即使用既往发表的相似的试验或超大型Meta分析中的数据（Turner et al., 2012; Rhodes et al., 2015），在这里就不过多论述了。

贝叶斯框架的优势

首先贝叶斯框架下的一个算不上优势的优势就是它的灵活性，这使得模型可以直接拓展到一些其他模型中，比如meta回归、协变量校正、偏倚校正等模型。当然，频率方法也并非不能实现。

事实上，如前述，我们使用了arm-based model而非通常方法中的effect-based model，也就是说对于每一个试验中的每一臂的结果均进行了建模，这样做的好处之一则是可以轻松处理0事件的结局。然而在频率框架中，需要对0事件结局进行“校正”处理，通常是直接加上0.5进行校正，而这会导致一个有偏估计。另一个好处则是可以简单的处理多臂试验的情形。

除此之外，贝叶斯框架的好处在于其“自动”考虑了研究间方差的不确定性，并将其考虑进效应量的后验分布中，形成一个“总”不确定性，即反映在95%可信区间中。然而对于频率框架，则通常简单将研究间方差设为常数，并不做过多的考虑（尽管理论上也是可以实现的）。

总结

简而言之，在贝叶斯框架下，网状Meta分析是Pairwise Meta分析的直接拓展，其所含假设也是在原假设基础上的自然拓展，因此，尽管这些假设依然是需要着重考虑的，但也并非不可论证，我们所需要做的是在每一个网状Meta分析研究开始前，首先考虑我们的研究问题是什么？针对什么人群？考虑哪些干预？是否在同一个临床场景下？这些问题清楚之后，实现网状Meta分析反而可能是研究中最简单的一部分。

附1：多臂试验

在前述框架中，多臂试验已经自然的被考虑进模型中，因此，在固定效应模型中无需任何额外的设定。然而，在随机效应模型，如何指定随机效应分布是一个单独的问题（Lu and Ades, 2004）。在可交换性的假定下，单个多臂试验中会产生多个效应量，如3臂试验产生2个效应量、4臂产生3个，同时这些效应量之间存在内部相关性，因此需要将其考虑进模型中。当我们假设其等方差时，有：

其中δi表示效应量的向量（即多个效应量），其满足多元正态分布。其后在多元正态分布中导出每一臂效应量的条件分布如下（Raiffa and Schlaiffer, 2000）：

附2：随机效应模型的BUGS过程

参考文献：

Bernardo, J. M. & Smith, A. F. M. 1994. Bayesian Theory, New York, John Wiley & Sons, Inc.

Dias, S., Ades, A. E., Welton, N. J., Jansen, J. P., & Sutton, A. J. (2018). Network meta-analysis for decision-making. John Wiley & Sons.

Higgins, J. P. T. & Whitehead, A. 1996. Borrowing strength from external trials in a meta-analysis. Statistics in Medicine, 15, 2733–2749.

Lu, G. & Ades, A. 2004. Combination of direct and indirect evidence in mixed treatment comparisons. Statistics in Medicine, 23, 3105–3124.

Lu, G. & Ades, A. E. 2009. Modelling between-trial variance structure in mixed treatment comparisons. Biostatistics, 10, 792–805.

Raiffa, H. & Schlaiffer, R. 2000. Applied Statistical Decision Theory, New York, Wiley Interscience.

Rhodes, K. M., Turner, R. M. & Higgins, J. P. T. 2015. Predictive distributions were developed for the extent of heterogeneity in meta-analyses of continuous outcome data. Journal of Clinical Epidemiology, 68, 52–60.