第一章
贝叶斯Meta分析基础知识

Meta分析可以定量、科学地合成研究结果，已在许多科学领域取得革命性的成果，在医学领域有助于建立循证实践、解决相互矛盾的研究结果问题。贝叶斯方法正被频繁地应用于卫生保健研究领域，包括Meta分析。本章将主要就经典的贝叶斯Meta分析方法基本知识作简要概述。

第一节　贝叶斯方法概述

一、贝叶斯起源与发展

贝叶斯统计学是以英国数学家Thomas Bayes的名字命名的。其生前著有《论有关机遇问题的求解》一文，该论文给出贝叶斯公式并采用归纳推理的方法以二项分布为基础对未知参数进行推断，但可能由于贝叶斯认为该理论尚存在不完善的地方，此论文在其生前并未发表。1763年，其朋友Richard Price于英国皇家学会宣读了该论文，该论文正式发表。1814年Laplace PC出版了《关于概率的哲学评述》重新推导了贝叶斯公式，并推广获得一系列新成果。此后很长一段时间，因其观点有悖经典统计学，同时因计算困难，长期未被普遍接受。现代贝叶斯统计方法的发展归功于Jeffreys（1939）、Savage（1954）、Ravaffa、Schleifer（1961）、Lindly（1972）及De Finetti（1974）等贝叶斯统计先驱的重要贡献。经后来统计学者的不断努力，发展为一种系统的统计推断方法，称为贝叶斯统计方法。采用这种方法作统计推断所得的全部结果，构成贝叶斯统计的内容。贝叶斯统计需要求解高维函数的积分，在20世纪80年代之前求解这些积分成为其发展的障碍，因此该方法一直停留在理论阶段。20世纪90年代起马尔可夫链蒙特卡罗（Markov Chain Monte Carlo，MCMC）方法被广泛应用于贝叶斯统计，成功地解决了限制贝叶斯统计发展的高维积分运算问题，为贝叶斯统计带来了革命性的突破，进而使贝叶斯统计在更多的领域得到应用。近期发展的哈密尔顿蒙特卡罗（Hamiltonian Monte Carlo，HMC）算法、积分嵌套拉普拉斯近似方法，丰富了贝叶斯计算方法，进一步促进了贝叶斯统计的发展。

二、贝叶斯统计学与经典统计学的区别

贝叶斯统计学与经典统计学在统计推断的理念和方法上有着本质的区别，主要可总结为以下三个方面：

第一，二者对信息的利用有本质差异。贝叶斯统计学派认为在观察到样本之前，对于任一未知的参数 θ 有一定的了解，即已经积累一些关于参数 θ 的信息——先验信息（prior information），在对未知参数进行统计推断时应综合总体信息、先验信息及样本信息。用统计学语言可描述为： θ 作为一个随机变量，有一定的先验分布，其分布为 θ ～ π （ θ ）。在获得样本之后（给定的样本信息）， θ 的后验分布 π （ θ | x ）应包含 θ 的综合信息，关于参数 θ 的统计推断均基于 θ 的后验分布（posterior distribution）进行。因此，贝叶斯统计方法的关键在于所作出的任何推断完全取决于后验分布 π （ θ | x ），而不再涉及样本 x 的分布。参数 θ 是否为随机变量、先验分布是否存在及如何选取，成为经典频率学派集中批评的靶点。频率学派的观点认为参数 θ 在每一个确定的问题中都有一个确定值，无随机性，因而也无分布而言。认为贝叶斯学派以主观概率的立场出发，引进先验分布，将先验分布看作为主观随意性的概念，进而认为贝叶斯统计问题的解也为主观随意性的解，因此无科学意义。迄今为止，贝叶斯学派仍未提出一种确定先验分布的方法，这成为其重大的弱点。但是也应看到，虽然贝叶斯学派采用主观概率的概念，并不是完全主观随意的选取先验分布，而是以从实践中获得的主观认识作为先验信息。虽然理论上尚无统一的、完整的、不失一般性的确定先验分布的方法，但是在实用范围内的常见问题所采用的先验分布，已经得到正确的验证。

第二，对未知参数 θ 解释的差异。频率学派对参数 θ 解释是概率的频率解释，尤其是对置信区间的求解和解释。经典统计学首先指定置信水平（1- α ），然后构建一个含有未知参数 θ 的枢轴量，通过枢轴量的分布求得参数 θ 置信水平为（1- α ）的置信区间。对于经典统计学来讲，参数是固定的而未知，无统计分布而言，因此对于所得到的置信区间的理解存在一定困难。正确的理解为若反复抽样多次，每个样本值可以确定一个区间，这个区间要么包含参数 θ ，要么不包含参数 θ ，在这么多的区间中，包含参数 θ 的占95%，未包含参数 θ 占5%，而不能理解为“有95%的概率使得参数 θ 落在置信区间中”，因为由经典统计学派求得的置信区间已不是随机区间。而贝叶斯学派认为参数 θ 为一随机变量，结合样本信息和先验信息可以构造一个区间，贝叶斯统计中称为可信区间（credit interval），应区别于经典统计学派的置信区间（confidence interval），使得未知参数 θ 以一定的概率落在这个区间中，因此贝叶斯学派可以陈述为“有95%的概率使得参数 θ 落在可信区间中”。

第三，统计推断的理念的差异。经典学派奠基人Fisher将经典统计推断总结为三个问题：选定模型、确定统计量和相应统计量的分布。即选定模型，构建一个分布已知含有未知参数的枢轴量，根据抽样分布来计算统计量的全部性质。而贝叶斯统计推断完全源于未知参数的后验分布，未知参数的所有的统计学性质均由后验分布决定。

第二节　贝叶斯统计学基础

一、贝叶斯公式与贝叶斯定理

贝叶斯统计学的基础是贝叶斯公式和贝叶斯定理。贝叶斯公式可由条件概率的定义和全概率公式推导而得，根据未知参数及样本的不同，贝叶斯公式有三种不同的形式：

（一）贝叶斯公式的事件形式

