Meta Regression

Published: April 17, 2026

8 元回归

来源页面：https://doing-meta.guide/metareg.html

上一章中，我们把亚组分析加入了 Meta 分析“工具箱”，作为一种新的方法。正如我们已经学到的，亚组分析会把分析重心从寻找单一的总体效应，转移到考察数据中的异质性模式及其成因上来。

我们也提到过，亚组分析其实是 元回归（meta-regression）的一种特殊形式。你很可能早就听说过“回归”这个术语。回归分析是最常见的统计方法之一，广泛应用于各个学科。在最简单的形式下，回归模型试图利用某个变量 $x$ 的取值来预测另一个变量 $y$ 的取值。通常，回归模型基于个体层面的数据建立，在这些数据中，$x$ 与 $y$ 的值都会被测量。

在元回归中，这一逻辑被应用到 整项研究 上。变量 $x$ 代表研究的某些特征，例如研究实施的年份。基于这些信息，元回归模型尝试预测 $y$，也就是该研究的效应量。不过，把效应量作为被预测变量会带来一些额外复杂性。

我们在第 3.1 节已经学过，观测到的效应量 $\hat\theta$ 作为研究真实效应的估计，其 精确性 会因标准误不同而有所差异。在“常规”的 Meta 分析中，我们通过赋予研究不同权重来处理这一点。而在元回归中，我们同样必须确保模型会更多关注抽样误差更小的研究，因为我们可以假定它们的估计值更接近“真实值”。

元回归正是通过假定 混合效应模型（mixed-effects model）来实现这一点的。该模型考虑到：观测研究之所以偏离真实总体效应，一方面是由于抽样误差，另一方面是由于研究间异质性。更重要的是，它还使用一个或多个变量 $x$ 来预测真实效应量之间的差异。我们在上一章已经提到，亚组分析同样建立在混合效应模型之上。本章将进一步深入讨论，解释为什么亚组分析与元回归在本质上是紧密相关的。

尽管元回归也有其局限，但它在 Meta 分析中仍然是一个非常有力的工具。它也非常灵活。例如，多元元回归（multiple meta-regression）允许我们不仅纳入一个预测变量，还能同时纳入多个预测变量以及它们之间的交互作用。因此，在本章后半部分，我们也将看看多元元回归，以及如何在 R 中实现它。

8.1 元回归模型

过去，你可能已经用原始研究数据做过回归分析，在那种情况下，参与者是分析单位。而在 Meta 分析中，通常拿不到每位参与者的个体数据，我们只能依赖汇总结果。因此，在元回归中，我们必须使用 研究层面 的预测变量。

这也意味着，尽管我们分析的总样本量往往比单项原始研究更大，但可用于元回归的数据点仍然可能不足，从而使元回归失去实际意义。我们在第 7.2 节已经提到，当 $K < 10$ 时，亚组分析通常没有太大意义。Borenstein 等人（2011，第 20 章）指出，这一经验法则同样可以应用于元回归模型，但它不应被视为一条绝对铁律。

在常规回归中，我们希望利用预测变量（或 协变量）$x_i$ 及其回归系数 $\beta$，来估计个体 $i$ 的 $y_i$ 值。因此，标准回归方程写作：

\[\hat{y_i} = \beta_0 + \beta_1x_i \tag{8.1}\]

在元回归中，我们要预测的变量 $y$ 是研究 $k$ 的观测效应量 $\hat\theta_k$。元回归 的公式与常规回归模型很相似：

\[\hat\theta_k = \theta + \beta x_{k} + \epsilon_k+\zeta_k \tag{8.2}\]

注意，这个公式中多了两项：$\epsilon_k$ 和 $\zeta_k$。这两项同样出现在随机效应模型的方程中（见第 4.1.2 节），表示两类彼此独立的误差。第一项 $\epsilon_k$ 是抽样误差，它使某项研究的效应量偏离其真实效应。

第二项 $\zeta_k$ 则表示，即便是研究的真实效应量，也只是从一个更高层次的效应量分布中抽取出来的。这意味着我们的数据中存在研究间异质性，而它由异质性方差 $\tau^2$ 来刻画。

由于上式同时包含固定效应（$\beta$ 系数）和随机效应（$\zeta_k$），因此元回归所使用的模型通常被称为 混合效应模型。从概念上说，这与我们在第 7.1.2 节解释亚组分析时介绍的混合效应模型是完全一致的。

8.1.1 带分类预测变量的元回归

正如前面所说，亚组分析本质上就是带有分类预测变量的元回归。这类分类变量可以通过 虚拟编码（dummy-coding）纳入，例如：

\[D_g=\begin{cases} 0: & \text{亚组 A}\\ 1: & \text{亚组 B} \end{cases} \tag{8.3}\]

如果要把亚组分析写成元回归形式，我们只需用 $D_g$ 替代协变量 $x_k$：

\[\hat\theta_k = \theta + \beta D_g + \epsilon_k + \zeta_k \tag{8.4}\]

要理解这个公式，我们需要从左往右读。元回归模型和其他统计模型一样，目标都是解释观测数据是如何生成的。对我们来说，这个观测数据就是 Meta 分析中某项研究 $k$ 的观测效应量 $\hat\theta_k$。上述公式就像一份“配方”，告诉我们要得到这个观测效应，需要哪些成分。

首先是 $\theta$，它在回归模型中充当截距。$\theta$ 的值等同于亚组 A 的真实总体效应量。要理解原因，就要看下一个“成分” $\beta D_g$。在这一项中，$\beta$ 表示亚组 A 与亚组 B 之间的效应量差异 $\theta_{\Delta}$。$\beta$ 与 $D_g$ 相乘，而 $D_g$ 只可能是 0 或 1，分别对应研究属于亚组 A（$D_g = 0$）或亚组 B（$D_g = 1$）。

由于乘以 0 的结果为 0，所以当研究属于亚组 A 时，$\beta D_g$ 这一项会完全消失。反之，当 $D_g = 1$ 时，这一项就等于 $\beta$，并与 $\theta$ 相加，从而给出亚组 B 的总体效应量。本质上，虚拟预测变量是一种把两个公式整合为一个公式的方法。把公式分别写成两个亚组的形式后，这一点会更清楚：

\[D_g=\begin{cases} 0: & \text{$\hat\theta_k = \theta_A + \epsilon_k+\zeta_k$}\\ 1: & \text{$\hat\theta_k = \theta_A + \theta_{\Delta} + \epsilon_k+\zeta_k$} \end{cases} \tag{8.5}\]

这样写以后，我们会更清楚地看到，这个公式实际上包含了两个模型：一个对应亚组 A，一个对应亚组 B。二者最主要的区别在于，第二个亚组的效应会根据 $\beta$ 的取值向上或向下“平移”（在上式中我们将其记作 $\theta_{\Delta}$）。

图 8.1：带分类预测变量（亚组分析）的元回归。

由此可见，亚组分析的工作方式和普通回归非常相似：它使用某个变量 $x$ 来预测 $y$ 的取值，而在这里，$y$ 是研究的效应量。特殊之处在于，$\beta x_k$ 不是连续变化的，而是一个固定加到预测值上的数，它取决于研究是否属于某个亚组。这个固定的 $\beta$ 值，就是两个亚组之间估计得到的效应量差异。

