Bayesian meta analysis in Brms

元分析

导航

• 元分析
• 会话信息
• 评论

如果你还没有被那些测量误差模型和缺失数据模型彻底震撼，那我们就继续深入。事实证明，元分析往往只是多层测量误差模型的一种特殊形式。因此，你同样可以使用 brms::brm() 来拟合贝叶斯元分析。

在继续之前，我需要说明，本节在很大程度上受到了 Matti Vourre 那篇优秀博文 Bayesian meta-analysis with R, Stan & brms 的影响。由于 McElreath 的教材并未直接讨论元分析，我们还需要从 Gelman、Carlin、Stern、Dunson、Vehtari 和 Rubin 的 Bayesian data analysis, Third edition 中借鉴一些内容。下面就让 Gelman 及其合作者为我们引入这一主题：

关于元分析的讨论有时会对分析中所关注的估计目标表述得不够精确，尤其是在主要关注点是检验拟合并研究中“无效应”原假设时更是如此。我们的关注点在于估计有意义的参数；在这一目标下，前提是接受一个总括性的假设，即这些研究在某种宽泛意义上是可比的。第一种可能性是，我们将这些研究视为彼此完全相同的重复，亦即认为所有研究中的个体都可被看作来自同一总体的独立样本，并且具有相同的结局测量等。第二种可能性是，这些研究之间差异巨大，以至于任意一项研究的结果都无法为其他任何研究的结果提供信息。第三种、也更一般的可能性是，我们将这些研究视为可交换的，但不必然是相同的，也不必然是完全无关的；换言之，我们允许研究之间存在差异，但这种差异在 a priori 意义上并不预期会以可预测的方式偏向某项研究而非另一项研究……第三种可能性代表着上述两个极端之间的一个连续体，而正是这种可交换模型（其总体分布由未知超参数加以刻画）构成了我们贝叶斯分析的基础……

元分析，或者更一般地说层级结构问题的第一个潜在估计目标，是效应量分布的均值，因为它代表了所有那些可被视为与已观测研究可交换的研究中的总体“平均”效应。其他可能的估计目标还包括任一已观测研究中的效应量，以及另一项可比的（即可交换的）未观测研究中的效应量。（第 125—126 页，原文 emphasis）

贝叶斯元分析的基本形式如下：

其中，$y_i$ 表示单个研究 $i$ 的效应量点估计，并假定它来自一个以 $\theta_i$ 为中心的正态分布。元分析中的数据通常是各项原始研究的统计汇总量。本章一个非常明确的启示是：这些估计值本身也带有误差，而这些误差应当在元分析模型中得到充分表达。我们正是这样做的。研究 $i$ 的标准误被记作 $\sigma_i$，它同时也可视为生成该点估计的正态分布的标准差。请注意，这里我们并不估计 $\sigma_i$。这些数值直接取自原始研究。

在这一模型基础上，我们进一步假定研究 $i$ 本身只是相关研究总体中的一次抽样，而这些研究各自都具有自身的效应量。因此，我们假定 $\theta_i$ 本身也服从如下分布：

其中，$\mu$ 是元分析效应（即总体均值），而 $\tau$ 是围绕该均值的变异程度，你也可以将其理解为 $\sigma_\tau$。

由于正文中没有提供元分析的例子，我们只好从别处获取数据。这里我们聚焦于 Gershoff 和 Grogan-Kaylor 的论文 Spanking and child outcomes: Old controversies and new meta-analyses。在其引言中，我们读到：

在世界范围内，大多数儿童（80%）都曾遭到父母打屁股或其他形式的身体惩罚。关于父母是否应当通过打屁股来纠正儿童不当行为的问题，在美国乃至全球范围内，处于伦理、宗教和人权等多种观点争论的交汇点上。已有数百项研究考察了父母使用打屁股或身体惩罚与儿童行为、情绪、认知以及身体结局之间的关联，因此，打屁股成为育儿领域中研究最为充分的主题之一。我们从这数百项研究中究竟学到了什么？（第 453 页）

我们的目标是通过回答这一问题中的一部分来学习贝叶斯元分析。我已经将 Gershoff 和 Grogan-Kaylor 文中的数值直接转录下来，并将其保存为名为 spank.xlsx 的文件。 你可以在 this project’s GitHub repository 中找到这些数据。下面让我们将其载入，并使用 glimpse() 查看其结构。

spank <- readxl::read_excel("spank.xlsx")
glimpse(spank)

在这篇论文中，关注的效应量是 Cohen’s $d$，其计算公式为

其中

为了使上述关于 $d$ 的公式在我们的例子中更为清晰，我们可以将其改写为

McElreath 在其教材中并未真正聚焦于效应量。如果你需要复习相关内容，可以参考 Kelley 和 Preacher 的 On effect size。简言之，Cohen’s $d$ 是两个组之间标准化后的均值差。

