Hierarchical Bayesian Model-Averaged Meta-Analysis

Example Data Set

We use the dat.konstantopoulos2011 data set from the metadat R package (Thomas et al., 2019) that is used for the same functionality in the metafor (Wolfgang, 2010) R package. We roughly follow the example in the data set’s help file, ?dat.konstantopoulos2011. The data set consists of 56 studies estimating the effects of modified school calendars on students achievement. The 56 studies were run in individual schools, which can be grouped into 11 districts. We might expect more similar effect size estimates from schools in the same district – in other words, the effect size estimates from the same district might not be completely independent. Consequently, we might want to adjust for this dependency (clustering) between the effect size estimates to draw a more appropriate inference.

First, we load the data set, assign it to the dat object, and inspect the first few rows.

data("dat.konstantopoulos2011", package = "metadat")
dat <- dat.konstantopoulos2011

head(dat)

In the following analyses, we use the following variables:

• yi, standardized mean differences,
• vi, sampling variances of the standardized mean differences,
• district, district id which distinguishes among the districts,
• and school, that distinguishes among different schools within the same district.

Frequentist Hierarchical Meta-Analysis with metafor

We follow the data set’s help file and fit a simple random effects meta-analysis using the rma() function from metafor package. This model ignores the dependency between effect size estimates. We use this simple model as our starting point and as a comparison with the later models.

fit_metafor.0 <- metafor::rma(yi = yi, vi = vi, data = dat)
fit_metafor.0

The model summary returns a small but statistically significant effect size estimate μ = 0.128 (se = 0.044) and a considerable heterogeneity estimate τ = 0.297.

We extend the model to account for the hierarchical structure of the data, i.e., schools within districts, by using the rma.mv() function from the metafor package and extending it with the random = ~ school | district argument.

fit_metafor <- metafor::rma.mv(yi, vi, random = ~ school | district, data = dat)
fit_metafor

We find that accounting for the hierarchical structure of the data results in (1) a slightly larger effect size estimate (μ = 0.187) and (2) larger standard error of the effect size estimate (se = 0.085). The larger standard error is a natural consequence of accounting for the dependency between the effect sizes. Since the effect sizes are dependent, they do not contribute independent information. Specifying the hierarchical model then accounts for the dependency by estimating similarity between the estimates from the same cluster (school) and discounting the information borrowed from each estimate. The estimate of the similarity among estimates from the same cluster is summarized in the \rho = 0.666 estimate.

Specifications of Hierarchical Meta-Analysis

We specify a simple hierarchical meta-analytic model (see Konstantopoulos (2011) for an example). Using distributional notation, we can describe the data generating process as a multi-stage sampling procedure. In a nutshell, we assume the existence of an overall mean effect μ. Next, we assume that the effect sizes in each district k = 1, …, K, γ_k, systematically differ from the mean effect, with the variance of the district-level effects summarized with heterogeneity τ_b (as between). Furthermore, we assume that the true effects θ_k, j of each study j = 1, …J_k systematically differ from the district-level effect, with the variance of the study effects from the district-level effect summarized with heterogeneity τ_w (as within). Finally, the observed effect sizes y_k, j that differ from the true effects y_k, j due to random errors se_k, j.

Mathematically, we can describe such a model as: $$ \begin{aligned} \gamma_k &\sim \text{N}(\mu, \tau_b^2),\\ \theta_{k,j} &\sim \text{N}(\gamma_k, \tau_w^2),\\ y_{k,j} &\sim \text{N}(\theta_{k,j}, \text{se}_{k,j}).\\ \end{aligned} $$ Where N() denotes a normal distribution with mean and variance.

Conveniently, and with a bit of algebra, we do not need to estimate the district-level and true study effects. Instead, we marginalize them out, and we sample the observed effect sizes from each district y_k, . directly from a multivariate normal distributions, MN(), with a common mean μ and covariance matrix S: $$ \begin{aligned} y_{k,.} &\sim \text{MN}(\mu, \text{S}),\\ \text{S} &= \begin{bmatrix} \tau_b^2 + \tau_w^2 + \text{se}_1^2 & \tau_w^2 & \dots & \tau_w^2 \\ \tau_w^2 & \tau_b^2 + \tau_w^2 + \text{se}_2^2 & \dots & \tau_w^2 \\ \dots & \dots & \dots & \dots \\ \tau_w^2 & \tau_w^2 & \dots & \tau_b^2 + \tau_w^2 + \text{se}_{J_k}^2 & \\ \end{bmatrix}. \end{aligned} $$ The random effects marginalization is helpful as it allows us to sample much fewer parameters from the posterior distribution (which significantly simplifies marginal likelihood estimation via bridge sampling). Furthermore, the marginalization allows us to properly specify selection model publication bias adjustment models – the marginalization propagates the selection process up through all the sampling steps at once (we cannot proceed with the sequential sampling as the selection procedure on the observed effect sizes modifies the sampling distributions of all the preceding levels).