设试验 E 的样本空间为 S ， A 为 E 的事件， B ₁ ， B ₂ ，… B _n 为样本空间 S 的一个划分，且 P （ A ）>0， P （ B _i ）>0（ i =1，2，…， n ），则由条件概率的定义和全概率公式可得： ， i =1，2，…， n ，该式称为贝叶斯公式的事件形式。

（二）贝叶斯公式的密度函数形式

设 x =（ x ₁ ， x ₂ ，…， x _n ）是来自某总体的一个样本，该总体的概率密度函数为 p （ x θ ）， θ =（ θ ₁ ， θ ₂ ，…， θ _k ），当给定一组观察值 x =（ x ₁ ， x ₂ ，…， x _n ）， θ 的条件概率分布，即后验分布为： ，该式即为贝叶斯公式的密度函数形式，即在样本 x =（ x ₁ ， x ₂ ，…， x _n ）下 θ 的后验分布，式中 π （ θ ）为参数 θ 的先验分布； 为样本 x =（ x ₁ ， x ₂ ，…， x _n ）的联合条件密度函数，即似然函数。

（三）贝叶斯公式的离散形式

当 θ 为离散型随机变量时，先验分布可以用分布列 π （ θ _i ）， i =1，2，…， n 表示，此时 θ 的后验分布也为离散形式： 。

当给定 x =（ x ₁ ， x ₂ ，…， x _n ）时，由于 p （ x ）不依赖于 θ ，其在计算 θ 的后验分布中仅起到正则化因子的作用，而且利用概率分布的正则性可以方便地求出这个因子。根据似然原理，如果两个似然函数成比例，且该比例常数与 θ 无关，则这两个似然函数所包含的关于 θ 的信息是相同，因此将 p （ x ）省略，则 θ 的后验分布可以表示为如下形式： ，该式右侧 p （ x θ ） π （ θ ）称为后验分布的核（kernel）。

二、似然函数和极大似然估计

似然函数是经典统计学和贝叶斯统计联系的纽带，极大似然估计是经典统计学的重要估计方法。若总体密度为 f （ x ； θ ）（当 X 为离散型 f （ x ； θ ）为分布列）， θ =（ θ ₁ ， θ ₂ ，…， θ _n ）为待估参数， Θ 为可能取值范围。设 x ₁ ， x ₂ ，…， x _n 是来自 X 的样本的一组样本值，则称 L （ θ ）= L （ x ₁ ， x ₂ ，…， x _n ； θ ）= 为样本的似然函数。若存在某个 ，使得 成立，则称 （ X ₁ ， X ₂ ，…， X _n ）为极大似然估计量， （ x ₁ ， x ₂ ，…， x _n ）为极大似然估计值。极大似然估计即是使得似然函数取得极大值的参数的取值，通常先将似然函数取对数，然后求导，使导函数等于0，求得参数的估计值。

［例1.1］ 设 X ～ N （ μ ， σ ² ）， μ ， σ ² 为未知参数， x ₁ ， x ₂ ，…， x _n 是来自 X 的一个样本值，求 μ ， σ ² 的最大似然估计值。

［解］

X 的概率密度函数为 ，则似然函数为 L （ μ ， σ ² ） 。对数似然函数为 ，分别对 μ ， σ ² 偏导数，并令其等 ，可得 μ ， σ ² 的极大似然估计值， 。

三、先验分布的选取

先验分布是贝叶斯统计的重要部分，也是最受经典统计学派批评的一点。先验信息主要指经验和历史资料。因此如何使用经验和过去的历史资料确定概率和先验分布是贝叶斯统计要研究的重要内容。统计学家提出多种方法，但至今仍未提出一种放之四海而皆准的确定后验分布的方法。参数的先验分布选取有多种方法，如利用先验信息确定先验分布、利用边缘分布来确定先验分布、无信息先验分布、共轭先验分布、Reference先验、最大熵先验等。贝叶斯Meta分析中最常用的先验分布为无信息先验分布和共轭先验分布，本节将对其做简要介绍。

（一）无信息先验分布

参数 θ 的无信息先验分布指除了参数 θ 的取值范围Θ和 θ 在总体中的地位之外，再无包含 θ 的任何信息的先验分布。无信息先验包括贝叶斯假设、位置参数的无信息先验（参数在平移变换群下的不变性）、尺度参数的无信息先验（参数在刻度变换群下的不变性）及一般情形下的无信息先验（Jeffreys先验）等。

贝叶斯假设对任何取值都没有偏爱，因此自然地将“均匀分布”作为先验分布，此即为贝叶斯假设。贝叶斯假设有三种情况，分别为离散均匀分布、有限区间上的均匀分布及广义先验分布。用数学公式表示为：

在贝叶斯假设下，似然函数 l （ θ x）即为后验密度的核，即 π （ θ x ）∝L（ x ）×1=L（ x）

［例1.2］ 设 X =（ X ₁ ， X ₂ ，… X _n ）为从 N （ θ ，1）总体中抽取的随机样本，若 θ 的先验密度为 π （ θ ）≡ c ，求 θ 的后验密度。

［解］ 由贝叶斯公式及贝叶斯假设可得：

结果显示参数 θ 的先验分布取贝叶斯假设时， θ 的后验分布为 N （ ，1 /n ），结果同极大似然估计，因此可以认为经典统计学是贝叶斯统计学取贝叶斯假设时的特殊情况。

Jeffreys对先验分布的确定做出了重大贡献，利用Fisher信息矩阵给出了确定无信息先验分布的一般方法。Jeffreys原则包括两部分：一是对先验分布应有一个合理的要求；二是给出一个具体的方法去求得合乎要求的先验分布。Jeffreys利用了Fisher信息矩阵的一个不变性质，发现 θ 的无信息先验分布应以信息阵 I （ θ ）的行列式的平方根为核。设 x =（ x ₁ ， x ₂ ，…， x _n ）是来自某总体的一个样本，该总体的概率密度函数为 p （ x θ ）， θ =（ θ ₁ ， θ ₂ ，…， θ _p ），为 p 维参数向量。首先求对数似然函数：

然后求Fisher信息矩阵： ，

则 θ 的无信息先验密度为： ，其中det I （ θ ）为Fisher信息矩阵 I （ θ ）的 p × p 阶行列式。

