广义线性混合模型与拉普拉斯近似、高斯求积与零事件

写点证据合并方法里面稍微有点难度的内容，希望可以拓展拓展值得合作的优质科研团队。本来很早就有所打算，但是考虑到多数人的专业差异，担心写出来能看懂的不多，也就一直拿不定主意。考虑再三，意识到它的一些积极作用，还是写下来。首先能给可能在摸索的童鞋一点参考资料，这方面的内容国内书籍鲜有介绍；其次，可以让大家重新认识循证医学这门学科，可能跟大家平时对它的了解不太一样。当然，鉴于我们水平有限，如有描述不当的地方，还请斧正！非常感谢！

友情提示，后面的内容如果不慎看不懂，请个别大侠不要烦躁、也不要怀疑自己、更不要写一些很低级没品味的评论。又如果您看了之后因为内容晦涩难懂而大为不爽，敬请移步适合您的社区，那里有大量跟您臭味相投志同道合的”灌水侠“，如真如此，也表示感谢！咱们希望有更多高水平、高学术品位的同仁关注。下面正文开始。文章所有涉及的分析代码，均在原文附件，可以直接下载，也可联系通讯作者获取。

广义线性混合模型与拉普拉斯近似、惩罚拟似然、高斯-埃尔米特求积与零事件

零事件一直是计数资料的统计推断里难搞的问题。这个问题从1964年之前（我目前能追溯到的最早的时间）开始，距今已经半个世纪了。1964年这个时间应该合理，因为大家广泛使用的Odds ratio (OR)也才刚刚（1951年）被就职于美国国家癌症研究所的Cornfield 发明出来，而我们所要说明的零事件问题，又刚刚最先与OR扯上关系。怎么与OR先扯上关系，可以参考Fisher的确切概率检验和Logistic回归的历史。（注：零事件与相对效应risk difference也有关联，但是早已被解决，此文不予讨论）

大量统计学家毕其精力力求解决零事件数据的统计推断问题。最开始（可能）由一位天才统计学家Plackett于1964年提出解决方案，通过对OR的对数使用泰勒级数进行展开，寻找其二阶近似下的最优解。结果发现，对四格表中每项元素添加0.5个虚拟样本时，能获得最好的二阶近似估计。这也是当前最被广泛使用的连续性校正法的来历。当然，之后又有众多统计学家对加0.5校正的方法进行了各种研究和改良，包括经验校正法，样本倒数校正法之类的，在此不多赘述。连续性校正法”一统江湖”数十载，直到1985年另一位天才统计学家，Peto，提出了另一种充满竞争力的方法，Peto OR。该方法是利用试验组实际观察到的发生例数与其期望发生例数的差值而构建的指数模型，用以估算OR值。实际上这个方法跟OR的定义扯不上任何关系，但是因为结果神接近OR，大家也习惯把Peto法获得的统计量称为OR。但是在2014年Brockhaus明确指出，Peto OR不是OR，应该当作一种新的统计量。Peto OR在处理四格表中出现单个零事件的统计性能明显优于连续性校正。以至于后面Cochrane官方指南明确指出不推荐使用连续性校正作为首选方法。

单个零事件的问题似乎解决了，但是出现两个零怎么办？也就是试验组和对照组都没有观察到事件的发生，后文我们称为双臂零事件。我们前面系列论文已经提到，双臂零事件的研究通常都被忽略不计了，一些统计学家认为它没有信息。从似然理论上讲这个观点是有道理的，就单个研究来看，零事件的似然值确实是零，无论是边际似然，还是条件似然。以至于这个观点后面极大的影响了Meta分析中对双臂零事件的处理：无视它的存在。我们已经通过多种角度论证了对于Meta分析来说，更大一点，从循证医学角度来看，双臂零事件是存在信息的。这里不再展开。那么，如何去处理双臂零事件呢？

对于单个研究来讲，似乎无解。部分研究人员天真的继续把另一个零也加上0.5，用来处理双臂零事件。这是非常错误的方法。单个零的时候，添加0.5已经会导致巨大的偏差，出现两个零继续加0.5，那还得了。就问结果你敢不敢信？但是对于多个研究一起，那是可以解决的。一群人里面有几个没钱，大家还可以给借点，但是如果只有一个人向谁借？所以双臂零事件的解决方案类似于”借钱“，从没有出现零事件的研究里面”借“。但这个”借“不是真的借，只是虚晃一枪，不带走一片云彩。也就是建立广义线性混合模型。