8.1.2 带连续预测变量的元回归

不过，当人们提到“元回归”时，通常首先想到的是把连续变量作为预测变量的模型。这就回到了公式 8.2 中所示的通用元回归公式。在这种情况下，我们刚才讨论过的回归项依然存在，但其作用稍有不同。$\theta$ 仍然表示截距，不过此时它代表的是当 $x = 0$ 时的预测效应量。

在截距的基础上，再加上 $\beta x_k$ 这一项。它形成了 回归斜率：连续变量 $x$ 与 回归权重 $\beta$ 相乘，从而使协变量不同取值对应的预测效应增大或减小。

元回归模型的目标，是找到合适的 $\theta$ 和 $\beta$，使预测效应量与研究真实效应量之间的差异最小（见图 8.2）。

图 8.2：带连续预测变量和四项研究的元回归。

仔细观察元回归公式会发现，其中包含两类项。有些项带有下标 $k$，有些则没有。带下标 $k$ 表示该值会随研究不同而变化；如果某一项不带下标 $k$，就表示它在所有研究中都保持不变。

（原页面此处展示了一个公式构成示意图，用于区分随研究变化与不变化的项。）

在元回归中，$\theta$ 和 $\beta$ 都是固定不变的。这个特征告诉我们元回归到底在做什么：它根据预测变量的变化以及观测到的效应，试图从数据中“提炼”出一个潜在的 固定模式，即一条 回归线。如果元回归模型拟合良好，那么估计得到的参数 $\theta$ 和 $\beta$ 就可以用来预测模型 从未见过 的研究的效应量（前提是我们知道该研究的 $x$ 值）。

因此，在同时考虑抽样误差 $\epsilon_k$ 和研究间异质性 $\zeta_k$ 的情况下，元回归试图找到一个能够良好泛化的模型；它不仅适用于已经观察到的效应量，也适用于所有我们感兴趣的潜在研究总体。

8.1.3 模型拟合的评估

关于元回归模型，一个重要细节是：它可以被看作我们用于合并效应量的“常规”随机效应模型的扩展。随机效应模型其实就是 没有斜率项 的元回归模型。因为它不包含斜率，所以随机效应模型对每项研究都只预测 同一个值，也就是合并效应量估计 $\mu$，这与截距是等价的。

因此，元回归的计算第一步与随机效应 Meta 分析非常相似：首先用前文介绍过的某种方法来估计研究间异质性 $\tau^2$（例如 DerSimonian-Laird 或 REML 方法）。第二步再估计固定权重 $\theta$ 和 $\beta$。常规线性回归使用 普通最小二乘法（ordinary least squares, OLS）寻找最佳拟合回归线；而元回归使用一种修正后的方法，即 加权最小二乘法（weighted least squares, WLS），以确保标准误更小的研究被赋予更高权重。

一旦找到最优解，我们就可以检查新加入的回归项是否解释了部分效应量异质性。如果元回归模型拟合良好，那么真实效应量相对于回归线的偏离，应当小于它们相对于合并效应 $\hat\mu$ 的偏离。如果确实如此，就说明预测变量 $x$ 解释了 Meta 分析中部分异质性方差。

（原页面此处展示了随机效应模型与混合效应模型异质性解释差异的示意图。）

因此，我们可以通过检查元回归模型解释了多少异质性方差，来评估模型拟合。混合效应模型中所包含的预测变量，应当尽量减小残余、也就是未被解释的异质性方差，我们将其记为 $\tau^2_{\text{unexplained}}$。

在回归分析中，通常使用 $R^2$ 指标量化模型解释的变异比例。对于元回归，也可以计算一个类似的指标 $R^2_{*}$。这里加上星号，是为了说明元回归中的 $R^2$ 与常规回归中的 $R^2$ 略有不同，因为这里处理的是 真实效应量，而不是观测数据点。$R^2_*$ 的公式为：

\[R^2_* = 1 - \frac{\hat\tau^2_{\text{unexplained}}}{\hat\tau^2_{\text{(total)}}} \tag{8.6}\]

$R^2_*$ 使用的是即使加入元回归斜率之后仍然无法解释的残余异质性方差，并将其与我们在最初 Meta 分析中发现的总异质性相比较。用 1 减去这一比值，就能得到预测变量解释的研究间异质性百分比。

$R^2_*$ 还可以换一种写法。我们也可以说，它表示与初始随机效应合并模型相比，混合效应模型使异质性方差 减少了多少百分比。于是可写为：

\[R^2_* = \frac{\hat\tau^2_{\text{REM}}-\hat\tau^2_{\text{MEM}}}{\hat\tau^2_{\text{REM}}} \tag{8.7}\]

在这个公式中，$\hat\tau^2_{\text{REM}}$ 表示随机效应合并模型中发现的研究间异质性，而 $\hat\tau^2_{\text{MEM}}$ 表示混合效应元回归模型中的（残余）方差，也就是相对于真实效应量的“预测误差”。

通常，我们不仅关心回归模型解释了多少异质性，也关心预测变量 $x$ 的回归权重是否显著。如果显著，我们就可以较为有把握地认为，$x$ 会影响研究的效应量。无论是常规回归还是元回归，回归权重显著性最常见的检验方法都是 Wald 型检验。它通过将 $\beta$ 的估计值除以其标准误来计算检验统计量 $z$：

\[z = \frac{\hat\beta}{SE_{\hat\beta}} \tag{8.8}\]

在原假设 $\beta = 0$ 下，这个 $z$ 统计量服从标准正态分布。于是我们就可以计算相应的 $p$ 值，并据此判断该预测变量是否显著。

不过，基于 $z$ 统计量的检验并不是唯一的做法。与常规 Meta 分析模型类似，我们也可以使用 Knapp-Hartung 调整，此时得到的是基于 $t$ 分布的检验统计量。正如前文所学，Knapp-Hartung 方法往往更值得推荐，因为它能够降低假阳性的风险。

8.2 在 R 中进行元回归

{meta} 软件包中包含一个名为 metareg 的函数，它允许我们进行元回归。metareg 函数只需要一个 {meta} 的 Meta 分析对象，以及一个协变量名称作为输入。

在这个例子中，我们再次使用基于 ThirdWave 数据集构建的 Meta 分析对象 m.gen（见第 4.2.1 节）。通过元回归，我们想检验研究的 发表年份 是否能够预测其效应量。默认情况下，ThirdWave 数据集中并没有专门存储发表年份的变量，因此我们必须自己创建一个包含这些信息的 numeric 变量。我们只需按照研究在 ThirdWave 数据集中出现的顺序，把所有研究的发表年份连接起来，并将这个变量命名为 year。

注：本例中使用的发表年份是虚构的，仅用于演示。

year <- c(2014, 1998, 2010, 1999, 2005, 2014, 
          2019, 2010, 1982, 2020, 1978, 2001,
          2018, 2002, 2009, 2011, 2011, 2013)