［例1.3］ 设 θ 为 Bernoulli 试验中的成功概率，则在 n 次对立的 Bernoulli 试验中，成功次数 x ～ Bernoulli （ n ， θ ），即 ，（ x =0，1，2…， n ），求 θ 的Jeffreys无信息先验分布。

［解］ ，（ x =0，1，2…， n ）的对数似然函数为 x ln θ +（ n - x ）ln（1- θ ），故有：

取 π （ θ ）∝ I （ θ ） ^{1 / 2} = θ ^{-1 / 2} （1- θ ） ^{-1 / 2} ，式中（0< θ <1），添加正则化因子即可得到先验密度函数，该先验分布恰好是贝塔分布 Beta （0.5，0.5）。

（二）共轭先验分布

设 θ 是总体分布中的参数， π （ θ ）是 θ 的先验密度函数，设F表示由 θ 的先验分布 π （ θ ）构成分布族。如果对任取 π ∈F及样本值 x ，后验分布 π （ θ x ）仍属于F，则称F是一个共轭先验分布族，也即 π （ θ ）与 π （ θ x ）属于同一类分布族。共轭先验分布有两个优点：①计算简便；②后验分布中的一些参数可以得到很好的解释。虽然共轭先验分布计算简便，但应以合理性为首要原则。共轭先验分布的计算见下文。常用的共轭先验分布见表1-1。

［例1.4］ 设 x ～ Bernoulli （ n ， θ ），若取 θ 的先验分布为 Beta （ α ， β ），证明 θ 的后验分布仍为贝塔分布，即 θ 的共轭先验分布为贝塔分布。

［解］ 若 θ ～ Beta （ α ， β ），则 ，（ x =0，1，2…， n ）由贝叶斯公式可得：

因此，参数 θ 的后验分布为 π （ θ x ）～ Beta （ x + α ， n - x + β ），即参数为（ x + α ， n - x + β ）的贝塔分布，所以 θ 的共轭先验分布为贝塔分布。共轭先验分布具有很好的可解释性， Beta （ α ， β ）分布的期望及方差分别为：

则在共轭先验分布下 θ 的后验分布数学期望（后验均值）为：

对上式进行等价变换，可以对参数 θ 的贝叶斯估计进行解释：

其中 可认为仅采用先验信息对参数 θ 的估计，而 为利用样本信息对参数 θ 的估计， γ 为权重因子，故可以认为 θ 的后验分布是上述两个估计的加权平均。

四、贝叶斯点估计及区间估计

贝叶斯包括三种常用的点估计，分别为后验众数估计（posterior mode estimator）、后验中位数（posterior median estimator）及后验期望（posterior expectation estimator）。使后验密度函数 π （ θ x ）取得最大值的 θ ，称为后验众数估计，记为： 满足 ；后验期望（后验均值）定义如下， π （ θ x ） dθ 。

［例1.5］ 从一批某疾病患者中随机抽取 n 人，其中有 x 人治疗无效，求治疗无效率 θ 的后验众数期望和后验期望估计。已知 Beta （ α ， β ）分布的众数及均值为 。

［解］ 由题意可知， X ～ Binomial （ n ， θ ），取 θ 的先验分布为 Beta （ α ， β ），则 θ 的后验分布为 ，由贝塔分布的众数及期望值公式可得， θ 的后验分布众数估计为 ；后验期望为 。

当 θ 的先验分布取 Beta （1，1），即是均匀分布 U （0，1）时， ，此时后验分布众数估计即为极大似然估计。

设参数 θ 的后验分布为 π （ θ x ），对于给定的样本 x 和概率1- α （0< α <1，通常 α 取较小的数），若存在两个统计量 ，满足 ，则称 为参数 θ 的可信水平为1- α 的可信区间（credible interval）。贝叶斯统计可以解释为参数 θ 属于上述区间的概率为1- α ，这一点应区别于置信区间（confidence interval）的解释。

设 θ 是总体 X 的待估参数， X ₁ ， X ₂ ，… X _n 是取自总体 X 的样本，对给定值0< α <1，若统计量 满足 ，则称随机区间 为 θ 的置信水平为1- α 的双侧置信区间， 分别称为置信下限和置信上限。 是一个随机区间，其端点及长度均为统计量，在获得 X 的一组观测值后，即可获得一个置信区间；由于样本的随机性，因此这个区间不一定包含 θ 的真值，包含真值 θ 的概率是1- α 。置信区间的求解，通常是构建一个枢轴量，根据枢轴量的分布来求得未知参数的置信区间。贝叶斯区间估计相比经典统计学的置信区间具有处理方便及含义清晰的特点。

贝叶斯可信区间估计 并不是唯一的，常用的可信区间为等尾可信区间，计算过程为取参数 θ 的后验分布 π （ θ x ）的 α/ 2和1- α/ 2分位数 θ _L 和 θ _U 作为可信区间的上下限。但等尾可信区间有时并不理想，尤其是对于偏态分布的数据。理想的可信区间的长度是最短的，把具有最大后验密度的点包含在区间内，而区间外的点后验密度值不超过区间内的后验密度值，这样的区间即为最大后验密度可信区间（highest posterior density，HPD），定义如下：设参数 θ 的后验分布为 π （ θ x ），对给定的概率1- α （0< α <1），若存在集合 C 满足如下条件：

（1） P （ θ ∈ C x ）=1- α

（2）对任意的 θ ₁ ∈ C 和 θ ₂ ∉ C ，总有

则称集合 C 为可信水平为1- α 的最大后验密度可信集。

［例1.6］ 不同先验分布对后验分布估计的影响。一项研究调查三年内医院早产婴儿的存活率，假设39个25周出生的婴儿，其中31个存活了至少6个月，问题：①采用（0，1）均匀分布或 Beta （1，1）为先验分布；②采用Jeffreys无信息先验分布；③设既往存活率为60%。计算三种不同先验分布下该新药的有效率的贝叶斯估计？

［解］ 根据前文可知，率的共轭先验分布为贝塔分布， Beta （1，1）即为（0，1）的均匀分布，即 Beta （1，1）= U （0，1）； Beta （0.5，0.5）为Jeffreys先验分布；既往存活率为60%，可选用 Beta （7，5）分布，此时为有信息先验分布。下述命令为自行编写计算二项分布后验分布的代码：