广义线性混合模型（Generalized Linear MixedModel，GLMM）最早是由Platt WR于1999年应用于Meta分析的数据合并。它属于广义线性模型的拓展，加入了随机效应。广义线性模型指的是一类能写成指数型分布族形式统计结构的回归模型，包括众所周知的普通线性模型、logistic回归模型、泊松回归模型等七类。在不考虑研究间异质性的前提下，logistic回归模型可用于粗略的Meta分析数据合并。但这种将不同研究直接看作一篇整合的大研究而直接采用logistic回归合并的方法并不可取，原因在于不同研究本身存在不同，包括不同的人群、不同的实施人员、不同的质控过程、不同的研究目的等。广义线性混合模型则能很好的处理研究间的差别以及潜在聚集性，它将研究内个体当作1水平，研究当作2水平，不仅考虑个体间的差异，也考虑研究间的差异。本质上是在多水平的框架内加入广义线性回归，因此也称为多水平广义线性模型，或分层广义线性模型，都行。处理异质性不是它最吸引人的性质。最吸引人的莫过于，它连双臂零事件研究也能处理。如何如理？

GLMM的参数估计需要迭代过程，迭代的初始值是不加入随机效应的普通广义线性模型计算的结果。如上所属，普通广义线性模型直接把所有纳入研究当作一个整体，样本量归公，发生例数归公，然后就只有一个吸收了所有研究的大研究。然后直接Logistic回归，获得初始参数（logOR等），这样零事件研究自然而然被一起归公纳入计算了，大树底下好乘凉。显然这个初始值很不准确，甚至都没有加权，有很大可能导致Simpson悖论。所以需要进一步剥离出随机效应了。考虑随机效应后，再次建立条件似然函数，使用上一步获得的初始值，请来牛顿-拉普森帮忙进行一通迭代，分离出随机效应，同时获得单篇研究的效应以及总效应，哦，零事件研究的效应也被分出来了。所以没什么是爵士搞不定的。

此时不得不提另一个方法，MH法。是统计学家Mantel N和Haenszel W的合体。目前大家了解到的用在Meta分析里面的MH法，不是经典的MH法。经典MH法，属于one-stage法，思路也类似于上述广义线性模型的”归公“，只是在归公的时候，对单个研究试验组的Odds和对照组的Odds进行了样本加权。达到合并不同分层的亚组结果，顺带把有零事件的亚组也解决了（只能解决有单个零的情况）。很显然，MH法也只能在把所有研究当作一个整体的时候才能发挥解决单臂零事件的作用。如果对单个研究进行一对一”辅导“，无法发挥应有性质。只能借助连续性校正法去处理零事件。所以，本文前面回顾处理单臂零事件的方法时，不提此方法。

言归正传。对Meta分析数据建立GLMM模型非常容易，无非就是在Logistic模型里的各个变量加个下标j，表示层次，个体水平-研究水平。但是参数的估算是及其复杂。大佬们通过发表的文献告诉我，GLMM建立的似然函数没有闭环解，简单理解为通过建立方程组根本无解。这个时候需要使用函数逼近或者数值优化技术了。我们可以想象，似然函数也是函数，因此可以建立一个函数去无限逼近这个似然函数的。另一种很好的思路是，对似然函数作无限切割，请高斯-埃尔米特过来帮忙求个积。逼近的思路，有两个常用的方法，惩罚拟似然函数，通过”阉割“掉一部分余项等复杂的处理建立拟似然函数，然后喊来泰勒把这个拟似然函数进行强行展开，再喊来拉普拉斯帮忙做个二阶逼近。第二种方法，拉普拉斯近似，拉普拉斯自己的主场。建立什么拟似然函数，忒复杂了。直接喊上泰勒把似然函数展开，然后自己给来个二阶逼近不就完事儿？求积的思路，主要介绍一种，也就是最常用的，高斯-埃尔米特求积。通过对似然函数建立n个n阶非线性方程组对似然函数进行极小极大逼近，没错，就是极小极大逼近。