现在，我们已经具备运行元回归所需的全部信息。在 metareg 函数中，我们把 Meta 分析对象 m.gen 作为第一个参数，把预测变量 year 作为第二个参数，并把结果保存为 m.gen.reg。

m.gen.reg <- metareg(m.gen, ~year)

现在来看结果：

m.gen.reg

## Mixed-Effects Model (k = 18; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.019 (SE = 0.023)
## tau (square root of estimated tau^2 value):             0.1371
## I^2 (residual heterogeneity / unaccounted variability): 29.26%
## H^2 (unaccounted variability / sampling variability):   1.41
## R^2 (amount of heterogeneity accounted for):            77.08%
## 
## Test for Residual Heterogeneity:
## QE(df = 16) = 27.8273, p-val = 0.0332
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 16) = 9.3755, p-val = 0.0075
## 
## Model Results:
## 
##         estimate     se   tval   pval    ci.lb    ci.ub 
## intrcpt   -36.15  11.98  -3.01  0.008  -61.551  -10.758  ** 
## year        0.01   0.00   3.06  0.007    0.005    0.031  ** 
## 
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

下面我们逐项看一下这些结果。第一行告诉我们，数据上确实拟合了一个混合效应模型，这与我们的预期一致。接下来的几行给出了模型解释异质性的程度。我们看到，残余异质性方差，也就是未被预测变量解释的方差，其估计值为 $\hat\tau^2_{\text{unexplained}} = 0.019$。

输出中还给出了一个相当于 $I^2$ 的指标，显示在纳入预测变量之后，我们数据中 29.26% 的变异仍可归因于剩余的研究间异质性。在常规随机效应 Meta 分析模型中，我们发现 $I^2$ 异质性为 63%，这意味着该预测变量能够“解释掉”真实效应量差异中的相当大一部分。

在最后一行，我们看到 $R^2_*$ 的值，在本例中为 77%。这意味着真实效应量差异中有 77% 可以由发表年份解释，这是一个相当可观的比例。

接下来是 Test for Residual Heterogeneity，它本质上就是我们之前已经接触过的 $Q$ 检验。只不过这里检验的是：预测变量未能解释掉的那部分异质性是否仍然显著。结果显示显著（$p = 0.03$）。不过，我们已经了解过 $Q$ 检验的局限，因此不应过度依赖这一结果。

下一部分是 Test of Moderators。我们看到这个检验同样显著（$p = 0.0075$），这意味着我们的预测变量，也就是发表年份，确实会影响研究的效应量。

最后一部分给出了估计回归系数的更多细节。第一行是截距（intrcpt）的结果。这表示当预测变量“发表年份”为 0 时的预期效应量（在这里即 Hedges’ $g$）。在本例中，这对应的是一个多少有些人为构造的情景：它表示公元 0 年开展的一项研究的预测效应，其值为 $\hat{g} = -36.15$。这再次提醒我们：好的统计模型并不需要完美地还原现实，它们只需要有用即可。

我们真正关心的是第二行中的系数。模型对 year 的回归权重估计为 0.01。这意味着每增加一年，研究的效应量 $g$ 预期会上升 0.01。因此，我们可以说，研究的效应量会随着时间推移而增加。其 95% 置信区间为 0.005 到 0.031，说明该效应是显著的。

重要的是，输出中还给出了每个回归系数对应的 $t$ 统计量（tval）。这说明模型使用了 Knapp-Hartung 方法来计算置信区间和 $p$ 值。由于我们在最初的 Meta 分析模型中也使用了这一调整，因此 metareg 在这里会自动继续使用它。否则，输出中展示的将会是 $z$ 值和 Wald 型置信区间。

气泡图（bubble plot）

{meta} 软件包允许我们用 bubble 函数对元回归结果进行可视化。它会生成一个 气泡图，其中展示了估计得到的回归斜率，以及每项研究的效应量。为了反映研究权重，气泡大小会有所不同，气泡越大表示该研究的权重越高。

如果要绘制气泡图，我们只需要把元回归对象传给 bubble 函数即可。由于我们还希望显示研究标签，所以把 studlab 设为 TRUE。

bubble(m.gen.reg, studlab = TRUE)

（原页面此处展示了对应的 bubble plot。）

为了完整起见，我们也可以把上一章的亚组分析（见第 7 章）放到元回归框架中重做一次。也就是说，把偏倚风险评估作为分类预测变量。由于 RiskOfBias 这个变量已经包含在 ThirdWave 数据集中，所以我们不需要再把这部分信息单独保存成一个对象。我们只需再次运行 metareg 函数，并把 RiskOfBias 作为第二个函数参数。

metareg(m.gen, RiskOfBias)

## [...]
## R^2 (amount of heterogeneity accounted for):            15.66%
## 
## Test for Residual Heterogeneity:
## QE(df = 16) = 39.3084, p-val = 0.0010
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 16) = 2.5066, p-val = 0.1329
## 
## Model Results:
## 
##               estimate    se   tval    pval  ci.lb ci.ub 
## intrcpt           0.76  0.15   5.00  0.0001   0.44  1.09  *** 
## RiskOfBiaslow    -0.29  0.18  -1.58  0.1329  -0.69  0.10      
## [...]

从输出中可以看到，$R^2_*$ 为 15.66%，显著小于 year 这一预测变量。与我们之前的结果一致，偏倚风险变量并不是显著的效应量预测变量（$p = 0.13$）。

在 Model Results 下，我们还可以看到，metareg 自动把 RiskOfBias 转换成了虚拟变量。截距估计值代表“高风险”亚组的合并效应，其值为 $g = 0.76$。代表低偏倚风险研究的回归系数估计值为 -0.29。

要得到这个亚组的效应量，我们需要把回归权重加到截距上，即 $g = 0.76 - 0.29 \approx 0.47$。这个结果与假定亚组间共享同一个 $\tau^2$ 估计值的亚组分析结果完全一致。

8.3 多元元回归

此前，我们只讨论了元回归模型中使用一个预测项 $\beta x_k$ 的情形。在前面的例子里，我们检验的是研究效应量是否取决于发表年份。现在，假设报告的效应量还取决于研究发表期刊的声望。我们认为，那些发表在高声誉期刊上的研究，可能会报告更高的效应。这也许是因为高声誉期刊筛选更严格，主要发表具有“突破性”发现的研究。

另一方面，也完全可能是高声誉期刊通常发表的是 质量更高 的研究。也许真正与更高效应量相关的，只是更好的研究质量。于是，要检验期刊声誉是否真的与更高效应有关，我们就必须确保这种关系没有受到混杂因素的影响，也就是没有被“高声誉期刊更可能发表高质量证据”这一事实所混杂。这意味着，在考察期刊声望与效应量的关系时，我们需要对研究质量进行控制。

这类问题，以及许多其他研究问题，都可以通过 多元元回归 来处理。在多元元回归中，我们使用多个预测变量，而不是只用一个预测变量，来解释效应量的变异。为了允许多个预测变量，我们需要把之前的元回归公式（见公式 8.2）修改为：

\[\hat \theta_k = \theta + \beta_1x_{1k} + \dots + \beta_nx_{nk} + \epsilon_k + \zeta_k \tag{8.10}\]