分别采用如下命令计算三种情况的存活率。

三个不同先验分布获得的后验分布主要结果如表1-2所示。

采用LearnBayes包的triplot命令绘制不同先验分布对应的后验分布曲线，如图1-1；命令如下：

结果解读：当取均匀分布及Jeffreys先验分布时，二者均为无信息先验，后验分布的贝叶斯估计非常相近，并且结果接近频率学派估计结果，说明频率学派是贝叶斯统计取无信息先验分布的特例；当取信息先验分布时，后验分布为先验分布与似然函数的折中。

第三节　贝叶斯统计的常见算法

一、MCMC算法

贝叶斯统计基本理论和方法虽然简单易懂，但是在实际应用中由于后验分布涉及多重积分，无法给出解析解及数值解，因此在20世纪90年代前贝叶斯统计多数模型停留在理论研究阶段，只有少数模型可以采用共轭先验分布及近似解获得结果。20世纪90年代后随着随机模拟技术及计算机技术的发展，贝叶斯统计领域被广泛应用，使得贝叶斯统计计算得到迅猛发展，尤其是马尔可夫链蒙特卡罗模拟方法的引进。

马尔可夫链蒙特卡罗模拟（Markov Chain Monte Carlo Simulation，MCMC），由两部分内容组成，即采用马尔可夫链（Markov Chain，MC）获得目标分布的抽样样本，然后通过蒙特卡罗（Monte Carlo，MC）积分获得目标分布的相应结果。

蒙特卡罗方法由美国在第二次世界大战中研制原子弹的“曼哈顿计划”的成员塔尼斯拉夫·乌拉姆和约翰·冯·诺依曼于20世纪40年代首先提出。数学家约翰·冯·诺伊曼用驰名世界的赌城——摩纳哥的Monte Carlo的名字来命名这种方法。蒙特卡罗方法是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。蒙特卡罗方法的解题过程可以归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。具体如下：

若随机变量 X 的密度函数为 f （ x ），假设 g 为一个可积函数，则随机变量 Y = g （ X ）的期望为： 。

如果从 X 的分布中获得一些随机数，则 E ［ g （ X ）］的无偏估计即为相应的样本均值。通俗地理解即当一个积分求解困难时，我们可以通过构建一个函数的期望，从随机变量所服从分布中抽样，计算每个样本值的函数值，求得均值即为积分值。如为估计 ，若 X ₁ … X _N 为均匀分布（0，1）的抽样样本，则由强大数定律可知 ，以概率1收敛到期望 E ［ g （ X ）］。

［例1.7］ 采用蒙特卡罗方法计算 。提示该积分解析解为 。

［解］ （0，1）的均匀分布的密度函数为 f （ x ）=1，从（0，1）的均匀分布中抽取 n =10 000样本，则有 ，可由R语言编写简单代码实现，如：

计算结果为theta.hat=0.631 417，即 ，而积分值 ，两者计算结果相近。

蒙特卡罗方法将积分问题转化为从目标概率密度 f （ x ）中产生随机样本，在实际计算中，后验分布较为复杂，无法直接抽样，必须解决如何方便得到各种复杂概率分布的对应采样样本集的问题，否则仍无法计算蒙特卡罗积分。马尔可夫链为解决抽样问题提供了很好的办法。在MCMC方法中，首先建立一个马尔可夫链，运行足够长时间后，使其收敛于平稳分布，产生的随机样本就是从达到平稳状态的马尔可夫链中产生的样本。

假设序列状态是…， X _{t -2} ， X _{t -1} ， X _t ， X _{t +1} ，…，那么在时刻 X _{t +1} 的状态的条件概率仅仅依赖于时刻 X _t ，即： 。

具有此性质的序列状态称为马尔可夫链性质。不可约、非周期、正常返的马尔可夫链在运行充分长时间后将收敛于平稳分布，与初始分布无关。马尔可夫链的实值函数的遍历均值依概率1收敛于极限分布下的均值。遍历均值做合适的变化将依分布收敛于标准正态分布。这三条性质从理论上保证了MCMC算法的有效性。马尔可夫链的理论基础较为复杂，读者有兴趣可参考相关著作。MCMC算法通过构建马尔可夫链，涉及一大类抽样方法，包括Metropolis抽样、Metropolis-Hastings抽样、Gibbs抽样、随机游走抽样等。最为常见的是Metrapolis-Hastings抽样和Gibbs抽样。

Metropolis-Hastings抽样简称M-H抽样，这个算法首先由Metropolis提出，后被Hastings改进，因此被称之为Metropolis-Hastings抽样，主要是解决Metrapolis抽样效率低的问题。具体抽样算法如下：

1.构造合适的建议分布 g （· X _t ），该分布满足不可约、正常返、非周期的正则化条件；

2.从分布 g （· X _t ）中产生 X ₀ ；

3.重复下列过程，直至马尔可夫链达到平稳状态：

（1）从 g （· X _t ）产生Y；

（2）从均匀分布U（0，1）中产生U；

（3）若 ，则接受 Y 并令 X _{t +1} = Y ，否则 X _{t +1} = X ；接受概率为 α （ X _t ， Y ）=

（4）增加 t ，返回到（1）。

M-H抽样虽然很好地解决了蒙特卡罗方法需要的任意概率分布的样本集问题。但是M-H抽样有两个缺点：一是需要计算接受率，在高维时计算量大，计算接受率导致算法收敛时间变长。二是有些高维数据，特征的条件概率分布好求，但是特征的联合分布不好求。Gibbs抽样很好地解决了M-H抽样的缺点。Gibbs抽样是M-H抽样的一个特例，重新构建了细致平衡条件，使得抽样的接受率等于1，即保留了所有的抽样样本。Gibbs抽样是基于条件概率分布进行的，其基本做法是从联合概率分布定义条件概率分布，依次对条件概率分布进行抽样得到样本的序列。可以证明这样的抽样过程是在一个马尔可夫链上的随机游走，每一个样本对应着马尔可夫链的状态，平稳分布就是目标的联合分布。预烧期之后的样本就是联合分布的随机样本。具体抽样算法如下：