我们已经知道，GLMM框架是完全能够处理包含零事件研究的Meta分析数据的。但它需要借助以上三种参数估算方法才行。这三种参数估算方法在遇到如雷贯耳的包含零事件研究的Meta分析数据时，表现怎么样？会不会怯场？谁知道呢。本文于是建立了5种GLMM模型，探讨这三种参数估计方法分别在这种GLMM模型下，处理包含零事件研究的Meta分析数据的统计性能，进行了全面考察。这五种GLMM模型分别是：模型1-研究效应随机模型，模型2-研究效应及基线风险双随机模型，模型3-重参数化研究效应随机模型，模型4-重参数研究效应及基线风险双随机模型，模型5-双变量研究效应随机模型。这里我们使用研究效应和基线风险，替代斜率和截距。别问为什么，Higgins他们都这么描述。或许这样能更好的与原始研究进行区分？例如：模型1也可以称为随机斜率模型（随机斜率模型是原始研究里面处理各个水平间的协方差的常用叫法）。模型3和模型4的重参数化，是对协方差结构的”类中心化”处理，有利于更好的反映研究间的异质性。

方法

数值模拟

本研究采用数值模拟（Numerical simulation）技术，对不同GLMM模型及参数估计方法的组合方案的统计性能进行评估及比较。具体思路如下：首先，根据实际应用情境设置不同的参数及真实值，通过模拟技术产生不同的情境下的数据；其次，将产生的数据分别采用上述15种组合方法进行合并，以获得模拟的Meta分析结果；最终将不同方法获得的结果与真实值进行比较，用以获取不同方法在不同情境下的统计性能表现。模拟参数设置部分较为简单，不进行详细展开。能看懂到这个地方的同仁，秒懂，也不需要展开。直接给出以下表格展示模拟参数的设置情况。本研究针对每种情景（5*5=25）均模拟20000次，也就是产生20000个Meta分析的数据。

统计性能指标

经典概率论与数理统计理论常将点估计的准确性和区间估计的精确性作为待评估方法的统计性能评价指标。其中点估计（即效应量）的准确性采用偏差百分比（Percentage bias，PB）来衡量，即Meta分析的观察值与真值的差值与真值的比例，假设观察值用θ表示，真值用表示，则偏差百分比PB的计算公式为。理论上，当观察值等于真实值时，偏差百分比为0；当观察值小于真实值时偏差百分比为负数，表示效应量存在低估；当观察值大于真实值时，偏差百分比为正数，表示效应量存在高估。区间估计采用均方标准误（Mean square error, MSE）来衡量，反映观察值方差的精确性，不妨以Var表示方差，均方标准误MSE的计算公式为。当MSE越小，说明精确度越高，区间将越窄，结果准确性的把握度越大；反之，MSE越大，说明精确度越低，区间将越宽，结果准确性的把握度越小。值得注意的是，并非MSE越小越好，因为当区间越窄，区间内包含真值的可能性会降低，而区间越宽，则其包含真值的可能性越大。因此，我们引入第三个指标，即区间覆盖率（Coverage possibility，简称Coverage），反映的是在模拟中产生的2万个Meta分析合并OR的区间有多少次覆盖了OR真值。由于待评估的5种方法包括经典频率学方法和贝叶斯方法，此处区间既表示置信区间（频率学派），也表示可信区间（贝叶斯学派），且取95%置信或可信区间作为置信域或可信域。理论上，区间覆盖率应保持在95%左右，覆盖率过低表明结果的准确度和精确度不够，过高则会增加I类错误的风险。

同时，考虑到算法优化中收敛性的问题，本研究还加入了收敛成功率这一指标，最大迭代次数设置为10次。这意味着，在迭代10次及以内能达到收敛的，我们认为收敛成功，反之则不成功。此处设置10次作为阈值是考虑到计算工作量的问题。一般的软件程序，会迭代1.5万次，这种迭代次数下一个Meta分析大约需要5分钟，每个Meta分析用15种方法就是5*15=75分钟，20000个Meta分析就是75*20000=1500000分钟=25000小时=1042天=2.85年。25个情景，每个情景2.85年，一共需要2.85*25=71.35年。这也是为什么我们设置迭代次数为10次的原因。其次，通常情况下，具有收敛性质的过程一般7-8次就能收敛。因此我们不考虑更高次数的收敛。

亚组分析