这个公式表明，我们可以向元回归模型中再加入 $n-1$ 个预测变量 $x$，从而将其扩展为多元元回归。公式中的省略号表示，理论上我们可以加入任意多个预测变量。但在现实中，事情通常要复杂得多。接下来，我们将讨论多元元回归中的几个重要陷阱，以及如何建立稳健而可信的模型。不过在此之前，我们先介绍多元元回归的另一个重要特征：交互作用。

8.3.1 交互作用

到目前为止，我们讨论的都是模型中存在多个预测变量 $x_1, x_2, \dots, x_n$ 的情况，它们及其回归权重 $\beta$ 被加总在一起。然而，多元元回归模型并不限于这种加性关系。它还可以刻画预测变量之间的 交互作用。所谓交互作用，是指某个预测变量（例如 $x_1$）与估计效应量之间的关系，会随着另一个协变量（例如 $x_2$）的不同取值而 发生变化。

设想我们要建模两个预测变量及其与效应量的关系：研究的发表年份（$x_1$）和研究质量（$x_2$）。研究质量可以这样编码：

\[x_2=\begin{cases} 0: & \text{低}\\ 1: & \text{中等}\\ 2: & \text{高} \end{cases} \tag{8.11}\]

如果我们假设发表年份与研究质量之间不存在交互作用，那么就可以给 $x_1$ 和 $x_2$ 各自赋予一个回归权重 $\beta$，并把它们相加到公式中：

\[\hat \theta_k = \theta + \beta_1x_{1k} + \beta_2x_{2k} + \epsilon_k + \zeta_k \tag{8.12}\]

但如果 $x_1$ 和 $x_2$ 之间的关系更复杂呢？正如前面的例子，较新的发表年份可能与更高的效应正相关，但并不是所有研究都一定遵循这一趋势。也许这种增加趋势在高质量研究中最为明显，而低质量研究的结果则基本没有随时间变化。我们可以把这种关于效应量（$\hat\theta_k$）、发表年份（$x_1$）与研究质量（$x_2$）之间关系的假设可视化如下：

（原页面此处展示了发表年份、研究质量与效应量之间交互关系的示意图。）

这张图展示了交互作用的经典例子。我们可以看到，回归斜率的陡峭程度取决于另一个预测变量的取值。高质量研究对应的斜率非常陡，表明年份与效应之间关系很强；而在低质量研究中，情况则不同。该亚组中的回归线几乎是水平的，这说明发表年份对结果没有影响，甚至可能略有负向影响。

这个例子体现了交互作用的一项优势：它能够帮助我们检验某个预测变量的影响是否在所有研究中都相同，还是会受到另一种研究特征的调节。

如果要通过元回归评估交互作用，我们就需要在模型中加入一个 交互项。在本例中，可以通过加入第三个回归权重 $\beta_3$ 来实现，它刻画的是我们要检验的交互项 $x_{1k}x_{2k}$。于是公式变为：

\[\hat \theta_k = \theta + \beta_1x_{1k} + \beta_2x_{2k} + \beta_3x_{1k}x_{2k}+ \epsilon_k + \zeta_k \tag{8.13}\]

尽管线性多元元回归模型只由这些简单的构件组成，但它们的应用非常广泛。不过，在开始用 R 拟合多元元回归之前，我们还应先考虑它的局限与潜在陷阱。

8.3.2 多元元回归中的常见陷阱

如果运用得当，多元元回归会非常有用，但它也伴随着一些重要注意事项。有学者认为，（多元）元回归在实践中经常被不恰当地使用和解释，从而导致结果效度偏低（Higgins and Thompson, 2004）。下面我们来讨论在拟合多元元回归模型时必须牢记的一些问题。

8.3.2.1 过拟合：在并不存在信号的地方看见信号

为了更好地理解（多元）元回归模型的风险，我们需要先理解 过拟合（overfitting）的概念。过拟合是指我们建立的统计模型对数据拟合得过于紧密。本质上，这意味着模型可以很好地预测手头的数据，却很难预测未来的数据。

当模型错误地把数据中的某些变异视为真实“信号”，而实际上它捕捉到的只是随机噪声时，就会发生这种情况（Iniesta, Stahl, and McGuffin, 2016）。结果是，模型会产生 假阳性 结果，也就是在实际上并不存在关系的地方“看见”关系。

图 8.3：过拟合模型与稳健拟合模型的预测对比。

在模型拟合中，回归依赖于优化技术，例如普通最小二乘法或极大似然估计。正如我们已经学到的，元回归使用的是加权版本的普通最小二乘法（见第 8.1.3 节），因此也不例外。

然而，这种“贪婪式”的优化会使回归方法容易发生过拟合（Gigerenzer, 2004）。不幸的是，从常规回归转向元回归后，构建不稳健模型的风险反而更高，原因有几个（Higgins and Thompson, 2004）：

在元回归中，数据点数量通常很少，因为我们只能使用纳入研究的汇总信息。
由于 Meta 分析的目标是尽可能全面地概括所有可得证据，我们没有额外的数据可以用来“测试”回归模型在预测未见数据时表现如何。
在元回归中，我们还必须处理效应量异质性可能存在的问题。想象这样一种情况：有两项研究的效应量不同，且置信区间彼此不重叠。那么，只要某个变量在这两项研究中的取值不同，它似乎都可能成为解释效应量差异的候选因素。但很显然，这些解释中的大多数都可能只是伪相关。
元回归，尤其是多元元回归，非常容易让人对预测变量“反复试探”。我们可以不断测试不同的元回归模型，加入更多预测变量，或移除某些变量，试图解释数据中的异质性。这种做法很诱人，而且实践中也很常见，因为 Meta 分析研究者总想找到效应量不同的解释（J. Higgins et al., 2002）。然而，这种行为已被证明会大幅提高伪发现的风险，因为我们可以无限制地修改模型，直到得到一个显著模型，而这个模型很可能已经过拟合了，也就是说，它主要拟合的只是统计噪声。

为了避免在构建元回归模型时产生过高的假阳性率，人们提出了一些建议：

尽量减少要考察的预测变量数量。在多元元回归中，这对应于 简约性（parsimony）的概念：在评估模型拟合时，我们更偏好那些用更少的预测变量就能达到良好拟合的模型。Akaike 信息准则和贝叶斯信息准则等指标有助于进行这类判断，我们将在后面的实操示例中展示如何解释这些指标。
预测变量的选择应基于我们在 Meta 分析中想要回答的、预先设定的、具有科学意义的问题。至关重要的是，应当在分析报告中（见第 1.4.2 节）事先明确，元回归模型将纳入哪些预测变量（组合）。如果我们最后决定运行一个分析计划中未提及的元回归，也并非不可接受，但这时应在 Meta 分析报告中如实说明：这个模型是在看到数据之后才拟合的。
当研究数量较少时（这在实践中很可能会发生），如果我们还希望检验某个预测变量的显著性，就应使用 Knapp-Hartung 调整，以获得更稳健的估计。
我们还可以使用置换（permutation）方法，在重抽样数据中评估模型的稳健性。稍后我们会介绍这一方法的细节。