1.随机初始化｛ x _i ： i =1，…， n ｝；

2.对 t =0，1，2…循环抽样：

（1）

（2）

（3）…

（4）

（5）…

（6）

Gibbs抽样大大提高了抽样的效率，成为众多软件的算法基础，如WinBUGS、JAGS等。下面简单演示正态分布参数的Gibbs抽样过程。

［例1.8］ 假设 X =（ X ₁ ， X ₂ ，… X _n ）为从 N （ θ ， σ ² ）总体中抽取的随机样本。取 θ 的先验分布为 N （ μ ， τ ² ）， σ ² 先验分布为 IG （ a ， b ），其中 μ ， σ ² ， a ， b 已知，采用Gibbs抽样对 θ ， σ ² 进行抽样。

［解］ （ θ ， σ ² ）的后验分布核为

化简得到参数的条件分布核如下：

继而化简可知：

由上述结果可知，参数 θ 的条件分布为含有 σ ² 的正态分布，而 σ ² 的条件分布为含有 θ 的 IG 分布；在给定初始值后，Gibbs抽样过程示意图如图1-2：

［例1.9］ 以Mackowiak等（1992）的数据为例，该数据记录了100个人的体温（华氏）、性别和每分钟的心跳数。实验的目的是检验Carl Wunderlich的观点——健康成年人的平均体温为37℃（98.6℉），试采用Gibbs抽样进行贝叶斯估计。

［解］ 采用N（0，100）及IG（0.01，0.01）弱信息先验分布。R语言代码如下：

体温均数、标准差抽样分布及迭代100路径分布图如图1-3：

二、MCMC抽样的收敛性诊断

MCMC方法通过构造一个概率转移矩阵，建立一个以分布π（ x ）为平稳分布的马尔可夫链来得到π（ x ）的样本，通过随机抽样得到的这些样本进行各种统计推断。马尔可夫链的极限分布是目标分布，但是如何判断马尔可夫链收敛于目标分布？这就涉及MCMC的收敛，目前尚无准确无误判断收敛的方法，但是有可以辅助判断是否收敛的手段。所谓算法的收敛性指的是所得到的马尔可夫链是否达到了平稳状态。如果达到平稳状态，则认为所获得的样本来自目标分布。一般情况下，并不知道算法需要运行多长时间才能到达平稳状态，因此对马尔可夫链收敛性进行监视。

MCMC产生的马尔可夫链的收敛性受到三个因素的影响：①初始值：用来初始化马尔可夫链。需要利用一切可得到的信息尽量使得初始值离参数真值更近。如果初始值远离后验密度的最高区域，且算法的迭代次数不足以消除初始值的影响，则会对后验推断产生影响。②后验分布密度：如果后验密度的形状是单峰的，问题相对简单。如果后验密度的形状是多峰的，则要小心伪收敛现象，也就是可能会陷入非最大值的其他极值点。③工具分布：好的工具分布的中心区域应当与后验密度的中心区域有较多的重合，尾部比较厚，支撑包含了后验密度的支撑。

收敛涉及相关的几个术语：①预烧期（burn-in period）：在MCMC机制中，用以保证马尔可夫链达到平稳状态所运行的时间称为预烧期。一般只要链运行时间足够长，预烧期对后验推断几乎没有影响。②筛选间隔或抽样步长（thinning interval or sampling lag）：马尔可夫链产生的样本并不是相互对立的，因此采用每间隔L个样本来抽取一个（近似）独立的样本。③蒙特卡罗误差（MC error）：其度量了估计因随机波动而导致的波动性。在计算感兴趣的参数时，其精度随着样本容量递增，与样本容量成反比。有组平均方法及窗口估计量计算方法。

收敛性的判断主要包括两大类方法，一是图示法，包括踪迹图、遍历均值图、自相关图、核密度估计图；二是收敛性指标，包括Gelman-Rubin法、Geweke法、Raftery-Lewis法、Heidelberg-Welch法等。

踪迹图本质是一个时间序列图，收敛的马尔可夫链应当是没有趋势与周期的，它会在一个水平线上下小幅度波动（形状上类似白噪声序列）。当马尔可夫链收敛时遍历均值将接近一条直线，只有小幅波动。马尔可夫链收敛后，抽样样本量越大，核密度估计图越光滑。下面继续以体温数据为例，绘制上述图形，如图1-4和图1-5，并判断马尔可夫链的收敛性。R语言代码如下：

由踪迹图可知，theta及sigmasqr围绕均值两端均匀分布；核密度估计图提示二者曲线光滑，theta呈现对称分布，sigmasqr呈左偏态分布；自相关图提示相关性低，最大自相关系数为0.01左右；遍历均值图提示，在迭代早期，遍历均值很快达到稳定状态。以上结果提示抽样分布收敛性好。

Gelman-Rubin法由Gelman及Rubin提出来，又称为方差比法。基本思想是当多条马尔可夫链平稳收敛后，马尔可夫链的统计特征应该是一样的，利用类似方差分析的方法构造出一个统计量诊断马尔可夫链的收敛性，分别计算链间的方差和链内的方差，然后计算二者的加权组合来估计方差，最后由加权方差估计量及链内方差估计量计算获得潜在尺度缩减因子（potential scale reduction factor，PSRF），简写为 R 。一般而言， R <1.1提示马尔可夫链收敛。具体计算过程可参考相关著作。R语言coda包可以给出Gelman-Rubin方法的结果。首先生成两组2 000个标准正态分布样本作为两条链；生成一组2 000个 N （2，1）正态分布样本为一条链，进行Gelman-Rubin检验，如图1-6、图1-7。具体如下：

X1由两组2 000个标准正态分布样本组成，可视为同一目标分布产生的两条链；X2由两组不同正态分布样本，可视为没有收敛到目标分布。由图1-7可知X1的两条马尔可夫链混杂在一起，均匀分布在0两侧，同时核密度估计图示两条链的样本呈现一条光滑的曲线。而X2的两条马尔可夫链分别来自 N （0，1）和 N （3，1），因此踪迹图并未混在一起；核密度估计图也呈现出双峰表现。由Gelman-Rubin图可知，X1在迭代500次后， R ≈1；而X2的 R ≈4；说明X2两条链并未收敛。