在大约上世纪90年代的时候，统计学家们发现使用似然估计的方法当发生例数较小的时候，会出现一种偏差，叫做小样本偏差。随后他们建议，当一项研究的发生例数少于10个的时候，不建议做回归分析。同时，当发生例数大于10个的时候，每10个例数最多只能对应一个自变量。也就是说，总例数是50的时候，模型最多只能放入5个自变量。这也是大家现在广泛使用的EPV(events per variable)原则，也称之为rule of 10。

本研究的15种组合方案，均基于似然估计，因此也需要考虑EPV对结果的影响。鉴于此，本研究将每个情景下模拟的每个Meta分析的试验和对照组总样本按照如下规则进行分组，共分为6组：1）试验组和对照组总样本量均≥10；2）两组中，一组总样本≥10，另一组小于10但同时≥5；3）两组均小于10，但两组均≥5；4）一组大于≥10，另一组小于5；5）一组≥5，另一组小于5；6）两组均小于5。

结果

收敛性

下图示15种组合方案在25种情景下迭代过程中的收敛成功率。很明显我们能看到，惩罚拟似然法在模型3、模型4、模型5框架下，收敛成功率非常低。因此，我们首先将这两种方案排除，剩余12种。当然，值得注意的是，这个收敛成功率理论上是低估的。正如我们之前所述，本研究以10次迭代为标准计算收敛成功率。部分情况可能迭代11次，或100次能收敛，但是在本研究视为不收敛。如果扩大这个限制，成功率理论上更高。

极大值比例

下图示OR=1以及Tau=0.2这一情景下15种方法在6个样本组别下偏差的分布情况（越偏离0，偏差越大）。此处我们只展示第一个情景原因是当OR和Tau持续变大的时候，可视化效果越来越差，也就是越往后，大部分小单元的分布曲线完全无法分辨。结果提示，在高斯-埃尔米特求积的方法下，5种广义线性混合模型均表现出较拉普拉斯近似、惩罚拟似然更大的偏差。我们统计了15种方法在6个样本组别、25个情景下出现极大值的比例(OR > 250)。结果发现，高斯-埃尔米特求积出现极大值的比例最高。注：极大值是由于小样本偏差导致。此结果符合预期，因为惩罚拟似然法的初衷之一，就是通过惩罚某些项来解决小样本偏差。

总事件数对GLMM估算的影响

我们比较了不同样本组别下15种方案估算的logOR的偏差，结果发现，当试验组和对照组两组样本量均大于或等于10的时候，GLMM的偏差中位数很小，几乎可以达到无偏估计，同时偏差百分比超过50%的比例也非常小。一组样本量大于等于10，一组大于5的情况下，GLMM的表现在一些情景下也较好，但是部分情景下表现不佳。当两组均小于10，或者一组小于5的剩下四种境况下，结果偏差较大，更高比例的Meta分析偏差百分比超过50%（如下图）。这一结果验证了EPV原则，同时也表明，EPV原则不仅是经典回归模型的原则，也适用于GLMM，也适用于采用GLMM模型的Meta分析。

剩余12种方案的比较

在淘汰三种收敛成功率较低的方案后，我们对剩余12种方案的偏差、均方标准误、覆盖率进行了全面比较，见图4。结果提示，综合偏差、均方标准误、区间覆盖率后，高斯-埃尔米特求积、惩罚拟似然下的GLMM模型搞不赢拉普拉斯近似下的GLMM。进一步比较拉普拉斯近似下的5种GLMM，我们发现，模型1和模型3具有更小的偏差以及均方标准误，但模型2、4、5具有更好的区间覆盖率。这意味着将基线风险视为随机项，可以获得更好的真实值所在范围的“预测”，但是损失了一定的准确性。同时我们发现，随着异质性增大，所有模型的统计性能都在“恶化”，并且当异质性非常明显的时候，所有模型的表现均较差。

讨论

主要结论

本研究获得了三个重要的结论：1. 广义线性混合模型是可以用来处理双臂零事件研究的，并且当纳入研究总事件数足够的情况下，能获得几乎无偏的估计；2.惩罚拟似然法以及高斯-埃米尔特求积尽管在有些情况下能体现出他们各自的优越性，但是在处理包含零事件的罕见事件Meta分析数据上性能并不优于拉普拉斯近似。3. 考虑到罕见事件Meta分析的研究间异质性通常较小，研究效应随机模型（随机斜率模型）在处理包含零事件的罕见事件Meta分析数据上可能性能优于研究效应及基线风险双随机模型（随机斜率及随机截距模型）。其余讨论内容请看原文。