8.3.2.2 多重共线性

多重共线性（multi-collinearity）是指回归模型中的一个或多个预测变量，可以被另一个模型预测变量以较高精度预测出来（Mansfield and Helms, 1982）。这通常意味着模型中的两个或更多自变量彼此高度相关。

多重共线性的多数风险都与过拟合问题相关。高共线性会使预测变量系数估计 $\hat\beta$ 的行为变得不稳定，并且在数据稍有变化时就发生明显波动。它还会限制模型可解释变异的大小，在这里也就是限制 $R^2_*$ 的数值。

多重共线性在元回归中很常见（Berlin and Antman, 1994）。尽管多元回归能够处理较低程度的共线性，但对于高度相关的预测变量，我们仍然应当检查，并在必要时加以控制。判断是否存在多重共线性，并没有一个统一的“是 / 否”规则。

一种粗略但往往有效的方法，是在拟合模型之前先检查预测变量之间是否存在很高的相关（例如 $r \geq 0.8$）。如果存在，多重共线性可以通过两种方式缓解：1）删除其中一个几乎冗余的预测变量；2）尝试把这些预测变量合并为一个单一变量。

8.3.2.3 模型拟合方法

在构建多元元回归模型时，选择和纳入预测变量的方法有很多。下面我们讨论其中最重要的几种，并概述它们各自的优缺点：

Forced entry。在 forced entry 方法中，所有相关预测变量会被同时强制纳入回归模型。对于 R 中的大多数函数来说，这就是默认做法。虽然这通常是一种推荐策略，但仍需注意，采用 forced entry 纳入的所有预测变量，也应当基于事先设定、由理论驱动的决定。
Hierarchical。分层多元回归（hierarchical multiple regression）是指根据明确的科学理由，按步骤把预测变量依次加入模型。首先，按照其重要性顺序纳入那些在既有研究中已经与效应量差异相关的预测变量；之后，再加入新的预测变量，以探索它们是否能够解释尚未被已知预测变量捕捉到的异质性。
Step-wise。逐步进入（step-wise entry）意味着把变量 / 预测变量一个接一个地加入模型。乍看之下，这与分层回归有些相似，但二者的关键差别在于：逐步回归是根据 统计准则 来选择预测变量的。在所谓 前向选择（forward selection）中，首先纳入的是最能解释数据变异的变量。之后对剩余变量重复这一过程，每次都选择最能解释剩余未解释变异的变量。还有一种方法叫 后向选择（backward selection），即一开始把所有变量都纳入模型，然后依据预设统计准则逐步删去变量。大量文献都不鼓励使用逐步法（Chatfield, 1995；Whittingham et al., 2006）。回想上文讨论的多元回归常见陷阱就会发现，这些方法很容易产生过拟合模型和伪发现。尽管如此，逐步法在实践中仍被频繁使用，因此了解这类方法的存在仍然很重要。不过，如果我们自己使用逐步法，最好仅将其作为探索性工具，并始终牢记其局限。
Multi-model inference。多模型推断与逐步法不同，它并不试图逐步构建出一个解释方差最多的“最佳”模型。相反，这种技术会对 所有可能的 预测变量组合都进行建模。这意味着会创建多个不同的元回归模型，并对它们逐一评估。这样，我们就能全面考察所有可能的预测变量组合及其表现。一个常见发现是：很多不同的模型设定都能产生良好的拟合。此时，可以把所有已拟合模型中的预测变量系数加总起来，以推断哪些变量总体上更重要。

8.3.3 在 R 中进行多元元回归

读完前面的理论内容之后，现在终于可以开始在 R 中拟合第一个多元元回归模型了。下面的示例将是本章第一次不再使用 {meta} 软件包，而改用 {metafor}（Viechtbauer, 2010）。这个软件包为 Meta 分析提供了大量高级功能，也有非常好的文档。开始之前，请确保已经安装并加载了 {metafor}。

library(metafor)

在下面的实操示例中，我们将使用 MVRegressionData 数据集。它是一个为了教学演示而模拟出的“玩具”数据集。

“MVRegressionData” 数据集

MVRegressionData 数据集直接包含在 {dmetar} 软件包中。如果你已经安装了 {dmetar} 并从库中加载它，那么运行 data(MVRegressionData) 就会自动把该数据集保存到你的 R 环境中。这样数据集就可以直接使用了。如果你没有安装 {dmetar}，也可以从网络上下载其 .rda 文件，保存到工作目录中，然后在 R Studio 窗口中点击它完成导入。

首先，我们来看一下这个数据框的结构：

library(tidyverse)
library(dmetar)
data(MVRegressionData)

glimpse(MVRegressionData)

## Rows: 36
## Columns: 6
## $ yi         <dbl> 0.09437543, 0.09981923, 0.16931607, 0.17511107, 0.27301641, ...
## $ sei        <dbl> 0.1959031, 0.1918510, 0.1193179, 0.1161592, 0.1646946, ...
## $ reputation <dbl> -11, 0, -11, 4, -10, -9, -8, -8, -8, 0, -5, -5, -4, -4, -3, ...
## $ quality    <dbl> 6, 9, 5, 9, 2, 10, 6, 3, 10, 3, 1, 5, 10, 2, 1, 2, 4, 1, 8, ...
## $ pubyear    <dbl> -0.85475360, -0.75277184, -0.66048349, -0.56304843, -0.4308...
## $ continent  <fct> 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, ...

我们可以看到，这个数据集中共有六个变量。yi 和 sei 两列分别存储某项研究的效应量和标准误。它们对应于我们之前使用过的 TE 和 seTE 两列。之所以这样命名，是因为这正是 {metafor} 的标准记法：yi 表示我们在（元）回归中想要预测的观测效应量 $y_i$，而 sei 表示研究 $i$ 的标准误 $SE_i$。

其余四个变量则是元回归中要使用的预测变量。首先是 reputation，即研究发表期刊的（均值中心化后的）影响因子。影响因子量化的是期刊文章被引用的频率，我们将其作为期刊声望的代理指标。

其他变量包括：quality，即 0 到 10 分的研究质量评分；pubyear，即经过中心化和标准化处理的发表年份；以及 continent，即研究实施所在的大洲。这些变量中，除了 continent 之外都是连续变量；continent 是一个具有两个水平的分类变量：Europe 和 North America。

8.3.3.1 检查多重共线性

正如前面提到的，我们需要检查预测变量之间是否存在多重共线性，以确保元回归系数估计是稳健的。一个快速的方法是计算所有连续变量之间的 相关矩阵。这可以通过 cor 函数来完成：

MVRegressionData[,c("reputation", "quality", "pubyear")] %>% cor()

##            reputation    quality    pubyear
## reputation  1.0000000  0.3015694  0.3346594
## quality     0.3015694  1.0000000 -0.1551123
## pubyear     0.3346594 -0.1551123  1.0000000

{PerformanceAnalytics} 软件包（Peterson and Carl, 2020）中有一个名为 chart.Correlation 的函数，可以把这个相关矩阵可视化。我们首先需要安装 PerformanceAnalytics 软件包，然后运行以下代码：