We can further re-parameterize the model by performing the following substitution, $$ \begin{aligned} \tau^2 &= \tau_b^2 + \tau_w^2,\\ \rho &= \frac{\tau_w^2}{\tau_b^2 + \tau_w^2}, \end{aligned} $$ and specifying the covariance matrix using the inter-study correlation ρ, total heterogeneity τ, and the standard errors se_.: $$ \begin{aligned} \text{S} &= \begin{bmatrix} \tau^2 + \text{se}_1^2 & \rho\tau^2 & \dots & \rho\tau^2 \\ \rho\tau^2 & \tau^2 + \text{se}_2^2 & \dots & \rho\tau^2 \\ \dots & \dots & \dots & \dots \\ \rho\tau^2 & \rho\tau^2 & \dots & \tau^2 + \text{se}_{J_k}^2 & \\ \end{bmatrix}. \end{aligned} $$ This specification corresponds to the compound symmetry covariance matrix of random effects, the default settings in the metafor::rma.mv() function. More importantly, it allows us to easily specify prior distributions on the correlation coefficient ρ and the total heterogeneity τ.

Hierarchical Bayesian Model-Averaged Meta-Analysis with RoBMA

Before we estimate the complete Hierarchical Bayesian Model-Averaged Meta-Analysis (hBMA) with the RoBMA package, we quickly repeat the simpler models we estimated with the metafor package in the previous section.

fit.0 <- RoBMA(d = dat$yi, v = dat$vi,
               priors_effect_null        = NULL,
               priors_heterogeneity_null = NULL,
               priors_bias               = NULL,
               parallel = TRUE, seed = 1)

summary(fit.0, type = "individual")

fit <- RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district,
             priors_effect_null        = NULL,
             priors_heterogeneity_null = NULL,
             priors_bias               = NULL,
             parallel = TRUE, seed = 1)

summary(fit, type = "individual")

par(mar = c(2, 4, 0, 0))
plot(fit, parameter = "rho", prior = TRUE)

fit_BMA <- RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district,
                 priors_bias = NULL,
                 parallel = TRUE, seed = 1)

summary(fit_BMA)

hierarchical_test <- RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district,
                           priors_heterogeneity_null = NULL,
                           priors_hierarchical_null = prior(distribution = "spike", parameters = list("location" = 0)),
                           priors_bias = NULL,
                           parallel = TRUE, seed = 1)

summary(hierarchical_test)

Summary

We illustrated how to estimated estimate a hierarchical Bayesian model-averaged meta-analysis using the RoBMA package. The hBMA model allows us to test for the presence vs absence of the effect and heterogeneity while simultaneously adjusting for clustered effect size estimates. While the the current implementation allows us to draw a fully Bayesian inference, incorporate prior information, and acknowledge model uncertainty, it has a few limitations in contrast to the metafor package. E.g., the RoBMA package only allows a simple nested random effects (i.e., estimates within studies, schools within districts etc). The simple nesting allows us to break the full covariance matrix into per cluster block matrices which speeds the already demanding computation. Furthermore, the computation complexity significantly increases when considering selection models as we need to compute an exponentially increasing number of multivariate normal probabilities with the increasing cluster size (existence of clusters with more than 4 studies makes the current implementation impractically). The current limitations are however not the end of the road as we explore other approaches (e.g., only specifying PET-PEESE style publication bias adjustment, and other dependency adjustments) in a future vignette.

Footnotes

¹ We could also model-average across the hierarchical structure assuming fixed effect models, i.e., τ ∼ f(.) and ρ = 1. However specifying such a model ensemble is a beyond scope of this vignette, see Custom ensembles vignette for some hints.