Multilevel Meta-Analysis
Published:
10 “多层次”Meta 分析
来源页面：https://doing-meta.guide/multilevel-ma
本页内容
- 10 “多层次”Meta 分析
- 10.1 Meta 分析的多层次本质
- 10.2 在 R 中拟合三层 Meta 分析模型
- 10.3 三层模型中的亚组分析
- 10.4 稳健方差估计
- 10.5 聚类 Wild Bootstrap
- 10.6 问答
- 10.7 小结
欢迎来到高级方法部分。在本指南前一部分中，我们已经深入讨论了那些几乎对每一项 Meta 分析都非常重要的主题。有了这些基础，我们现在可以继续学习一些更高级的方法。
我们之所以把下面这些方法称为“高级”，是因为它们的数学基础更复杂，或者因为它们在 R 中的实现更讲究一些。不过，只要你已经认真学完前面几章，就完全有能力理解并应用接下来要介绍的内容。下面许多主题其实都足以单独写成一本书，因此我们在这里给出的只能算是简要导论；在合适之处，我们也会提供进一步阅读的文献线索。
高级方法部分的第一章讨论的是“多层次”Meta 分析（multilevel meta-analysis）。你也许会好奇，为什么我们要给 “multilevel” 加上引号。把某项研究称为“多层次”Meta 分析，似乎暗示着它相对于“标准”Meta 分析而言是一种特别或非同寻常的方法。
然而事实并非如此。为了合并研究结果，任何 Meta 分析模型都预设我们的数据具有多层次结构。在前面几章中，我们其实已经在不知不觉间多次拟合过多层次的 Meta 分析模型。
当人们谈到 multilevel meta-analysis 时，他们通常真正指的是 三层 Meta 分析模型（three-level meta-analysis models）。这类模型确实与我们已经熟悉的固定效应模型和随机效应模型略有不同。因此，本章首先会说明：为什么 Meta 分析天然就蕴含多层次结构，以及我们如何把传统的 Meta 分析扩展为三层模型。和以往一样，我们也会通过一个实操示例看看如何在 R 中拟合这类模型。
10.1 Meta 分析的多层次本质
要理解为什么 Meta 分析天然具有多个层级，我们先回到第 4 章讨论过的随机效应模型公式：
\[\hat\theta_k = \mu + \epsilon_k + \zeta_k \tag{10.1}\]我们已经讨论过，在随机效应模型中之所以要引入 \(\epsilon_k\) 和 \(\zeta_k\)，是因为我们假定数据存在两类变异来源。第一类是单个研究的抽样误差（\(\epsilon_k\)），它会使效应量估计偏离真实效应量 \(\theta_k\)。
第二类是 \(\zeta_k\)，它代表研究间异质性。之所以会有这种异质性，是因为某项研究 \(k\) 的真实效应量本身也只是从一个更高层次的 真实效应量分布 中抽取出来的。因此，在随机效应模型中，我们的目标是估计这个真实效应量分布的均值，也就是 \(\mu\)。
两个误差项 \(\epsilon_k\) 与 \(\zeta_k\)，恰好对应了 Meta 分析数据中的两个层级：一是“参与者”层级（第 1 层），二是“研究”层级（第 2 层）。下图形象地展示了这一结构。
图 10.1：传统随机效应模型的多层次结构。
在最低层（第 1 层），我们有参与者（或者病人、样本、标本等，具体取决于研究领域）。这些参与者又隶属于更大的单位，也就是被纳入 Meta 分析的各项研究；这一更高层次就构成了第 2 层。
当我们开展 Meta 分析时，第 1 层的数据通常已经以“汇总”的形式到达我们手中。例如，论文作者给出的往往是样本均值和标准差，而不是原始个体数据。但第 2 层，也就是研究层面的合并，则需要由 Meta 分析本身来完成。传统上，这样的数据结构被称为 嵌套（nested）：也就是说，参与者“嵌套”在研究之中。
再回到式（10.1）。实际上，这个公式已经隐含描述了 Meta 分析数据的多层次结构。为了让这一点更加直观，我们可以把它拆成两个公式，每个公式对应一个层级：
第 1 层（参与者）模型：
\[\hat\theta_k = \theta_k + \epsilon_k \tag{10.2}\]第 2 层（研究）模型：
\[\theta_k = \mu + \zeta_k \tag{10.3}\]你可能已经注意到，我们可以把第二个方程中对 \(\theta_k\) 的定义代入第一个方程。这样得到的结果，就与前面的随机效应模型公式完全相同。固定效应模型也可以写成同样的形式，只不过把 \(\zeta_k\) 设为 0 即可。由此可见，传统的 Meta 分析模型本身就“内建”了多层次属性；出现这种属性，是因为我们假定数据中的参与者嵌套于研究之内。
这就说明，Meta 分析天然拥有多层次结构。为了更好地刻画生成数据的某些机制，我们还可以把这种结构进一步扩展。这时就轮到 三层模型 出场了。
在 Meta 分析中合并效应量时，统计独立性是一项核心假设。如果效应量之间存在依赖关系，也就是效应量彼此相关，那么这会人为降低异质性，并可能导致假阳性结果。这类问题被称为 分析单位错误（unit-of-analysis error）。效应量之间的依赖可能来自不同来源：
来自原始研究作者引入的依赖。 例如，研究者可能在多个研究点收集数据，把多个干预组与同一个对照组合并比较，或者使用不同问卷测量同一结局。在这些情况下，我们都可以认为报告出的数据内部存在某种依赖性。
来自 Meta 分析者自身引入的依赖。 例如，你正在做一项关于某种心理机制的 Meta 分析，这项分析纳入了来自世界不同文化区域的研究（如东亚文化与西欧文化）。根据所研究机制的不同，来自同一文化区域的研究结果，可能会比来自不同文化区域的研究结果更为相似。
我们可以通过在 Meta 分析模型结构中加入第三层，来处理这类依赖。例如，可以建立一个模型，把基于不同问卷得到的效应量视为嵌套在研究之内；也可以建立一个模型，把研究视为嵌套在不同文化区域之内。这样就形成了一个三层 Meta 分析模型，如下图所示。
我们看到，三层模型包含三个汇总步骤。首先，原始研究者先把单个参与者的结果进行“汇总”，并报告聚合后的效应量。随后，在第 2 层，这些效应量嵌套在若干 簇（cluster）之中，记作 \(\kappa\)。这些簇既可以是单个研究（也就是一项研究内嵌套多个效应量），也可以是研究的某种分组（也就是一个分组内嵌套多项研究，而每项研究只贡献一个效应量）。
最后，再把各个簇的聚合效应汇总起来，得到总体真实效应 \(\mu\)。从概念上看，这个平均效应与固定效应模型或随机效应模型中的合并真实效应 \(\mu\) 非常接近。不同之处在于，三层模型明确考虑了数据中存在相互依赖的效应量。
我们也可以沿用前面的层级记号，把三层模型写成公式。与之前最大的不同在于，这次我们需要定义三个方程，而不是两个：
第 1 层模型：
\[\hat\theta_{ij} = \theta_{ij} + \epsilon_{ij} \tag{10.4}\]第 2 层模型：
\[\theta_{ij} = \kappa_j + \zeta_{(2)ij} \tag{10.5}\]第 3 层模型：
\[\kappa_j = \mu + \zeta_{(3)j} \tag{10.6}\]其中，\(\hat\theta_{ij}\) 是真实效应量 \(\theta_{ij}\) 的估计。这里的下标 \(ij\) 可以读作“嵌套在簇 \(j\) 中的某个效应量 \(i\)”。参数 \(\kappa_j\) 表示簇 \(j\) 内的平均效应量，而 \(\mu\) 则表示总体平均效应。和前面一样，我们可以把这些公式拼接起来，最终把模型写成一行：
\[\hat\theta_{ij} = \mu + \zeta_{(2)ij} + \zeta_{(3)j} + \epsilon_{ij} \tag{10.7}\]与随机效应模型相比，这个公式现在包含了 两个 异质性项。其一是 \(\zeta_{(2)ij}\)，它代表第 2 层的 簇内 异质性，也就是说，簇 \(j\) 内各个 真实 效应量围绕均值 \(\kappa_j\) 呈分布。其二是 \(\zeta_{(3)j}\)，也就是第 3 层的 簇间 异质性。因此，拟合三层 Meta 分析模型不再只是估计一个异质性方差参数 \(\tau^2\)，而是要估计两个 \(\tau^2\)：一个对应第 2 层，另一个对应第 3 层。
{metafor} 包非常适合用于拟合 Meta 分析中的三层模型，它采用（限制）最大似然方法来完成估计。此前我们主要使用 {meta} 包的函数开展 Meta 分析，是因为这个包在技术细节上更简洁，也更适合初学者。不过，正如我们在前面章节已经看到的那样，只要数据整理得当，{metafor} 其实也并不难用。下一节我们就会具体看看，如何借助 {metafor} 在 R 中拟合三层模型。
补充说明：较新的 {meta} 包版本也支持三层 Meta 分析模型。在第 4 章介绍过的所有 Meta 分析合并函数中，现在都新增了一个 cluster 参数，用来定义数据集中每个效应量所属的第 3 层簇变量。如果指定了 cluster，函数就会自动拟合一个分层三层模型。例如，我们可以把前面某章中的 Meta 分析直接改写为三层模型：metagen(TE, seTE, cluster = InterventionType, data = ThirdWave)。不过，学习如何用 {metafor} 拟合三层模型依然很有价值：第一，{meta} 在后台也是借助 {metafor} 来拟合这类模型；第二，本章将要介绍的 rma.mv 函数用途非常广，不仅可以处理“简单”的分层三层模型，在后面的稳健方差估计部分我们还会继续看到它的扩展用途。
10.2 在 R 中拟合三层 Meta 分析模型
正如前面所提到的，我们需要 {metafor} 包来拟合三层 Meta 分析模型。因此，第一步就是先从库中加载它。
library(metafor)
在下面的实操示例中，我们将使用 Chernobyl 数据集。这个数据集大体上基于一项真实的 Meta 分析：该研究考察了 1986 年切尔诺贝利反应堆灾难所造成的电离辐射（“核沉降”）与人类突变率之间的相关关系。
“Chernobyl” 数据集
Chernobyl数据集是 {dmetar} 包的一部分。如果你已经安装并加载了 {dmetar}，运行data(Chernobyl)就会自动把该数据集保存到当前 R 环境中，之后即可直接使用。如果你尚未安装 {dmetar}，也可以从网络下载该数据集对应的
.rda文件，将其保存到工作目录中，然后在 RStudio 窗口中点击导入。
# 从 'dmetar' 加载数据集
library(dmetar)
data("Chernobyl")
为了查看数据的基本结构，我们可以使用 head 函数。它会打印出刚刚加载到全局环境中的数据框前六行。
head(Chernobyl)
## author cor n z se.z var.z radiation es.id
## 1 Aghajanyan & Suskov (2009) 0.20 91 0.20 0.10 0.01 low id_1
## 2 Aghajanyan & Suskov (2009) 0.26 91 0.27 0.10 0.01 low id_2
## 3 Aghajanyan & Suskov (2009) 0.20 92 0.20 0.10 0.01 low id_3
## 4 Aghajanyan & Suskov (2009) 0.26 92 0.27 0.10 0.01 low id_4
## 5 Alexanin et al. (2010) 0.93 559 1.67 0.04 0.00 low id_5
## 6 Alexanin et al. (2010) 0.44 559 0.47 0.04 0.00 low id_6
这个数据集包含 8 列。第一列 author 显示研究名称；cor 列给出辐射暴露与突变率之间未经变换的相关系数；n 表示样本量。z、se.z 与 var.z 三列分别是 Fisher-\(z\) 变换后的相关系数，以及它们的标准误和方差。radiation 列是一个调节变量，用于把效应量划分为总体辐射暴露低、中、高三个亚组。es.id 列则只是为每个效应量（也就是数据框中的每一行）提供一个唯一 ID。
这个数据集有一个很显著的特点：author 中存在重复条目。这是因为该 Meta 分析中的大多数研究都贡献了不止一个观测效应量。有些研究使用了多种突变测量方法，或者比较了不同类型的指标人群（例如暴露父母与其后代），这些情况都会导致同一研究产生多个效应。
从这种结构可以很清楚地看出，我们数据集中的效应量并不是彼此独立的。它们呈现出一种嵌套结构：多个效应量嵌套在同一项研究之内。因此，为了恰当地处理这些依赖关系，拟合一个三层 Meta 分析模型是非常合理的。
10.2.1 模型拟合
在 {metafor} 中，可以使用 rma.mv 函数来拟合三层 Meta 分析模型。下面列出这个函数最重要的一些参数，以及它们应如何指定：
yi：数据集中存放已计算效应量的那一列名称。在本例中，这一列是z，因为 Fisher-\(z\) 变换后的相关系数在数学性质上优于未经变换的相关系数。V：数据集中存放效应量 方差 的那一列名称。在这里是var.z。也可以直接使用效应量标准误的 平方，因为 \(SE_k^2 = v_k\)。slab：数据集中存放研究标签的那一列，作用类似于 {meta} 中的studlab。data：数据集名称。test：用于检验回归系数的方法。可选"z"（默认）和"t"（更推荐，使用与 Knapp-Hartung 方法相似的检验）。method：估计模型参数的方法。可以选择"REML"（推荐，限制最大似然）或"ML"（最大似然）。需要注意的是，其他类型的研究间异质性估计量（例如 Paule-Mandel）在这里并不适用。
不过，最重要的参数其实是 random，它也是最容易让人困惑的部分。在这个参数中，我们需要写一个公式，用来定义模型中的（嵌套）随机效应。对于三层模型，这个公式始终以 ~ 1 开头，然后跟上一条竖线 |。在竖线之后，我们把一个 随机效应 赋给某个分组变量（例如研究、测量工具、地区等）。这个分组变量常被称为 随机截距（random intercept），因为它告诉模型：要为每个分组假定不同的效应，也就是不同的截距。
在三层模型中，我们有两个分组变量：一个位于第 2 层，另一个位于第 3 层。我们假定这两个分组变量是嵌套关系，也就是多个第 2 层单位共同构成一个更大的第 3 层簇。
rma.mv 有一种专门的写法可以表达这种嵌套随机效应：用斜杠 / 把高层级与低层级分组变量分开。斜杠左边放第 3 层（簇）变量，右边放嵌套在其中的低层级变量。因此，这个公式的一般形式是：~ 1 | cluster/effects_within_cluster。
在我们的例子中，我们假定单个效应量（第 2 层，由 es.id 定义）嵌套在研究（第 3 层，由 author 定义）之内。因此，对应的公式是 ~ 1 | author/es.id。完整的 rma.mv 调用如下：
full.model <- rma.mv(yi = z,
V = var.z,
slab = author,
data = Chernobyl,
random = ~ 1 | author/es.id,
test = "t",
method = "REML")
我们把这个模型对象命名为 full.model。要打印结果概览，可以使用 summary 函数。
summary(full.model)
## Multivariate Meta-Analysis Model (k = 33; method: REML)
## [...]
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1788 0.4229 14 no author
## sigma^2.2 0.1194 0.3455 33 no author/es.id
##
## Test for Heterogeneity:
## Q(df = 32) = 4195.8268, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.5231 0.1341 3.9008 0.0005 0.2500 0.7963 ***
## [...]
先看 Variance Components。这里展示了模型各个层级对应的随机效应方差。第一项 sigma^2.1 表示第 3 层的 簇间 方差。在本例中，它等价于传统 Meta 分析中的研究间异质性方差 \(\tau^2\)，因为模型中的簇就是研究本身。
第二个方差分量 sigma^2.2 表示簇 内部 的方差，也就是第 2 层方差。nlvls 一列显示了每个层级上的分组数。第 3 层有 14 个分组，也就是纳入的 14 项研究；而这 14 项研究总共贡献了 33 个效应量，这一点在第二行可以看到。
在 Model Results 下，我们可以看到合并效应的估计值，也就是 \(z=0.52\)（95%CI：0.25–0.80）。为了便于解释，最好把这个效应量再转换回普通相关系数。这可以通过 {esc} 包中的 convert_z2r 函数完成：
library(esc)
convert_z2r(0.52)
## [1] 0.4777
可以看到，这对应的相关系数约为 \(r \approx 0.48\)，属于较大的相关程度。这说明，切尔诺贝利辐射暴露与突变率之间似乎存在相当强的关联。
输出中的 Test for Heterogeneity 表明，我们的数据中确实存在真实效应差异（\(p < 0.001\)）。不过，这个结果本身信息量并不算大。我们更关心的是：模型中每一个层级究竟捕捉到了多少异质性方差。换句话说，我们想知道有多少异质性来自研究 内部 差异（第 2 层），又有多少来自 研究之间 的差异（第 3 层）。
10.2.2 各层级方差的分布
我们可以通过计算多层次版本的 \(I^2\) 来回答这个问题。在传统 Meta 分析中，\(I^2\) 表示不归因于抽样误差的变异比例，也就是研究间异质性。在三层模型中，这部分异质性被拆分成两部分：一部分归因于簇 内部 的真实效应差异，另一部分归因于 簇之间 的变异。因此，这里会有两个 \(I^2\) 值，分别量化第 2 层和第 3 层所对应的总变异百分比。
var.comp函数{dmetar} 中的
var.comp函数可以用来计算多层次 \(I^2\)。如果你已经安装并加载了 {dmetar}，那么这个函数就可以直接使用。若你 没有 安装 {dmetar}，可以按下面步骤操作：
- 在线查看该函数的源代码。
- 把完整源代码复制并粘贴到 R 控制台（RStudio 左下窗格）中，然后按回车，让 R “学会”这个函数。
- 确保 {ggplot2} 已经安装并加载。
var.comp 函数只需要一个已经拟合好的 rma.mv 模型作为输入。这里我们把输出保存为 i2，然后再用 summary 函数打印结果。
i2 <- var.comp(full.model)
summary(i2)
## % of total variance I2
## Level 1 1.254966 ---
## Level 2 39.525499 39.53
## Level 3 59.219534 59.22
## Total I2: 98.75%
从输出可以看到，总方差中有多少比例归属于三个层级中的每一个。第 1 层的抽样误差方差非常小，大约只占 1%。\(I^2_{\text{Level 2}}\) 表示簇内异质性所占比例，大约为 40%。而最大的份额落在第 3 层：簇间异质性（在本例中也就是研究间异质性）约占总变异的 59%。
总体来看，这说明第 3 层存在相当显著的研究间异质性。与此同时，我们也看到，总方差中有相当大一部分，超过三分之一，可以由 研究内部 的差异来解释。
也可以把这种总方差分布可视化。只需把 var.comp 的输出传给 plot 函数即可。
plot(i2)
10.2.3 模型比较
只有当三层模型比两层模型更好地反映了数据变异时，拟合三层模型才真正有意义。如果我们发现两层模型的拟合程度与三层模型相当，那么就应当遵循 奥卡姆剃刀原则：在解释能力相同的前提下，优先选择复杂度更低的两层模型。
幸运的是，{metafor} 允许我们把三层模型与“删去一个层级”的模型进行比较。做法是再次使用 rma.mv，但这次把某个层级的方差分量固定为 0。这可以通过 sigma2 参数实现。我们需要提供一个一般形式为 c(level 3, level 2) 的向量：如果某个方差分量需要固定为 0，就写 0；如果该参数应由数据估计，则写 NA。
在本例中，我们关心的是：把单个效应量嵌套在研究之内，是否真的改善了模型。因此，我们拟合一个把第 3 层方差固定为 0 的模型，也就是把研究间异质性设为零。这相当于拟合一个简单的随机效应模型，并假定所有效应量彼此独立，而我们明知事实并非如此。既然第 3 层被固定为 0，那么 sigma2 的输入就是 c(0, NA)。我们把结果保存为 l3.removed。
l3.removed <- rma.mv(yi = z,
V = var.z,
slab = author,
data = Chernobyl,
random = ~ 1 | author/es.id,
test = "t",
method = "REML",
sigma2 = c(0, NA))
summary(l3.removed)
## [...]
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 14 yes author
## sigma^2.2 0.3550 0.5959 33 no author/es.id
##
## Test for Heterogeneity:
## Q(df = 32) = 4195.8268, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.5985 0.1051 5.6938 <.0001 0.3844 0.8126 ***
## [...]
从输出可以看出，sigma^2.1 确实被成功固定为 0。总体效应的估计值也发生了变化。那么，这个结果是否比三层模型更好呢？要评估这一点，可以使用 anova 函数比较两个模型。
anova(full.model, l3.removed)
## df AIC BIC AICc logLik LRT pval QE
## Full 3 48.24 52.64 49.10 -21.12 4195.82
## Reduced 2 62.34 65.27 62.76 -29.17 16.10 <.0001 4195.82
我们看到，相比只有两层的 Reduced 模型，Full 三层模型的拟合确实更好。该模型的 AIC 与 BIC 都更低，说明它表现更优。用于比较两模型的似然比检验（LRT）也显著（\(\chi^2_1 = 16.1, p < 0.001\)），同样支持这一结论。
换言之，尽管三层模型多引入了一个参数，也就是自由度从 2 增加到 3，但这部分额外复杂度看起来是值得的。对嵌套数据结构进行建模，很可能改善了我们对合并效应的估计。
不过请注意，即便三层模型 没有 显著优于两层模型，也常常仍然有充分理由保留三层结构。尤其是当我们相信三层结构具有坚实的理论依据时，更是如此。
例如，当数据中包含一项研究贡献多个效应量的情况时，我们 知道 这些效应不可能彼此独立。因此，保留嵌套模型仍然是合理的，因为它更贴近数据真正的“生成机制”。假如我们这个例子中的 anova 结果支持两层模型，我们大概会得出这样的结论：同一研究内部的效应总体上 大体 是同质的。但即便如此，我们很可能仍会报告三层模型的结果，因为它更准确地反映了数据生成过程。
当第 3 层簇变量的重要性并不明确时，情况就有所不同。比如说，在某个三层模型中，第 3 层簇代表不同文化区域。如果我们发现所研究的现象在文化之间根本没有差异，那么删除第三层、改用两层模型就是完全合理的。
10.3 三层模型中的亚组分析
一旦三层模型建立好，我们也可以进一步考察总体效应的潜在调节变量。前面我们已经学到，亚组分析可以表达为带有虚拟变量预测项的元回归模型。类似地，我们也可以把回归项加入“多层次”模型中，从而得到一个 三层混合效应模型：
\[\hat\theta_{ij} = \theta + \beta x_i + \zeta_{(2)ij} + \zeta_{(3)j} + \epsilon_{ij} \tag{10.8}\]其中，\(\theta\) 是截距，\(\beta\) 是预测变量 \(x\) 的回归系数。当我们把 \(x_i\) 替换为虚拟变量时，就得到一个可用于亚组分析的模型；若 \(x\) 是连续变量，上式则表示一个三层元回归模型。
在 rma.mv 中，分类或连续预测变量都可以通过 mods 参数指定。这个参数需要一个以波浪线 ~ 开头的公式，后面接预测变量名称。也可以通过同时提供多个预测变量（例如 ~ var1 + var2）来拟合多元元回归。
在我们的 Chernobyl 示例中，我们想考察：相关系数是否会随着样本总体辐射暴露水平的不同而改变（低、中或高）。这些信息就存放在数据集的 radiation 列中。下面的代码可以拟合一个三层调节变量模型：
mod.model <- rma.mv(yi = z, V = var.z,
slab = author, data = Chernobyl,
random = ~ 1 | author/es.id,
test = "t", method = "REML",
mods = ~ radiation)
summary(mod.model)
## [...]
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 28) = 0.4512, p-val = 0.6414
##
## Model Results:
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.58 0.36 1.63 0.11 -0.14 1.32
## radiationlow -0.19 0.40 -0.48 0.63 -1.03 0.63
## radiationmedium 0.20 0.54 0.37 0.70 -0.90 1.31
## [...]
首先要看的，是 Test of Moderators。我们看到 \(F_{2, 28} = 0.45\)，对应的 \(p = 0.64\)，这说明不同亚组之间并不存在显著差异。
Model Results 的呈现方式属于元回归框架，因此我们不能直接从这里读出各亚组内部的合并效应量。
第一个值，也就是截距（intrcpt），表示当总体辐射暴露水平为高时的 \(z\) 值（\(z = 0.58\)）。低暴露组和中暴露组的效应，则需要把它们各自的 estimate 加到截距上。因此，低暴露组的效应是 \(z = 0.58 - 0.19 = 0.39\)，中等暴露组的效应则是 \(z = 0.58 + 0.20 = 0.78\)。
如何报告三层（调节）模型的结果
报告三层模型结果时，除了合并效应外，至少还应报告估计得到的方差分量。
rma.mv用 \(\sigma^2_1\) 和 \(\sigma^2_2\) 分别表示第 3 层和第 2 层的随机效应方差。不过在正式撰写结果时，使用 \(\tau^2_{\text{Level 3}}\) 与 \(\tau^2_{\text{Level 2}}\) 往往更合适，因为这样更明确地表明我们讨论的是 真实（研究）效应 的方差，也就是异质性方差。如果同时报告多层次 \(I^2\) 值，也会更便于读者理解，当然前提是你先解释清楚这些指标所表示的含义。
如果你使用了
anova做模型比较，那么至少应报告似然比检验的结果。调节变量分析的结果，则可以像前面亚组分析章节中那样整理成表格。对于我们这个例子，一种可能的写法如下：“基于三层 Meta 分析模型得到的合并相关为 \(r = 0.48\)（95%CI：0.25-0.66；\(p < 0.001\)）。估计得到的方差分量分别为 \(\tau^2_{\text{Level 3}} = 0.179\) 与 \(\tau^2_{\text{Level 2}} = 0.119\)。这意味着，总变异中的 \(I^2_{\text{Level 3}} = 58.22\%\) 可归因于簇间异质性，而 \(I^2_{\text{Level 2}} = 31.86\%\) 可归因于簇内异质性。我们发现，三层模型相比于将第 3 层异质性限制为 0 的两层模型，拟合显著更优（\(\chi^2_1 = 16.10; p < 0.001\)）。”
10.4 稳健方差估计
前面几节中，我们介绍了三层 Meta 分析模型，以及如何利用它处理数据中效应量之间的依赖。刚才拟合的分层模型显然比把所有效应量都看作完全独立的“传统”Meta 分析，更能反映数据的真实结构。但它仍然只是对现实的 简化。在实际研究中，效应量之间的依赖往往比我们的嵌套模型所能捕捉的关系 更复杂。
回到 Chernobyl 数据集就能看出这一点。大多数研究都报告了多个效应量，但这些重复效应量出现的 原因 在不同研究之间并不相同。有的研究比较了不同目标人群中的辐射效应，因此报告了多个效应量；另一些研究则是在同一样本上使用了不同测量方法，因此同样产生了多个效应量。
当同一研究中的多个效应量基于同一样本时，我们通常预期它们的抽样误差会彼此 相关。但这一点还没有被前面的三层模型纳入。此前的模型实际上假定：在同一簇或研究内部，抽样误差之间的相关（从而协方差）为 0。换句话说，它假定同一研究内部的效应量估计彼此 独立。
图 10.2：原始三层（分层）模型假定研究或簇内部的效应量估计彼此独立。
因此，在这一节中，我们将介绍一种扩展的三层结构，即所谓的 相关且分层效应（CHE）模型（Correlated and Hierarchical Effects model）。和前面的分层三层模型一样，CHE 模型同样允许我们基于某种共同特征，把多个效应量归并为更大的簇，例如它们来自同一研究、同一工作团队、同一文化区域等。
但除此之外，这个模型还显式考虑到：同一簇内部有些效应量可能基于同一样本（例如进行了多次测量），因此它们的抽样误差会相关。在许多实际情形下，CHE 模型都应当被视为一个 很好的起点，尤其当数据中的依赖结构比较复杂，或者我们只知道其中的一部分时。
补充说明：Pustejovsky 和 Tipton 还提供了一棵决策树，用来判断 CHE 模型在什么情况下合适，你可以参考其论文中的 Figure 1。这棵决策树可以帮助你判断：对于当前数据，CHE 模型是否是最合理的工作模型，或者是否应选择别的依赖结构假设。
除了 CHE 模型之外，我们还会讨论在 Meta 分析中广泛使用的 稳健方差估计（Robust Variance Estimation, RVE）。这是一组在过去经常用于处理 Meta 分析中依赖效应量的方法。RVE 的核心围绕一个被称为 三明治估计量（Sandwich estimator）的工具展开。无论是结合 CHE 模型，还是结合其他 Meta 分析模型，它都可以帮助我们得到稳健的置信区间和 \(p\) 值；即使我们选择的模型并没有完美刻画数据中复杂的依赖结构，也仍然如此。
因此，在拟合第一个 CHE 模型之前，我们先整体了解一下 Meta 分析中的 RVE 和三明治估计量，并看看它为什么会有这样一个“令人垂涎”的名字。
10.4.1 三明治型方差估计量
在已发表的 Meta 分析研究中，“稳健方差估计”这个术语有时会被以一种有点奇怪的方式使用，仿佛它是一种 只适用于 带有依赖效应量的 Meta 分析数据的特殊方法。事实恰恰相反。稳健方差估计量最初是为 常规回归模型 开发的方法，用于计算回归系数 \(\hat\beta\) 的方差。
之所以称它为“稳健”估计量，是因为即便线性模型的常规假设不完全成立，它仍然能给出渐近标准误的 一致估计。这些常规假设之一就是残差方差齐性，也就是 同方差性（homoskedasticity）。而回归系数方差的稳健估计非常关键：方差估计的平方根，就是回归系数的 标准误，即 \(\sqrt{V_{\hat\beta}} = SE_{\hat\beta}\)；它直接影响到置信区间和 \(p\) 值的计算，因此会影响我们最终从模型中做出的推断。
我们这里介绍的稳健方差估计量，只是常规回归中原始方法的一个 特殊版本。Hedges、Tipton 和 Jackson 提出了一种适用于 带有依赖效应量的元回归模型 的 RVE 方法，而这一方法近年来又得到了进一步扩展。
要理解它，我们首先需要再次看一下元回归模型的公式。从概念上说，这个公式与前面章节中介绍的元回归方程非常相似。这里只是把它改写成 矩阵记号：
\[\boldsymbol{T}_j = \boldsymbol{X}_j \boldsymbol{\beta} + \boldsymbol{u}_j + \boldsymbol{e}_j \tag{10.9}\]这个公式的含义很简单：\(\boldsymbol{T}\) 中的若干效应量，是由协变量矩阵 \(\boldsymbol{X}\) 与回归系数 \(\boldsymbol{\beta}\) 共同预测得到的；除此之外，模型还包含抽样误差 \(\boldsymbol{e}_j\)，以及每项研究的随机效应 \(\boldsymbol{u}_j\)，因此它构成了一个混合效应元回归模型。
这里特别值得注意的是公式中的下标 \(j\)，以及变量全部用 粗体 表示。这说明：数据集中每个研究或簇 \(j\) 都可能贡献不止一个效应量。设某项研究 \(j\) 中共有 \(n_j\) 个效应量，那么该研究中的效应量可以写成如下列向量：\(\boldsymbol{T}_j = (T_{j,1}, \dots, T_{j,n_j})^\top\)。类似地，\(\boldsymbol{X}_j\) 是研究 \(j\) 的 设计矩阵（design matrix），其中包含了该研究内所有效应量对应的协变量取值：
\[\boldsymbol{X}_j = \begin{bmatrix} x_{1,1} & \cdots & x_{1,p} \\ \vdots & \ddots & \vdots \\ x_{n_j,1} & \cdots & x_{n_j,p} \end{bmatrix} \tag{10.10}\]其中，\(p-1\) 是协变量总数。向量 \(\boldsymbol{\beta} = (\beta_1, \dots, \beta_p)^\top\) 不带下标 \(j\)，因为我们假定这些回归系数在所有研究之间都是固定的。
补充说明：在线性回归中，设计矩阵（或模型矩阵）包含了用于估计回归系数的全部协变量取值。最简单的理解方式，是把它看作一个“协变量数据框”，只不过第一列额外加了一列全为 1 的数，用来拟合回归截距。假设我们的元回归中有三个协变量，而数据集中第 4 项研究贡献了 3 个效应量，那么它的设计矩阵可能写成：
\[\boldsymbol{X}_4 = \begin{bmatrix} 1 & 4.5 & 0 & 2 \\ 1 & 7.3 & 1 & 2 \\ 1 & 2.4 & 0 & 2 \end{bmatrix}\]总体来看，这种记号强调了一个事实：当研究可以贡献多个效应量时，我们的数据其实像是若干个较小的数据集 纵向堆叠 在一起；其中 \(J\) 表示数据中研究或簇的总数：
\[\begin{bmatrix} \boldsymbol{T}_1 \\ \boldsymbol{T}_2 \\ \vdots \\ \boldsymbol{T}_J \end{bmatrix} = \begin{bmatrix} \boldsymbol{X}_1 \\ \boldsymbol{X}_2 \\ \vdots \\ \boldsymbol{X}_J \end{bmatrix} \boldsymbol{\beta} + \begin{bmatrix} \boldsymbol{u}_1 \\ \boldsymbol{u}_2 \\ \vdots \\ \boldsymbol{u}_J \end{bmatrix} + \begin{bmatrix} \boldsymbol{e}_1 \\ \boldsymbol{e}_2 \\ \vdots \\ \boldsymbol{e}_J \end{bmatrix}. \tag{10.11}\]基于这个公式，我们可以估计元回归系数 \(\boldsymbol{\hat\beta}\)。而要计算这些系数的置信区间并进行显著性检验，我们就需要估计它们的方差 \(\boldsymbol{V_{\hat\beta}}\)。这可以通过稳健抽样方差估计量来实现，其公式如下：
\[\boldsymbol{V}^{\text{R}}_{\boldsymbol{\hat\beta}} = \left(\sum^J_{j=1}\boldsymbol{X}_j^\top\boldsymbol{W}_j\boldsymbol{X}_j \right)^{-1} \left(\sum^J_{j=1}\boldsymbol{X}_j^\top\boldsymbol{W}_j \boldsymbol{A}_j\Phi_j \boldsymbol{A}_j \boldsymbol{W}_j \boldsymbol{X}_j \right) \left(\sum^J_{j=1}\boldsymbol{X}_j^\top\boldsymbol{W}_j\boldsymbol{X}_j \right)^{-1} \tag{10.12}\]这个公式看起来确实相当复杂，而且你并不需要理解其中每一个细节。眼下真正重要的，是它的 结构 以及其中几个关键“配料”。
首先，我们看到这个公式是一个三段式结构：左边和右边的括号内容完全相同，中间夹着另一部分。它看起来就像一个三明治，外层是“面包”，中间是“夹心”，这也正是它被称为 三明治估计量 的原因。公式中真正关键的“配料”，则是矩阵 \(\boldsymbol{\Phi}_j\)、\(\boldsymbol{W}_j\) 和 \(\boldsymbol{A}_j\)：
第一个矩阵 \(\boldsymbol{\Phi}_j = \text{Var}(\boldsymbol{u}_j + \boldsymbol{e}_j)\) 是一个 \(n_j \times n_j\) 的方差-协方差矩阵。它描述的是某项研究 \(j\) 中效应量之间 真实的依赖结构。遗憾的是，我们很少真正知道同一研究内部效应量之间究竟以何种方式、在何种程度上相互相关；更不用说对 Meta 分析中的 所有 研究都了如指掌。因此，模型中通常必须做若干简化假设。举例来说，在 Hedges、Tipton 和 Jackson 的原始做法中，\(\boldsymbol{\Phi}_j\) 会被模型残差的叉积 \(\boldsymbol{e}_j \boldsymbol{e}_j^\top\) 所替代，因为这部分信息是现成可得的。它当然只是对真实依赖结构的一个粗略估计，但当 Meta 分析中的研究数足够多时，可以作为一种“最佳猜测”。而 CHE 模型则假定：同一研究内效应量之间存在一个 已知相关系数 \(\rho\)，而且这一 \(\rho\) 在各研究内、各研究之间都保持相同，这就是所谓的 “constant sampling correlation” 假设。
\(\boldsymbol{W}_j\) 矩阵包含每个效应量的 权重。前面几章中我们已经学过，只有先考虑效应量估计的精确性，才可能合理地对它们进行合并。最优的做法是取方差的倒数，也就是令 \(\boldsymbol{W}_j = \boldsymbol{\Phi}^{-1}_j\)。但正如前面提到的，\(\boldsymbol{\Phi}_j\) 的真实值几乎从来不可得，因此实践中使用的是基于模型得到的估计值 \((\boldsymbol{\hat\Phi}_j)^{-1}\)。补充说明：Hedges、Tipton 和 Jackson 的原始方法使用的是另一套近似有效的、经过简化的对角权重矩阵。
最后一部分 \(\boldsymbol{A}_j\) 是一个 调整矩阵。它的作用是：即便 Meta 分析中研究数量较少（比如 40 项或更少），估计量仍能给出有效结果。通常推荐采用基于 bias-reduced linearization 的 CR2 方法。CR2 调整矩阵可写为：
这个表达式包含了对权重矩阵 \(\boldsymbol{W}_j\) 的对称平方根计算。
10.4.2 在稳健方差估计下拟合 CHE 模型
现在，是时候在 R 中拟合我们的第一个相关且分层效应模型，并利用稳健方差估计来降低模型设定错误的风险了。和前面一样，我们仍然使用 {metafor} 中的 rma.mv 函数来拟合模型。不过，这次还需要 {clubSandwich} 包提供的一些附加函数，因此请先安装并加载该包。
library(clubSandwich)
正如前面所说，CHE 模型假定：同一研究或簇内部的效应量彼此 相关，并且这种相关在研究内部与研究之间都是相同的。
因此，我们首先要 指定一个相关系数 供模型使用。在 Chernobyl 数据中，我们暂时假定这种相关较高，因此设定 \(\rho = 0.6\)。这只是一个猜测，实际分析中强烈建议针对不同的 \(\rho\) 值进行 敏感性分析。
# 常数抽样相关假定
rho <- 0.6
在这个相关系数的基础上，我们就可以为每一项研究计算一个假定的方差-协方差矩阵。这可以通过 {clubSandwich} 中的 impute_covariance_matrix 函数实现：
vi参数用于指定数据集中存放各效应量 方差 的变量名，也就是标准误的平方。cluster参数用于指定把每个效应量归属于某个 研究 或 簇 的变量。在Chernobyl数据集中，这个变量是author。r参数则接收我们所假定的、效应量之间 恒定不变的相关系数。
# 常数抽样相关工作模型
V <- with(Chernobyl,
impute_covariance_matrix(vi = var.z,
cluster = author,
r = rho))
准备好 V 中的方差-协方差矩阵后，就可以拟合 rma.mv 模型了。假设我们要分析与第 10.3 节相同的元回归模型，也就是把 radiation 作为协变量。
这一次，参数写法稍有变化：第一个参数是一个 formula 对象，我们在其中声明效应量 z 由截距（1）和协变量 radiation 来预测；V 参数接收刚刚建立好的方差-协方差矩阵列表；sparse 参数设为 TRUE 可以加快计算速度。
只有 random 和 data 这两个参数与前面保持一致。我们把结果保存为 che.model。
che.model <- rma.mv(z ~ 1 + radiation,
V = V,
random = ~ 1 | author/es.id,
data = Chernobyl,
sparse = TRUE)
要计算元回归系数的 置信区间，可以使用 {clubSandwich} 中的 conf_int 函数。我们只需要提供已拟合的模型，并在 vcov 中指定所采用的小样本调整方式。按照推荐做法，这里使用 "CR2" 调整（见第 10.4.1 节）。
conf_int(che.model,
vcov = "CR2")
## Coef. Estimate SE d.f. Lower 95% CI Upper 95% CI
## intrcpt 0.584 0.578 1.00 -6.76 7.93
## radiationlow -0.190 0.605 1.60 -3.52 3.14
## radiationmedium 0.207 0.603 1.98 -2.41 2.83
可以看到，Estimate 下的点估计与我们在第 10.3 节得到的结果大体相似；但估计的标准误与置信区间却大得多。若要进一步计算回归系数的 \(p\) 值，可以使用 coef_test 函数：
coef_test(che.model,
vcov = "CR2")
## Coef. Estimate SE t-stat d.f. (Satt) p-val (Satt) Sig.
## intrcpt 0.584 0.578 1.010 1.00 0.497
## radiationlow -0.190 0.605 -0.315 1.60 0.789
## radiationmedium 0.207 0.603 0.344 1.98 0.764
我们看到，在使用稳健方差估计之后，这些系数都不再显著。默认情况下，conf_int 和 coef_test 使用的是经过 Satterthwaite 修正 的自由度，这通常也是推荐保留的默认设置。
稳健方差估计与模型设定错误
有些读者可能会问：为什么我们如此强调要为模型使用稳健方差估计？最主要的原因在于，多变量模型和多层模型都 很容易设定错误。我们已经看到，即便是 CHE 模型，也只是通过假定研究内部与研究之间的相关完全相同，来对现实进行近似；而在很多时候，我们其实并不清楚所建立的模型是否 足够合理地逼近了 数据中复杂的依赖结构。
稳健方差估计的价值正在于此：即便模型本身可能存在某种程度的设定错误，它仍然能够在一定程度上保护我们的 统计推断，也就是保护我们计算出的置信区间与 \(p\) 值不至于因为模型失配而过于乐观。
{robumeta} 包
在这一节里，我们讨论的是稳健方差估计与相关且分层效应模型的结合。这一模型连同其他一些扩展，是由 Pustejovsky 和 Tipton 提出的。
Hedges、Tipton 与 Jackson 最初提出的 “原始版”RVE 方法，以及一些小样本扩展，则可以通过 {robumeta} 包来实现。该包允许我们拟合 Hedges、Tipton 和 Jackson 最早提出的两类模型：分层效应模型（hierarchical effects model）和 相关效应模型（correlated effects model），但不能把二者同时结合起来。
10.5 聚类 Wild Bootstrap
在上一节中，我们已经看到如何拟合相关且分层效应模型，以及如何利用稳健方差估计来计算置信区间和系数检验。
另一种有时更有优势的系数检验方式，是基于 bootstrap 的程序；其中一种特殊形式就叫做 聚类 Wild Bootstrap。当 Meta 分析中的研究总数 \(J\) 较少 时，这种方法尤其合适；相较之下，RVE 在小样本情况下可能会过于保守，这一点我们在 Chernobyl 例子里已经见过。
当我们需要检验所谓的 多重对比假设（multiple-contrast hypotheses）时，这种方法也非常有用。例如，如果我们想检验某个经过虚拟编码的分类协变量的整体效应，就需要使用这样的假设检验。
聚类 Wild Bootstrap 算法
Wild bootstrap 是一种基于零模型残差的重抽样方法。这里所说的零模型，是指不含额外协变量的模型。在聚类 wild bootstrap 中，为了处理依赖效应量，需要先用某个调整矩阵 \(\boldsymbol{A}_j\)（例如基于 CR2 的调整矩阵；见第 10.4.1 节）来变换残差。其一般算法如下：
- 基于原始数据拟合完整模型，并计算感兴趣的检验统计量（例如 \(t\) 值或 \(F\) 值）。
- 基于原始数据拟合零模型，并提取其残差 \(\boldsymbol{e}\)。
- 对每个研究或簇 \(j\)，从某个分布中抽取一个随机值；再用这个随机值去乘该簇的残差。在本章使用的 {wildmeta} 包中，这里采用的是 Rademacher 分布。
- 将变换后的残差加到零模型基于原始数据得到的预测值上，从而生成新的 bootstrap 效应量。
- 用这些 bootstrap 效应量重新拟合完整模型，并再次计算检验统计量。
第 3 步到第 5 步会重复 \(R\) 次。最终的 bootstrap \(p\) 值，就是 bootstrap 检验统计量比原始数据对应统计量 更极端 的比例。
为了用 bootstrap 检验多重对比假设，我们可以使用 {wildmeta} 包。这个包需要先安装并加载。除此之外，我们还会用到 {tidyverse} 中的函数，在 Chernobyl 数据集中新增一个变量，用来保存每项研究的发表年份。
# 确保已经加载 {wildmeta} 和 {tidyverse}
library(wildmeta)
library(tidyverse)
# 新增年份变量
Chernobyl$year <- str_extract(Chernobyl$author,
"[0-9]{4}") %>% as.numeric()
接下来，我们把这个变量作为 新的预测变量 加入 rma.mv 元回归模型。做法很简单：在第一个参数的公式里把 year 加进去，并用 scale 函数对它做中心化与标准化。这里还有一个小变化：我们把截距中的 1 改成了 0。这意味着模型 不再包含截距项，从而使 year 的预测在不同 radiation 水平下 分层呈现。这样做的结果是：我们不再得到“以某一组为参照”的回归系数，而是得到三个分别对应不同 radiation 水平的合并效应估计。我们把结果保存为 che.model.bs。
che.model.bs <- rma.mv(z ~ 0 + radiation + scale(year),
V = V,
random = ~ 1 | author/es.id,
data = Chernobyl,
sparse = TRUE)
在正式开始 bootstrap 之前，我们还需要为想要检验的假设定义一个 线性对比。假设我们想检验 radiation 变量的 整体调节效应。这时，可以使用 {clubSandwich} 中的 constrain_equal 函数，为检验创建一个约束矩阵。原假设是：radiation 的三个水平对应的效应彼此相等，因此把 constraints 设为 1:3。同时，我们还要通过 coefs 参数传入刚刚拟合模型的系数。结果保存为 rad.constraints。
rad.constraints <- constrain_equal(constraints = 1:3,
coefs = coef(che.model.bs))
rad.constraints
## [,1] [,2] [,3] [,4]
## [1,] -1 1 0 0
## [2,] -1 0 1 0
现在，我们就可以使用 {wildmeta} 中的 Wald_test_cwb 函数来计算这一多重对比假设的 bootstrap \(p\) 值了。我们需要提供完整模型、约束矩阵、想要使用的小样本调整类型，以及 R（bootstrap 重复次数）。通常建议使用 较高的重复次数（例如 1000 次或更多），因为这样可以提高检验效能。在这个例子里，我们使用 2000 次重复，并把结果保存为 cw.boot。请注意，根据迭代次数的多少，这个过程 可能需要几分钟 才能完成。
cw.boot <- Wald_test_cwb(full_model = che.model.bs,
constraints = rad.constraints,
adjust = "CR2",
R = 2000)
cw.boot
## Test Adjustment CR_type Statistic R p_val
## 1 CWB Adjusted CR2 CR0 Naive-F 2000 0.3595
可以看到，我们对调节效应所做检验的 \(p\) 值为 0.36，结果 不显著。这与我们在第 10.3 节中对辐射强度所做的调节变量分析结论是一致的。
此外，还可以用 plot 函数把所有 bootstrap 重复中的 检验统计量密度分布 可视化出来。
plot(cw.boot,
fill = "lightblue",
alpha = 0.5)
10.6 问答
检验一下你的掌握程度！
- 为什么说“三层模型”比“多层次模型”这个说法更准确？
- 在什么情况下，三层 Meta 分析模型特别有用？
- 请说出效应量依赖的两个常见来源。
- 多层次 \(I^2\) 统计量应如何解释？
- 如何把三层模型扩展为能够纳入调节变量效应的模型？
这些问题的答案列在本书末尾的 附录 A 中。
10.7 小结
所有随机效应 Meta 分析都建立在一个多层次模型之上。当再增加一层时，我们就得到三层 Meta 分析模型。这类模型尤其适合处理 成簇（clustered）的效应量数据。
三层模型可用于处理 相互依赖 的效应量。例如，当一项研究贡献多个效应量时，我们通常不能再假定这些结果彼此独立。三层模型通过假定效应量 嵌套 在更大的簇中（例如研究）来处理这一问题。
与传统 Meta 分析不同，三层模型需要估计两个异质性方差：一个是簇 内部 的随机效应方差，另一个是 簇之间 的异质性方差。
我们同样可以使用三层模型来检验分类或连续预测变量。这样得到的就是三层混合效应模型。