library(PerformanceAnalytics)

MVRegressionData[,c("reputation", "quality", "pubyear")] %>% 
  chart.Correlation()

（原页面此处展示了相关矩阵图。）

从结果来看，这些变量确实存在相关，但大概还没有高到需要把其中某个变量剔除的程度。

8.3.3.2 拟合多元元回归模型

现在，我们可以使用 {metafor} 来拟合第一个元回归模型了。此前，我们想探究的是：高期刊声誉是否会预测更高的效应量，还是这只是因为高声誉期刊中的研究质量通常更高，从而形成的一种假象。

假设我们已经非常清楚地知道，例如来自既有研究，研究质量确实能够预测效应量。那么在这种情况下，做一个分层回归就是合理的：我们先纳入已知预测变量 quality，再检验 reputation 是否能在此基础上进一步解释异质性。如果结果成立，我们就可以说，即使控制了高声誉期刊更可能发表高质量研究这一事实，期刊声誉仍然与更高的效应相关。

为此，我们要使用 {metafor} 中的 rma 函数。这个函数原本用于运行随机效应 Meta 分析；当加入调节变量之后，它就扩展为混合效应元回归模型。rma 函数可以接受大量参数，我们可以通过在 R 控制台中运行 ?rma 来查看。通常情况下，我们只需要指定其中几个：

yi。数据框中存储每项研究效应量的那一列。
sei。数据框中存储每项研究效应量标准误的那一列。
data。包含全部 Meta 分析数据的数据框名称。
method。我们想使用的 $\tau^2$ 估计量。这个参数可用的代码与 {meta} 中一致（例如 "REML" 表示限制性极大似然）。如果后面想比较不同的元回归模型，使用 "ML" 往往更方便。
mods。这个参数定义元回归模型。先用 ~ 指定模型，再把想纳入的预测变量用 + 分隔开（例如 variable1 + variable2）。两个变量之间的交互作用用星号表示（例如 variable1 * variable2）。
test。用于检验回归系数的方法，可以选择 "z"（默认）或 "knha"（Knapp-Hartung 方法）。

首先，我们只把 quality 作为预测变量来做一次元回归。结果保存到对象 m.qual 中，然后查看输出。

m.qual <- rma(yi = yi,
              sei = sei,
              data = MVRegressionData,
              method = "ML",
              mods = ~ quality,
              test = "knha")

m.qual

## Mixed-Effects Model (k = 36; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.066 (SE = 0.023)
## tau (square root of estimated tau^2 value):             0.2583
## I^2 (residual heterogeneity / unaccounted variability): 60.04%
## H^2 (unaccounted variability / sampling variability):   2.50
## R^2 (amount of heterogeneity accounted for):            7.37%
## 
## Test for Residual Heterogeneity:
## QE(df = 34) = 88.6130, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 34) = 3.5330, p-val = 0.0688
## 
## Model Results:
## 
##          estimate      se    tval    pval    ci.lb   ci.ub 
## intrcpt    0.3429  0.1354  2.5318  0.0161   0.0677  0.6181  * 
## quality    0.0356  0.0189  1.8796  0.0688  -0.0029  0.0740  . 
## 
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

在输出中，我们可以在 Model Results 部分查看预测变量 quality 的结果。可以看到，这个回归权重并不显著（$p = 0.069$），不过在趋势水平上是显著的（$p < 0.1$）。总体而言，该模型解释了 $R^2_* = 7.37\%$ 的异质性。

接下来，我们再把 reputation 加入为预测变量。只需在 mods 中加入 + reputation，并把输出保存为 m.qual.rep。

m.qual.rep <- rma(yi = yi, 
                  sei = sei, 
                  data = MVRegressionData, 
                  method = "ML", 
                  mods = ~ quality + reputation, 
                  test = "knha")

m.qual.rep

## Mixed-Effects Model (k = 36; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0238 (SE = 0.01)
## tau (square root of estimated tau^2 value):             0.1543
## I^2 (residual heterogeneity / unaccounted variability): 34.62%
## H^2 (unaccounted variability / sampling variability):   1.53
## R^2 (amount of heterogeneity accounted for):            66.95%
## 
## Test for Residual Heterogeneity:
## QE(df = 33) = 58.3042, p-val = 0.0042
## 
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 33) = 12.2476, p-val = 0.0001
## 
## Model Results:
## 
##             estimate      se    tval    pval    ci.lb   ci.ub 
## intrcpt       0.5005  0.1090  4.5927  <.0001   0.2788  0.7222  *** 
## quality       0.0110  0.0151  0.7312  0.4698  -0.0197  0.0417      
## reputation    0.0343  0.0075  4.5435  <.0001   0.0189  0.0496  *** 
## 
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

现在，在 Model Results 中出现了新的一行，显示预测变量 reputation 的结果。模型估计该回归权重为 0.034，而且高度显著（$p < 0.001$）。

我们还可以看到，整个元回归模型解释了显著的一部分异质性，具体来说是 $R^2_* = 66.95\%$。这意味着，即使在控制研究质量之后，期刊声誉仍然与更高的效应量相关。

但是，我们的第二个模型是否真的比第一个模型拟合得更好呢？为了评估这一点，我们可以使用 anova 函数，并把想要比较的两个模型都传进去。需要注意的是，这之所以可行，是因为我们这两个混合效应模型都使用了极大似然（"ML"）而不是限制性极大似然（"REML"）进行拟合。

anova(m.qual, m.qual.rep)

##         df   AIC   BIC  AICc logLik   LRT   pval    QE tau^2    R^2 
## Full     4 19.86 26.19 21.15  -5.93              58.30  0.03          
## Reduced  3 36.98 41.73 37.73 -15.49 19.11 <.0001 88.61  0.06 48.32% 

这个函数执行的是模型比较，并提供多个统计量，帮助我们评估 m.qual.rep 是否比 m.qual 拟合更好。这里，我们把完整模型 m.qual.rep（包含 quality 和 reputation）与简化模型 m.qual（只包含 quality）进行了比较。

anova 函数执行的是 似然比检验（likelihood ratio test），结果可以在 LRT 列中看到。这个检验高度显著（$\chi^2_1 = 19.11, p < 0.001$），说明完整模型的拟合确实更好。

另一个重要统计量是 AICc 列。它表示经过小样本修正后的 Akaike 信息准则（AIC）。正如前面提到的，AICc 会对包含更多预测变量的复杂模型施加惩罚，从而避免过拟合。

需要注意的是，AIC 值越低，模型表现越好。在这里，尽管完整模型包含更多参数，但它的 AICc 值（21.15）仍然低于简化模型（37.73）。这些结果都表明，我们的多元回归模型确实对数据有更好的拟合。

8.3.3.3 建模交互作用

假设我们现在想为两个额外预测变量 pubyear（发表年份）和 continent 建立一个交互模型。我们假设：发表年份与效应量之间的关系，在欧洲研究与北美研究中并不相同。要在 rma 函数中表示这一假设，我们需要在 mods 参数里用 * 将预测变量连接起来。由于这一次我们不打算用 anova 直接比较模型，因此这里改用 "REML"（限制性极大似然）来估计 $\tau^2$。