因此，如果你回头查看 glimpse(spank) 的结果，就会注意到其中的 d 列；它确实是一个 Cohen’s $d$ 效应量向量。最后两列 ll 和 ul 分别是对应 95% 频率学置信区间的下限和上限。但我们并不想要这些 d 值的置信区间；我们想要的是它们的标准误。幸运的是，我们可以用下式计算：

其代码如下。

spank <- spank %>%
  mutate(se = (ul - ll) / 3.92)

glimpse(spank)

现在数据已经准备就绪，我们可以用以下公式表达第一个贝叶斯元分析模型：

当然，最后两行给出的就是我们的先验分布。在心理学中，Cohen’s $d$ 的绝对值大于 $\pm 1$ 的情况相当少见。因此，在缺乏更具体领域知识的情况下，而我确实也没有这样的知识，$\operatorname{Normal}(0, 1)$ 看起来是一个合理的起点。并且，就像 McElreath 将 $\operatorname{HalfCauchy}(0, 1)$ 用作组层标准差的默认先验一样，在这里将其用于 我们元分析中的 $\tau$ 参数也是合理的。

下面是第一个模型的代码。

b14.5 <- brm(
  data = spank,
  family = gaussian,
  d | se(se) ~ 1 + (1 | study),
  prior = c(prior(normal(0, 1), class = Intercept),
            prior(cauchy(0, 1), class = sd)),
  iter = 4000, warmup = 1000, cores = 4, chains = 4,
  seed = 14,
  file = "fits/b14.05")

你可能会注意到，我们的 se(se) 函数中没有包含 sigma 参数。如果你还记得前文关于测量误差的部分，在测量误差模型里我们曾指定 sigma = T。而 brms 的默认设置是，在 se() 内部 sigma = FALSE。因此，在这里我们并没有像使用各研究原始数据进行分析时那样再去估计 sigma。希望这很好理解。每个研究 $i$ 的 d 值所对应的不确定性，已经通过数据中的 se 编码进来了。

这也引出了另一个问题。我们通常是在数据汇总层面上进行元分析。在我的领域中，也许在你的领域中也是如此，这在一定程度上源于一个历史性的偶然：研究者过去并不常规地公开其原始数据。因此，效应量汇总往往就是我们通常所能得到的最好材料。然而，时代正在改变。如果你的元分析所纳入研究的全部原始数据都可获得，那么你完全可以直接拟合一个多层模型，使数据嵌套在研究之中。甚至，你还可以通过采用分布式建模方法让研究之间在 $\sigma$ 上也存在差异，并指定诸如 sigma ~ 0 + study，甚至 sigma ~ 1 + (1 | study) 这样的模型。

不过，技术性的讨论就先到这里。让我们看看模型结果。

print(b14.5)

因此，在这个简单的贝叶斯元分析中，我们得到的总体 Cohen’s $d$ 大约为 r round(fixef(b14.5)[1], 2)。而我们对 $\tau$ 的估计，即 r round(posterior_summary(b14.5)[2, 1], 2)，表明研究之间存在相当大的变异。你可能会问：这些 Cohen’s $d$ 到底测量的是什么？ 我们已经将这一点编码在 spank 数据中的 outcome 向量里。

spank %>%
  distinct(outcome) %>%
  knitr::kable()

这里有几点值得注意。首先，除了 Adult support for physical punishment 可能是一个例外之外，其余所有结局都是负向的。对于诸如 Child aggression 和 Adult mental health problems 之类的指标，我们更偏好其取值较低的情形。其次，按照数据的编码方式，更大的效应量被解释为与儿童曾被打屁股相关的更负面的结局。也就是说，我们的分析表明，打屁股与儿童较差的结局相关。还有一点可能并不那么显而易见：尽管数据中共有 111 个案例，但实际上只有 76 项不同的研究。

spank %>%
  distinct(study) %>%
  count()

换言之，有些研究包含多个结局。为了更好地容纳 study 层面和 outcome 层面的方差，我们来拟合一个交叉分类的贝叶斯元分析模型，这与前文关于交叉分类黑猩猩模型的小节思路相似。

b14.6 <- brm(
  data = spank,
  family = gaussian,
  d | se(se) ~ 1 + (1 | study) + (1 | outcome),
  prior = c(prior(normal(0, 1), class = Intercept),
            prior(cauchy(0, 1), class = sd)),
  iter = 4000, warmup = 1000, cores = 4, chains = 4,
  seed = 14,
  file = "fits/b14.06")

print(b14.6)

现在我们有了两个 $\tau$ 参数。我们可以把它们画出来，以便感受方差大致分布在什么位置。

# we'll want this to label the plot
label <- tibble(
  tau   = c(.12, .3),
  y     = c(15, 10),
  label = c("sigma['outcome']", "sigma['study']"))

