Meta分析中随机效应模型的Gibbs抽样及其应用，本文主要内容关键词为：及其应用论文,效应论文,模型论文,Meta论文,Gibbs论文，此文献不代表本站观点，内容供学术参考，文章仅供参考阅读下载。
Meta分析效应合并时应选用固定效应模型还是随机效应模型，一直存有争议。尽管不少学者认为用随机效应模型得到的meta结果更“保守”，但定量刻划随机效应的方法很少，而且基于统计量Q的异质性检验的功效较低[1]。近年来，meta分析的贝叶斯方法逐渐受到重视。但贝叶斯方法计算参数的后验分布常面临多重积分问题，除了少数特定假设下可以得到解析表达式的情形。Gibbs抽样使贝叶斯方法中许多看起来复杂困难的计算变得简单直观。应用Gibbs抽样可以很方便地解决经典meta分析方法不能解决的一些问题。本文探讨分层贝叶斯框架下，随机效应模型meta分析的Gibbs抽样及其应用。
一、原理与方法
Gibbs抽样的基本原理是通过构造马尔可夫链模拟参数的后验分布，关键是得到所有待估参数的完全条件分布，即每一个待估参数在假设其他待估参数已知时的边缘后验分布[2]。
（一）基于分层贝叶斯的随机效应模型meta分析的基本形式
在分层贝叶斯框架下，随机效应模型meta分析可表示为
这里，Y[,i]为纳入meta分析的第i个研究的效应值，μ[,i]是第i个研究的“真正效应”(true effect)，μ[,1]，μ[,2]，…，μ[,k]相互独立。μ即为我们所关心的效应合并值（效应的平均水平或总体水平），τ[2]为研究间变异，即随机效应。为计算方便，可假定σ[2,i]已知，考虑一般情形，超参数μ、
Gibbs抽样时，先对参数μ[,i]、μ、τ[2]赋初值，再依次从式(1)、(2)、(3)各参数的完全条件分布中产生Gibbs抽样值。得到t次抽样值后，前m次用于“退火”(burn in)，后t-m次的抽样值用于计算各参数的边缘后验密度。
在模型构造中，我们给出的是最常见的情形。实际应用时，先验分布的选取、似然函数的构造可以方便、灵活地进行调整。如果我们对待估参数没有任何特别的先验信息，μ和τ[2]可以用无信息先验，这时，参数的完全条件分布的形式将更简单。
（三）WinBUGS简介
WinBUGS(Bayesian Inference Using Gibbs Sampling)是英国剑桥公共卫生研究所的MRC Biostatistics Unit推出的用马尔可夫链－蒙特卡罗(Markov chain Monte Carlo,MC-MC)方法进行贝叶斯推断的专用软件包[3]，可免费下载(http://www.mrcbsu.cam.ac.uk/bugs/)。用WinBUGS可以很方便地对许多常用的模型和分布进行Gibbs抽样。WinBUGS用有向图模型方式(directed graphical model)对模型进行直观的描述，并给出参数的Gibbs抽样动态图、用Smoothing方法得到的后验分布的核密度估计图、抽样值的自相关图及均数和置信区间的变化图等，使抽样结果更直观、可靠。Gibbs抽样收敛后，可以得到参数的后验分布的均数、标准差、95%置信区间和中位数等信息。
二、应用实例
（一）资料来源
表1的第(1)、(2)、(6)、(10)、(14)栏给出了1993～1998年间发表的有关中国人HBV、HCV及其双重感染与原发性肝癌关系的16个病例－对照研究的基本数据（未感染用上标[00]表示，HBV感染用上标[10]表示，HCV感染用上标[01]表示，双重感染用上标[11]表示）[4～16]。所有资料经检索国内外主要生物医学文献数据库，并严格按照meta分析前制订的文献纳入和剔除标准得到。OR的自然对数Y[10,i]、Y[01,i]、Y[11,i]为纳入meta分析的第i个病例－对照研究HBV、HCV及其双重感染与肝癌关系的效应值，Var(Y[10,i])、Var(Y[01,i])、Var(Y[11,i])分别为Y[10,i]、Y[01,i]、Y[11,i]的方差。基本计算结果见表1。
（二）建立模型及Gibbs抽样
由表1的原始数据用二项分布建模，通过拟合三个logistic模型，在WinBUGS软件包上通过Gibbs抽样得到HBV、HCV及其重叠感染与肝细胞癌的关系的meta分析效应合并值、研究间方差等参数的后验分布。以HBV感染为例，假设第i个研究对照组和病例组HBV感染的概率分别是p[00,i]、p[10,i]，在n[00,i]个对照和n[10,i]个病例中，HBV感染人数r[00,i]和r[10,i]服从二项分布。令
抽样过程：共进行55000次迭代，前5000次用于“退火”，消除初值的影响。后50000次迭代时每隔5000交记录一次抽样结果，每个参数有10个遍历均值。经比较，这些遍历均数的波动范围小于10[-2]，表明抽样结果基本趋向稳定，由此判断Gibbs抽样收敛。用5000次“退火”后的20000次抽样结果作为参数后验分布的估计值。
（三）结果与分析
主要参数（效应合并值μ、研究间方差τ[2]、OR值）的Gibbs抽样结果见表2。
WinBUGS可以给出参数的很多信息，特别是τ[2]作为刻划随机效应大小的指标，经典方法只能给出τ[2]的矩估计，而Gibbs抽样能得到τ[2]的后验均数及95%CI。此外，WinBUGS还可以给出每个纳入研究的后验均数β[10,i]、β[01,i]、β[11,i]，进一步得到每个纳入研究HBV、HCV及其重叠感染与肝癌关系的OR值的后验分布，限于篇幅，我们未报告这些结果。
为作对比，我们把表1中原始数据为0的格子一律用0.5代替，用经典的随机效应模型DerSimonian-Laird法进行meta分析。计算结果见表3。
表2 HBV、HCV及其重叠感染与肝癌关系的Gibbs抽样结果
表3 随机效应模型DerSimonian-Laird法的meta分析结果
从结果的比较看，关于HBV、HCV感染，两种方法得到参数μ[10]、μ[01]的点估计和区间估计基本一致，经典方法估计的效应合并值总的来说略小于GibbS抽样得到的效应合并值。对双重感染，因为病例－对照研究中对照组HBV、HCV双重感染率很低，四格表资料中一些格子的实际频数很小，甚至为0（本例纳入的16项研究中有4项研究的原始数据包含0，见表1）。经典方法或者舍弃这部分包含0的数据，或者用0.5代替0进行近似计算。本例用0.5代替0作近似计算得到的效应合并值为4.026，Gibbs抽样的效应合并值为4.489，差异相当大。因此，当纳入研究中数据极端值较多时，经典方法得到的效应合并值是不可靠的。
对随机效应的估计，Gibbs抽样和经典方法得到的HBV、HCV感染的(τ[10])[2]、(τ[01])[2]也基本一致，而且Gibbs抽样还可得到研究间方差的95%CI。但经典方法估计双重感染的研究间方差(τ[11,2])[2]为0，而Gibbs抽样估计(τ[11,2])[2]为0.237，这说明数据中极端值较多、而随机效应又相对较小时，经典方法几乎不能识别随机效应。
三、讨论
1.Meta分析中对模型选择问题的争议焦点在于效应合并时是否应考虑随机效应，如何定量刻划。尽管经典的随机效应模型meta分析基于Q检验可得到研究间方差的矩估计，但很难计算其95%CI，且检验功效较低。而且，经典方法将效应合并值、研究间方差等参数视为未知的某个数，忽略了参数估计的不确定性。基于分层贝叶斯的随机效应模型将未知参数视为服从某一分布的随机变量，由先验和样本信息，得到参数的后验分布。
2.经典方法很难对复杂模型进行meta分析。对小样本资料，如果不符合正态近似条件，经典方法也无法处理。而且，当纳入研究中数据的极端值较多时，经典方法很难识别随机效应。本文的实例应用结果表明，用Gibbs抽样不但能灵活地构造模型对复杂问题进行meta分析，而且当纳入meta分析的原始数据极端值较多，经典方法不再适用时，用Gibbs抽样仍然可以很方便地得到效应合并值、研究间方差等参数的后验均数和95%CI。
3.Gibbs抽样建模时要求由先验和似然函数，给出每个参数的完全条件分布。当参数的完全条件分布不能表示为显式形式时，不能直接进行Gibbs抽样。这时，可以用Metropolis-Hastings进行抽样。Metropolis-Hastings是比Gibbs抽样更一般化的MCMC方法，其应用范围也更广。
4.应用WinBUGS软件包，可以大大提高Gibbs抽样的速度和效率，程序的编写也相当简单。本文的实例应用中，我们直接利用原始数据用二项分布建模，拟合logistic模型，同时，还可以在模型中加入其他参数，如特异性协变量，探讨其对效应合并的影响。在先验分布的选取上，可以假设效应合并值的先验分布是t分布，而不正态分布。