为了便于解释结果，我们在运行模型之前，先给 MVRegressionData 中的 continent 变量加上因子标签。

# Add factor labels to 'continent'
# 0 = Europe
# 1 = North America
levels(MVRegressionData$continent) = c("Europe", "North America")

# Fit the meta-regression model
m.qual.rep.int <- rma(yi = yi, 
                      sei = sei, 
                      data = MVRegressionData, 
                      method = "REML", 
                      mods = ~ pubyear * continent, 
                      test = "knha")

m.qual.rep.int

## Mixed-Effects Model (k = 36; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0 (SE = 0.01)
## tau (square root of estimated tau^2 value):             0
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability):   1.00
## R^2 (amount of heterogeneity accounted for):            100.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 32) = 24.8408, p-val = 0.8124
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 32) = 28.7778, p-val < .0001
## 
## Model Results:
## 
##                       estimate    se  tval    pval  ci.lb ci.ub 
## intrcpt                    0.38  0.04  9.24  <.0001   0.30  0.47  *** 
## pubyear                    0.16  0.08  2.01  0.0520  -0.00  0.33    . 
## continentNorth America     0.39  0.06  6.05  <.0001   0.26  0.53  *** 
## pubyear:continent          0.63  0.12  4.97  <.0001   0.37  0.89  *** 
##    North America
## [...]

最后一行 pubyear:continentNorth America 就是交互项的系数。注意，{metafor} 会自动把交互项以及 pubyear 和 continent 这两个“普通”的低阶预测变量一并包含进模型中，这也是规范做法。

还要注意，由于 continent 是一个因子变量，rma 检测到它是一个虚拟编码的预测变量，因此把 “Europe” 作为 $D_g = 0$ 的基准组，并用它来与 “North America” 进行比较。我们看到，交互项系数为正（0.63），而且高度显著（$p < 0.001$）。

这表明，近年来效应量确实在增大，而且这一趋势在北美开展的研究中更强。我们还看到，该模型解释了 $R^2_* = 100\%$ 的异质性。

这是因为我们的数据是为演示目的而模拟出来的。在真实研究中，你几乎不可能把数据中的异质性全部解释掉。事实上，如果在真实数据中看到这样的结果，反而应该提高警惕，因为这可能意味着模型已经过拟合。

8.3.3.4 置换检验

置换（permutation）是一种数学操作：从一个包含若干数字或对象的集合中，反复抽取元素并把它们重新排成序列。如果我们原本已经有一个有序数字集合，那么这就相当于重新排列、也就是打乱数据顺序的过程。

例如，设集合 $S$ 包含三个数字：$S = \{1,2,3\}$。这个集合的一种置换是 $(2,1,3)$，另一种则是 $(3,2,1)$。我们看到，置换后的结果仍然包含原来的三个数字，只是顺序不同而已。

置换还可以用来进行 置换检验（permutation test），它是一种特殊的重抽样方法。广义上说，重抽样方法通过向统计模型提供来自同一来源或生成过程、但略有差异的数据样本，来检验模型的稳健性（Good, 2013，第 3.1 章）。这有助于我们判断：模型中的系数是否真的捕捉到了数据背后的真实模式，还是仅仅因为过拟合而把统计噪声误当成了规律。

置换检验并不要求我们额外准备一个“测试”数据集，用来评估元回归在预测未见效应量时的表现。正因如此，以及其他一些原因，人们建议在评估元回归模型稳健性时使用置换检验（Higgins and Thompson, 2004）。

这里我们不深入讨论元回归置换检验的具体执行细节。最关键的一点是：我们基于原始数据集所有可能的、或大量随机抽取的置换结果所得到的检验统计量，重新计算模型的 $p$ 值。

这里最重要的指标是：在置换后的数据中，我们得到的检验统计量 有多频繁 会 大于或等于 原始检验统计量。例如，如果在 1000 个置换数据集中，有 50 个的检验统计量大于或等于原始统计量，那么得到的 $p$ 值就是 $p = 0.05$。

如果要对元回归模型执行置换检验，我们可以使用 {metafor} 内置的 permutest 函数。举例来说，我们重新计算前面拟合过的 m.qual.rep 模型结果。我们只需把 rma 对象传给 permutest 函数即可。需要注意的是，置换检验的计算代价较高，尤其是在数据集较大时，因此运行可能需要一些时间。

permutest(m.qual.rep)

## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 33) = 12.7844, p-val* = 0.0010
## 
## Model Results:
## 
##             estimate      se    tval   pval*    ci.lb   ci.ub 
## intrcpt       0.4964  0.1096  4.5316  0.2240   0.2736  0.7193      
## quality       0.0130  0.0152  0.8531  0.3640  -0.0179  0.0438      
## reputation    0.0350  0.0076  4.5964  0.0010   0.0195  0.0505  *** 
## 
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我们再次看到了熟悉的输出，其中包含所有预测变量的结果。观察 pval* 这一列，可以看到 reputation 预测变量的 $p$ 值从 $p < 0.001$ 变成了 $p_* = 0.001$。不过，这一结果依然高度显著，说明该预测变量的效应是稳健的。

人们已经建议，在报告元回归模型结果之前，最好总是先做这种置换检验（Higgins and Thompson, 2004）。

小数据集中的置换检验

请注意，当模型中纳入的研究数量 $K$ 很小时，常规使用的统计显著性阈值（即 $p < 0.05$）实际上可能无法达到。

对于元回归模型，使用 permutest 的置换检验只有在 $K > 4$ 时才有可能达到统计显著性（Viechtbauer et al., 2015）。

8.3.3.5 多模型推断

前面我们已经提到，还可以通过一种叫做 多模型推断（multi-model inference）的程序，对所有可能的预测变量组合都进行建模。这样可以帮助我们检验：哪一种预测变量组合的拟合最好，以及总体上哪些预测变量最为重要。要实现多模型推断，我们可以使用 multimodel.inference 函数。

注：如果你想进一步了解这个主题，也可以参考 Wolfgang Viechtbauer 作为 {metafor} 文档一部分撰写的一篇说明性 vignette。

“multimodel.inference” 函数

multimodel.inference 函数包含在 {dmetar} 软件包中。一旦 {dmetar} 已经安装并在你的计算机上加载完成，这个函数就可以直接使用了。如果你没有安装 {dmetar}，可以按以下步骤操作：

在线访问这个函数的源代码。
将完整源代码复制粘贴到 R 控制台（R Studio 左下角窗格）中，让 R “学会”这个函数，然后按下 Enter。
确保 {metafor}、{ggplot2} 和 {MuMIn} 软件包都已经安装并加载。

在这个函数中，需要指定以下参数：