# wrangle
as_draws_df(b14.6) %>%
  select(starts_with("sd")) %>%
  gather(key, tau) %>%
  mutate(key = str_remove(key, "sd_") %>% str_remove(., "__Intercept")) %>%

  # plot
  ggplot(aes(x = tau)) +
  geom_density(aes(fill = key),
               color = "transparent") +
  geom_text(data = label,
            aes(y = y, label = label, color = label),
            size = 5, parse = T) +
  scale_fill_viridis_d(NULL, option = "B", begin = .5) +
  scale_color_viridis_d(NULL, option = "B", begin = .5) +
  scale_x_continuous(expression(tau), limits = c(0, NA)) +
  scale_y_continuous(NULL, breaks = NULL)

因此，到目前为止，最主要的结论是：研究之间的变异大于不同结局之间的变异。不过，我仍然希望进一步把握各个结局本身的情况。这里我们将使用 tidybayes::geom_halfeye() 来帮助绘制我们版本的 forest plot，并使用 tidybayes::spread_draws() 来完成初始数据整理。

library(tidybayes)

b14.6 %>%
  spread_draws(b_Intercept, r_outcome[outcome,]) %>%
  # add the grand mean to the group-specific deviations
  mutate(mu = b_Intercept + r_outcome) %>%
  ungroup() %>%
  mutate(outcome = str_replace_all(outcome, "[.]", " ")) %>%

  # plot
  ggplot(aes(x = mu, y = reorder(outcome, mu), fill = reorder(outcome, mu))) +
  geom_vline(xintercept = fixef(b14.6)[1, 1], color = "grey33", linewidth = 1) +
  geom_vline(xintercept = fixef(b14.6)[1, 3:4], color = "grey33", linetype = 2) +
  stat_halfeye(.width = .95, size = 2/3, color = "white") +
  scale_fill_viridis_d(option = "B", begin = .2) +
  labs(x = expression("Cohen's "*italic(d)),
       y = NULL) +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank())

图中背景里的白色实线和虚线分别标示了总体均值（即元分析效应）及其 95% 区间。不过，无论如何，不同结局之间的变异并不算大。让我们在模型上再向前推进一步。回到 Gelman 及其合作者的论述，我们读到：

在假定可交换性时，我们实际上假定不存在能够构成更复杂模型基础的重要协变量，而这一假定在元分析中被广泛采用，尽管它或许并不完全恰当。如果在一项元分析中，除了数据 $(n, y)$ 之外，还有其他信息可以用来区分这 $J$ 项研究，从而使可交换模型不再适用，那么该怎么办？在这种情形下，我们可以扩展模型框架，使其在观测数据和协变量的条件下保持可交换，例如使用层级回归模型。（第 126 页）

Gershoff 和 Grogan-Kaylor 在其元分析中处理的一个重要协变量，是研究类型。其元分析所依据的 76 篇论文同时包含组间设计和被试内设计。在 spank 数据中，我们用 between 和 within 向量对这一信息进行了哑变量编码。两者都是虚拟变量，并且 $\text{within} = 1 - \text{between}$。下面是它们的计数。

spank %>%
  count(between)

在模型中使用哑变量时，我倾向于让多数类别作为参照组。因此，我通常会用少数类别来命名这些变量。在这里，大多数情形都基于被试内设计。因此，我们接下来将 between 变量加入模型。顺便，我们也练习一下 0 + Intercept 这种语法。

b14.7 <- brm(
  data = spank,
  family = gaussian,
  d | se(se) ~ 0 + Intercept + between + (1 | study) + (1 | outcome),
  prior = c(prior(normal(0, 1), class = b),
            prior(cauchy(0, 1), class = sd)),
  iter = 4000, warmup = 1000, cores = 4, chains = 4,
  seed = 14,
  file = "fits/b14.07")

下面来看摘要输出。

print(b14.7)

让我们更仔细地看看 b_between 参数。

color <- viridis_pal(option = "B")(7)[5]

as_draws_df(b14.7) %>%
  ggplot(aes(x = b_between, y = 0)) +
  stat_halfeye(.width = c(.5, .95), color = "white", fill = color) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("Overall difference for between- vs within-participant designs")

这一差异并没有我原本预期的那么大。不过话说回来，我也并不是研究打屁股问题的专家。所以我又能知道多少呢？

对于这些数据，你还可以做其他事情。例如，你可以检验随年份变化的趋势，或者像作者在论文中所做的那样，区分不同严重程度的体罚。不过，我想我们已经讲得足够多，可以帮助你入门了。

如果你还想进一步了解这些方法，不妨看看 Vourre 的 Bayesian meta-analysis with R, Stan & brms。从他的博客中，你还可以学到一些额外技巧，例如如何绘制更传统风格的 forest plot，以及我们的贝叶斯 brms 方法与通过 metafor package 进行频率学元分析相比有何不同。你也可以参考前面提到的 Williams 及其同事的论文 Bayesian meta-analysis with weakly informative prior distributions，以获得在元分析 $\tau$ 参数中使用 half-Cauchy 先验的经验依据。

会话信息

sessionInfo()

评论