Geweke检验法，包括了数量诊断与图形诊断两个部分。这种方法的基本思想是：如果马尔可夫链是收敛的，那么这个马尔可夫链前一部分与后一部分的均值应该是相等的。在马尔可夫链前一部分与后一部分渐进独立的假设条件下，Geweke构造的检验统计量渐进服从标准正态分布。数量诊断中，Geweke检验对每个分量马尔可夫链算出检验统计量的值，当值小于1.96或2的时候，马尔可夫链被认为是收敛的。图形诊断中，采用同样的方法算出一个统计量的值，然后，把最前边的一小段马尔可夫链切除，再用同样的方法算出第二个统计量的值，以此类推。一共算出了若干个统计量值之后，最后在坐标系上画出他们的位置，同时用虚线画出某一置信度的置信区间（一般认为是0.95）形成的置信带。如果绝大多数统计量值在可信带内，那么认为马尔可夫链是收敛的。R语言coda包的geweke.diag命令可实现，具体略。

三、贝叶斯统计的其他算法

贝叶斯后验分布因涉及多个参数，计算时面临多重积分的挑战，在早期限制了其在多元统计建模中的使用。MCMC方法的出现很好地解决了计算困难的问题，但仍存在一些问题，比如M-H算法接受率低下、高维参数空间Gibbs抽样效率低下及所有基于马尔可夫链抽样样本收敛性的问题。随着大型数据集和复杂模型的出现，出现了更先进的算法，如期望最大化（EM）法、积分嵌套拉普拉斯逼近（integrated nested Laplace approximation，INLA）、哈密尔顿蒙特卡罗（Hamiltonian Monte Carlo，HMC）、变分逼近（variational approximation，VA），这些算法有的可以避免收敛性问题，有的可以提高计算效率。这些先进的算法往往涉及复杂的数学及统计原理，本小节只简单介绍INLA及HMC的基本思想，读者感兴趣可参考相关著作。

INLA是综合贝叶斯推断、潜变量高斯模型、高斯马尔可夫随机场及拉普拉斯近似的方法。首先按条件概率重新编写后验分布，采用拉普拉斯近似值替代边际分布 p （ x θ ， y ），再计算分母边际分布的拉普拉斯近似值，采用高斯马尔可夫随机场的方法。其中拉普拉斯近似是基于泰勒级数展开式，将概率密度函数对数化后，按泰勒级数展开，分别计算二阶导数后，可以与正态分布关联，将求解的概率密度函数逼近特定均数及方差的正态分布。Finn Lindgren等人在2011年进行的研究提出了一种针对INLA的随机偏微分方程（SPDE）方法，该方法可以在较大的地理区域内进行更实际的空间建模，这些地理区域的曲面具有点数据，从而具有连续的协方差函数。因此，该方法在处理空间数据上具有极大的优势。R语言的inla包可以实现贝叶斯统计的非迭代抽样计算。

HMC方法又称为混合蒙特卡罗（Hybrid Monte Carlo），其采用系统动力学而不是概率分布来建立马尔可夫链的未来状态。哈密尔顿动力学是经典力学的重要理论，采用位置、动量及时间来刻画物体的状态。按能量守恒定律，物体总能量等于势能和动能之和，当物体在整个系统中运动时，势能和动能之间相互转化，可通过Hamiltonian方程的微分方程来描述。描述运动系统的当前状态需要有物体所处的位置式（position）和物体的动量。位置对应了统计中参数的值，而动量则是人为添加的一个变量，也被称为辅助变量。通过统计力学的正则分布与目标分布关联起来，通过欧拉法或蛙跳法进行模拟。HMC法的关键在于选取轨迹的长度及步长，当选择不当时会增加拒绝率，对于高维复杂模型可能耗时。No-U-Turn抽样（NUTS）是HMC的扩展，可解决上述问题。Stan是贝叶斯推断的C++库，基于No-U-Turn采样器（NUTS）进行抽样，进而对后验分布进行参数估计。不同统计软件有不同的Stan包，如rstan、pystan等。

四、贝叶斯统计常用软件简介

尽管贝叶斯统计涉及复杂的计算，但是有不同软件及程序包可以方便地进行计算。这些软件包括专用软件及通用软件。专用软件主要包括WinBUGS、OpenBUGS、JAGS及Stan等。通用软件包括R、SAS、Matlab、Python等。通用软件中可以自行编程进行计算，但多数是依靠相应的包，如R语言中的R2WinBUGS、rjags、mcmcglmm、mcmcpack等；SAS软件中的proc mcmc包；Stata 16.0软件中的bayes系列命令；Python软件中的pymcmc、pymc3等。读者只须根据自己的使用习惯进行选择。推荐读者采用R语言结合专用软件进行贝叶斯Meta分析。

BUGS是Bayesian Inference Using Gibbs Sampling的缩写，指基于Gibbs抽样的贝叶斯统计推断。BUGS软件最初于1989年由位于英国剑桥的生物统计学研究所（Biostatistics the Medical Research Council，Cambridge，United Kingdom）研制，现在由这个研究所和位于伦敦的S.t Mary’s皇家学院医学分院（the Imperial College School of Medicine）共同开发。BUGS所使用的编程语言类似R语言，通过设定贝叶斯模型的似然函数及先验分布，并将其自动转换成马尔可夫链的转移核，迭代后获得后验分布（目标分布）的抽样数据。BUGS主要包括WinBUGS和OpenBUGS，WinBUGS目前已暂停更新。OpenBUGS作为WinBUGS的开源版本，在持续更新中，目前的版本为OpenBUGS323，可于Windows、Linux等系统。尽管Open-BUGS在持续更新，但其基本代码、数据、初始值及操作与WinBUGS无明显差异。主要差异可参考：http：//www.openbugs.net/w/OpenVsWin。WinBUGS下载地址为https：//www.mrcbsu.cam.ac.uk/wp-content/uploads/2018/11/winbugs143_unrestricted.zip；OpenBUGS下载地址为http：//www.openbugs.net/w/OpenBUGS_3_2_3？action=AttachFile&do=get&target=Open-BUGS323setup.exe。