TE。每项研究的效应量。必须以字符串形式提供数据集中的效应量列名（例如 TE = "effectsize"）。
seTE。效应量的标准误。必须提供标准误列名，同样需要用引号括起来（例如 seTE = "se"）。
data。一个包含效应量、标准误以及元回归预测变量的数据框。
predictors。一个由字符串拼接而成的向量，用来指定进行多模型推断时要使用的预测变量。预测变量名称必须与传给 data 的数据框列名完全一致。
method。用于合并效应量的 Meta 分析模型。"FE" 表示固定效应模型；随机效应模型还可选择 "DL"、"SJ"、"ML" 或 "REML" 等。如果使用 "FE"，则 test 参数会自动被设置为 "z"，因为 Knapp-Hartung 方法不适用于固定效应模型。默认值是 "REML"。
test。用于计算检验统计量和置信区间的方法。默认值是 "knha"，即使用 Knapp-Hartung 调整；若要计算常规 Wald 型检验，则把该参数设为 "z"。
eval.criterion。用于评估拟合模型的准则。可以是 "AICc"（默认值，即小样本修正后的 Akaike 信息准则）、"AIC"（Akaike 信息准则）或 "BIC"（贝叶斯信息准则）。
interaction。若设为 FALSE（默认），则不考虑预测变量之间的交互作用；若设为 TRUE，则会把所有交互项都纳入建模。

现在，我们使用 MVRegressionData 中的全部预测变量来进行多模型推断，但不包括交互项。请注意，运行 multimodel.inference 可能需要一些时间，尤其是在预测变量数量较多时。

multimodel.inference(TE = "yi", 
                     seTE = "sei",
                     data = MVRegressionData,
                     predictors = c("pubyear", "quality", 
                                    "reputation", "continent"),
                     interaction = FALSE)

## Multimodel Inference: Final Results
## --------------------------
##  - Number of fitted models: 16
##  - Full formula: ~ pubyear + quality + reputation + continent
##  - Coefficient significance test: knha
##  - Interactions modeled: no
##  - Evaluation criterion: AICc 
## 
## 
## Best 5 Models
## --------------------------
## [...]
##    (Intrc) cntnn  pubyr   qulty   rpttn df logLik AICc delta weight
## 12       +     + 0.3533         0.02160  5  2.981  6.0  0.00  0.536
## 16       +     + 0.4028 0.02210 0.01754  6  4.071  6.8  0.72  0.375
## 8        +     + 0.4948 0.03574          5  0.646 10.7  4.67  0.052
## 11       +       0.2957         0.02725  4 -1.750 12.8  6.75  0.018
## 15       +       0.3547 0.02666 0.02296  5 -0.395 12.8  6.75  0.018
## Models ranked by AICc(x) 
## 
## 
## Multimodel Inference Coefficients
## --------------------------
##                          Estimate  Std. Error   z value  Pr(>|z|)
## intrcpt                0.38614661 0.106983583 3.6094006 0.0003069
## continentNorth America 0.24743836 0.083113174 2.9771256 0.0029096
## pubyear                0.37816796 0.083045572 4.5537402 0.0000053
## reputation             0.01899347 0.007420427 2.5596198 0.0104787
## quality                0.01060060 0.014321158 0.7402055 0.4591753
## 
## 
## Predictor Importance
## --------------------------
##        model importance
## 1    pubyear  0.9988339
## 2  continent  0.9621839
## 3 reputation  0.9428750
## 4    quality  0.4432826

（原页面此处展示了模型平均后的预测变量重要性图。）

这里有很多信息，我们一步一步来看。

Multimodel Inference: Final Results。这部分输出提供了拟合模型的总体信息。我们可以看到，总共拟合了 $2^4 = 16$ 个可能模型。还可以看到，该函数使用修正后的 AIC（aicc）来比较模型。
Best 5 Models。这里展示了 AICc 最低的五个模型，按从低到高排序。表格的列表示预测变量，行表示模型。数字（即权重）或 + 号（对于分类预测变量）表示某个预测变量 / 交互项被纳入该模型，空白单元格则表示该预测变量被省略。我们看到，TE ~ 1 + continent + pubyear + reputation 的拟合最好（AICc = 6.0）。不过，其他几种预测变量组合与它的差距也非常小，因此很难断言哪一个才是真正“最佳”的模型。但前五个模型全部包含 pubyear，这提示该变量可能尤其重要。
Multimodel Inference Coefficients。这里给出了所有预测变量的系数，这些系数是把各个变量在所有包含它的模型中的结果汇总后得到的。我们看到，pubyear 的系数估计最大（$\hat\beta = 0.378$），这与我们前面的发现一致。近似置信区间可以通过在 Estimate 上减去和加上 Std.Error × 1.96 得到。
模型平均后的预测变量重要性图。在图中，展示了每个预测变量在所有模型中的平均重要性。我们再次看到，pubyear 是最重要的预测变量，其次依次是 reputation、continent 和 quality。

多模型推断的局限

这个例子说明，多模型推断可以是一种有用的方法，用来全面考察哪些预测变量对解释效应量差异最为重要。

尽管它避免了逐步回归方法的一些问题，但仍要注意：这种方法应当被看作 探索性 的，适用于我们对某一研究领域中预测变量与效应量之间关系缺乏先验知识的时候。

如果你决定根据多模型推断的结果来建立一个元回归模型，那么务必在报告中明确说明这一点。因为这样的模型并不是建立在 先验假设 基础上的，而是根据样本中的统计性质构建出来的。

8.4 问答

检验一下你的掌握程度！

原始研究中使用的常规回归分析，与元回归之间有什么区别？
亚组分析与元回归关系密切。那么，元回归公式应如何针对亚组数据进行调整？
在元回归中，使用哪种方法来赋予不同研究不同的权重？
一个能够很好拟合我们数据的元回归模型应具有哪些特征？可以用哪个指标来检验？
当我们用元回归技术来计算亚组分析时，在各个亚组中是假定 $\tau^2$ 分别不同，还是共享同一个值？
（多元）元回归有哪些局限和陷阱？
请说出两种可以提高（多元）元回归模型稳健性的方法，并说明它们为什么有帮助。

这些问题的答案列在本书末尾的附录 A 中。

8.5 小结

在元回归中，我们把常规回归技术应用到研究层面的数据上。亚组分析可以看作带有分类预测变量、并假定共享 $\tau^2$ 估计值的元回归的一种特殊情况。
元回归模型的目标，是解释数据中真实效应量差异的全部或部分来源，也就是研究间异质性方差 $\tau^2$。当模型拟合良好时，真实效应相对于回归线的偏离，应小于它们相对于合并效应的初始偏离。在这种情况下，未被解释的、也就是残余异质性就会较小。这一点由 $R^2_*$ 指标来刻画，它告诉我们模型解释了多少比例的异质性变异。
在 多元元回归 中，同一个元回归模型会同时使用两个或更多预测变量。我们还可以通过加入交互项，检验某个变量的预测作用是否会随着另一个变量的不同取值而改变。
尽管（多元）元回归非常灵活，但它并非没有局限。多元元回归很容易导致 过拟合，也就是模型拟合的不是真实关系，而是随机噪声。预测变量之间的多重共线性也可能威胁模型的有效性。
有多种方法可以帮助我们确保元回归模型的稳健性。例如，我们可以只拟合那些基于 预先设定 理论依据的模型，或者使用置换检验。多模型推断也可以作为一种探索性方法，帮助我们识别潜在的重要预测变量，并据此提出未来研究中需要检验的假设。

Chaofan Chen