JAGS为Just Another Gibbs sampler的缩写，是由Martyn Plummer创作的基于Gibbs抽样的软件，可在不同的操作系统内使用。目前主要是与R语言相关的包联合应用，如R2jags、runjags、rjags等。JAGS语法命令与BUGS相似，多数代码是通用的；最新的版本是JAGS4.3.0，下载地址为https：//sourceforge.net/projects/mcmc-jags/。

Stan是用于贝叶斯统计建模和高性能统计计算的最新平台。可在社会学、生物和物理科学、工程和商业领域进行统计建模、数据分析和预测。用户使用Stan的概率编程语言指定对数密度函数，使用No-U-Turn抽样器进行贝叶斯统计推断、带变分推理（ADVI）的近似贝叶斯推理、带有优化的惩罚最大似然估计（L-BFGS）。更加方便地进行复杂的贝叶斯模型的计算。

第四节　贝叶斯Meta分析模型

一、贝叶斯Meta分析基本概念

循证医学是近些年来国际临床医学领域迅速发展的一种新的医学模式，其核心思想是任何医疗决策都应该基于客观的临床研究的证据。Meta分析是循证医学重要的技术和工具，是证据汇总的定量分析方法，它是将独立研究的结果结合起来，提供一个共同效应估计的方法，即合并不同研究的效应值。近些年来贝叶斯统计方法也被广泛应用于Meta分析，特别是复杂模型如纵向数据Meta分析、网络Meta分析等。贝叶斯Meta分析是采用贝叶斯原理对多个同类研究结果进行合并汇总的统计分析体系。随着计算方法的发展，其具有处理复杂模型以及解决传统方法无法解决的问题的能力。

贝叶斯Meta分析将估计参数视为随机变量，而非固定值。结合纳入研究的样本信息及先验信息，获得后验信息，并基于后验信息进行完整的参数估计及统计推断。因所有参数均有分布，因此所有参数均可获得点估计及区间估计。例如频率学派Meta分析中异质性方差 τ ² 是由矩估计获得，并无置信区间。

在经典Meta分析过程中经常会遇到以下情况：①罕见事件的Meta分析，即纳入研究数据呈明显稀疏性；②纳入研究数量较少。因此频率学派Meta分析会导致一系列有偏估计。贝叶斯统计框架下，通过引入弱信息先验分布（weakly informative priors，WIPs）来调整（或稳定）估计值。稀疏数据的Meta分析主要的问题在于效应量方差或标准差的估计，贝叶斯方法可针对方差项引入弱信息先验分布。基本思想是当采用无信息先验时，方差项后验抽样会出现过度离散现象，而采用弱信息先验分布可适当限制过度离散，进而获得稳健的后验估计。

二、经典Meta分析模型：正态-正态层次模型

经典Meta分析是典型的二步过程，其基本原理是：第一步，计算纳入Meta分析的每个研究的统计量。用相同方法来描述每个研究干预的观测效应量。第二步，通过对每个研究干预的观测效应量进行加权取平均数来获得总的合并干预效应。

从贝叶斯统计的角度来讲，Meta分析的本质为层次模型，根据各项研究的真实效应值是否相同可以分为固定效应模型和随机效应模型，如图1-8。虽然固定效应模型可以在不考虑异质性的情况下用于估计共同效应，但有可能夸大综合效应的精度。层次随机效应模型的贝叶斯方法可对异质性和未知的处理或研究效应进行估计。层次贝叶斯方法在处理Meta分析中出现的问题具有明显优势，例如在第二阶段先验分布进行选择，用于评估小型研究或非高斯效应的稳健推断，以及在假设的新研究中进行预测。本节将以正态-正态层次模型为例进行简单介绍。

固定效应模型下，假设所有纳入Meta分析的研究来自同一个总体，拥有相同的效应值（量级和方向均相同），因效应值固定，所以称为固定效应（fixed-effect，FE）模型。经典的Meta分析数据类型为连续型数据（如比值比对数、相对危险比对数、风险比对数、率差、均数差、标化均数差等效应量），假设有 n 个研究，每个研究的效应值为 y _i ，对应的标准误为 s _i ，固定效应模型为： y _i ～ Normal （ μ ， ）。式中， μ 为 y _i 的加权均值，权重由1/ 计算获得，此方法对应于频率学派的倒方差法。

随机效应（random-effects，RE）模型假设每个研究效应量来自各自的正态分布，而这些正态分布的均值服从更高水平的正态分布，随机效应模型即估计更高水平正态分布的均值，具体模型如下： y _i ～ Normal （ μ _i ， ）， μ _i ～ Normal （ θ ， τ ² ）。式中 θ 为合并的均值， τ ² 为研究间的方差，二者为模型待估计参数，在贝叶斯分析中，须指定先验分布。通常 θ 选用无信息正态分布， τ ² 选用无信息Gamma分布。

经典Meta分析的统计模型对于计算和解释Meta分析的结果非常重要，但由于固定效应模型和随机效应模型采用相似的公式计算统计量，有时可能得到相似的结果，以至于被误认为两个模型可以相互替换使用。但实际上，不同的模型基于不同的假设，并且提供不同的参数估计值。Meta分析中如何选择统计模型，历来存有争议。不同的人偏爱不同的事物，不同的统计学家和临床研究人员可能偏爱不同的统计模型，即使是第6版《考克蓝干预措施系统评价员手册》也未能提供权威统一的推荐意见。笔者认为，应从统计模型假说、Meta分析目的、纳入Meta分析的研究数量和样本量、研究间异质性、抽样框架等不同方面综合考虑来选择合适的统计模型，并建议在制定系统评价/Meta分析研究方案时就应该考虑选择合适的模型。基于随机效应模型的假说和抽样框架更符合实际、统计推断目的对研究者而言更有吸引力、从数学角度而言固定效应模型是随机效应模型的特例等方面来考虑，除了使用随机效应模型不可能（如只有一个研究）、不合理（异质性参数估计不可靠）等情况外，在Meta分析时应首先选用随机效应模型。但应牢记《考克蓝干预措施系统评价员手册》的重要观点：“决不应该根据异质性统计检验做出使用固定效应模型或随机效应模型的选择”。

［例1.10］ 以阿司匹林及安慰剂预防心梗死亡的数据为例，采用R语言rjags包演示贝叶斯Meta分析过程。数据如表1-3所示：

［解］ 对于研究MRC-1， OR =［（615-566）×557］ / ［566×（624-557）］=0.718， y ₁ = ln （ OR ）=-0.328 9， = Var （ ln （ OR ））=1 / （615-566）+1 / 566+（1 / 624-557）+（1 / 557）=0.038 9，其余类似。 OR 值为偏态分布，对数化后近似服从正态分布，每个研究的效应量及相应标准误结果如表1-3。

首先进行固定效应模型的计算，具体过程如下：

主要计算结果如下：

由上述结果可知， d =ln（ OR ）=-0.109， OR =0.891，95% CI （0.841，0.956）。结果中Rhat为潜在尺度缩减因子，均接近1，说明收敛良好。三个参数的核密度估计图及踪迹图如图1-9。

随机效应模型运行代码如下：

主要计算结果如下：

由上述结果可知， d =ln（ OR ）=-0.146， OR =0.868，95% CI （0.706，1.029），研究间异质性方差点估计及95% CI 为0.039（0.000，0.195）；结果中Rhat为潜在尺度缩减因子，均接近1，说明收敛良好。四个参数的核密度估计图及踪迹图如图1-10。

三、贝叶斯Meta分析模型拟合与比较

任何统计分析的一个重要部分都是评估来自特定模型的预测与观测数据的吻合程度。如果所选模型的预测与数据不太吻合，那么该模型的任何输出，如参数估计、预期净收益、结果的最佳处理策略均具有不确定性，都不能很好地反映证据基础。在贝叶斯Meta分析过程中会有不同的模型，如固定效应模型及随机效应模型等，如何合理的选择模型至关重要。

拟合优度的标准统计量，如残差、残差平方和（ RSS ）和残差偏差（通常称为似然比统计量）都可以在贝叶斯框架内计算，以评估模型的拟合。贝叶斯统计最常用的是残差偏差（residual deviance）。偏差统计量是用似然函数测量由特定模型作出的预测与观测数据的拟合度。定义为对数似然函数（Loglik _model ）的-2倍，即：

对应给定的模型，对数似然函数越大拟合越好，因此 D _model 就越小。偏差统计量越小，模型拟合越好。但是多小才算小呢？使用原始偏差统计量（ D _model ）的缺点是对这个问题没有明确的答案。残差偏差（ D _res ）通过提供一个参考点，帮助衡量模型拟合的好坏。定义为：

D _res = D _{mod el} - D _sat

其中 D _sat 为饱和模型的偏差。对于一个能很好地拟合数据的模型，期望个体对残差偏差度的贡献大致是服从自由度为1的卡方分布。如果观测数据正态分布的似然，这是精确服从卡方分布。如果对N个无约束数据点进行求和，那么残差偏差的期望近似服从自由度为N的卡方分布。

有效参数个数（effective number of parameters，pD）用来度量模型的复杂程度，定义为：

其中 为模型偏差的后验均值， D （ ）为参数为 时的偏差值。

Spiegelhalter等提出用偏差信息准则（Deviance Information Criteria，DIC）来比较贝叶斯模型。它利用和AIC类似的思想，使用DIC来综合表示模型的拟合度及模型的复杂程度，即：

DIC判断准则为如果两个模型的DIC相差大于5，说明DIC较小的模型更好；如果两个模型的DIC相差小于3，说明两个模型无太大差别。

以前文阿司匹林及安慰剂预防心梗死亡的数据为例，固定效应模型DIC=-2.7，随机效应模型DIC=-2.8，相差非常小。这种情况下，固定效应模型更为简单，且便于解释，因此本例更适合选用固定效应模型。

第五节　R语言函数小结

本章使用的R语言主要函数见表1-4。

参考文献

［1］GUREVITCH J，KORICHEVA J，NAKAGAWA S，et al.Meta-analysis and the science of research synthesis［J］.Nature，2018，555（7691）：175-182.

［2］SUTTON A J，ABRAMS K R.Bayesian methods in meta-analysis and evidence synthesis［J］.Stat Methods Med Res，2001，10（4）：277-303.

［3］茆诗松，汤银才.贝叶斯统计［M］.2版.北京：中国统计出版社，2012.

［4］SAMUEL KOTZ，吴喜之.现代贝叶斯统计学［M］.北京：中国统计出版社，2000.

［5］张天嵩，董圣杰，周支瑞.高级Meta分析方法：基于Stata实现［M］.上海：复旦大学出版社，2015.

［6］张天嵩，钟文昭，李博.实用循证医学方法学［M］.2版.长沙：中南大学出版社，2014.

［7］BROOKS S，GELMAN A，JONES G，et al.Handbook of Markov Chain Monte Carlo［M］.Boca Raton：Chapman and Hall/CRC Press，2011.

［8］GRILLI L，METELLI S，RAMPICHINI C.Bayesian estimation with integrated nested Laplace approximation for binary logit mixed models［J］.J Stat Comput Sim，2015，85（13-15）：2718-2726.

［9］GOMEZ-RUBIO V，RUE H.Markov chain Monte Carlo with the Integrated Nested Laplace Approximation［J］.Stat Comput，2018，28（5）：1033-1051.

［10］GRANT R L，CARPENTER B，FURR D C，et al.Introducing the StataStan interface for fast，complex Bayesian modeling using Stan［J］.Stata J，2017，17（2）：330-342.

［11］WARN D E，THOMPSON S G，SPIEGELHALTER D J.Bayesian random effects meta-analysis of trials with binary outcomes：methods for the absolute risk difference and relative risk scales［J］.Stat Med，2002，21（11）：1601-1623.

［12］NIKOLAKOPOULOU A，MAVRIDIS D，SALANTI G.Demystifying fixed and random effects meta-analysis［J］.Evid Based Ment Health，2014，17（2）：53-57.

［13］张天嵩.经典Meta分析统计模型的合理选择［J］.中国循证医学杂志，2020，20（12）：1477-1481. IrbbrXp0zDfK8cctJrDaYklt/ANFoYinlea5FFbc9eU3gFLB9t5YZcteVMgXyg0v