Bayesian Model Averaging

In meta-analysis, analysts often have to balance different assumptions about the data-generating process and decide what kind of model to fit, which covariates to include, and what publication-bias adjustment method to use. Bayesian model averaging (BMA) tries to alleviate this problem by combining competing meta-analytic models based on their prior predictive performance. Models that predicted the data better receive higher posterior model probability and more weight in the overall inference. The model-averaged estimates allow us to propagate uncertainty about the selected model and inclusion Bayes factors quantify evidence across sets of models (Hoeting et al., 1999).

In this vignette we illustrate Bayesian model-averaged meta-analysis using BMA() on the dat.baskerville2012 dataset to highlight the basics of Bayesian model-averaging using the RoBMA R package (Bartoš & Maier, 2020). In a later vignette Robust Bayesian Meta-Analysis, we illustrate robust Bayesian meta-analysis that extends this concept to model-averaging across publication-bias adjustment models. All models use the default prior distributions; see the Prior Distributions vignette for the prior definitions and customization options.

We start by estimating simple fixed-effect and random-effects Bayesian meta-analytic models, turn the model-comparison into a model-averaged fit, and finally also average over the presence vs absence of the effect. Finally, we extend the model into a Bayesian model-averaged meta-regression. Note that we do not demonstrate all of the diagnostics and summary features available to all models fitted via the RoBMA R package (see the Bayesian Meta-Analysis vignette for more details).

Dataset

The example dataset is distributed with the metadat package as dat.baskerville2012 and contains 23 randomized and controlled trials of practice-facilitation interventions in primary care. The effect size smd is the standardized mean difference for adherence to evidence-based primary-care guidelines and se is its standard error. The binary tailor indicator records whether the intervention was tailored to the participating practice (1) or applied as a uniform protocol (0); we treat it as a factor and use it as a moderator later in the vignette.

library("RoBMA")

data("dat.baskerville2012", package = "metadat")
dat <- dat.baskerville2012
dat$tailor <- as.factor(dat$tailor)
head(dat[, c("year", "smd", "se", "tailor")])

Comparing Fixed and Random Effects

metafor::rma.uni() fits the conventional random-effects model on the standardized mean differences.

res <- metafor::rma.uni(
  yi = smd, sei = se,
  data = dat, method = "REML"
)
res

The Bayesian counterpart uses brma() with the default prior distributions: a Normal(0, 0.71) unit-information prior distribution on the effect and a Normal(0, 0.35)[0, Inf] half-normal prior distribution on the heterogeneity.

fit_RE <- brma(
  yi = smd, sei = se, measure = "SMD",
  data = dat, seed = 1
)
fit_RE <- add_marglik(fit_RE)

summary(fit_RE, include_mcmc_diagnostics = FALSE)

Setting the heterogeneity prior distribution to a spike at zero turns the random-effects model into a fixed-effect model.

fit_FE <- brma(
  yi = smd, sei = se, measure = "SMD",
  prior_heterogeneity = prior("spike", list(location = 0)),
  data = dat, seed = 1
)
fit_FE <- add_marglik(fit_FE)

summary(fit_FE, include_mcmc_diagnostics = FALSE)

With marginal likelihoods attached to the models via add_marglik(), bf() returns the Bayes factor between the two fits.

bf(fit_RE, fit_FE)

The Bayes factor is close to one and neither model is decisively preferred. Rather than picking one of the two models, we can model-average over them and weight the inference by the posterior model probabilities.

From Two Fits to Model Averaging

BMA() fits a mixture of the two models in one call. Each top-level prior-distribution argument has a paired _null version; when we supply both, the model averages over that component. Setting prior_effect_null = NULL drops the null component on the effect (spike at zero) and leaves a 2-model average over presence vs absence of heterogeneity, which is the BMA analogue of the comparison above.

fit_BMA_he <- BMA(
  yi = smd, sei = se, measure = "SMD",
  prior_effect_null = NULL,
  data = dat, seed = 1
)

summary(fit_BMA_he)

The Components table reports the posterior probability and inclusion Bayes factor for heterogeneity. This BF closely matches bf(fit_RE, fit_FE) above up to MCMC noise. The two routes compute the same quantity differently: bf() bridge-samples marginal likelihoods from two separate fits, while BMA() runs a single MCMC over the joint model space with a model indicator as a latent variable and recovers the posterior model probabilities from the indicator frequencies (Lodewyckx et al., 2011). Adding a third or fourth component to this product-space sampler costs one more indicator instead of a set of full model fits and bridge-sampling runs, which is what makes the larger ensembles in RoBMA() tractable.

Adding Effect Averaging

The default BMA() call also averages over the presence vs absence of the effect, giving a 4-model ensemble (effect \(\times\) heterogeneity, each present or absent).

fit_BMA <- BMA(
  yi = smd, sei = se, measure = "SMD",
  data = dat, seed = 1
)

summary(fit_BMA)

Both inclusion Bayes factors come from the same MCMC. The effect BF is overwhelming on this dataset; the heterogeneity BF is the same as before. The model-averaged estimates are spike-and-slab mixtures: a point mass at the null value with weight equal to the posterior probability that the parameter is zero, and a slab carrying the remaining weight.

plot() draws both components together. The arrow at zero is the spike’s posterior probability (right axis), the curve is the slab’s density, and the grey shapes show the corresponding visualization for the prior distribution.

par(mar = c(4, 4, 1, 4))
plot(fit_BMA, parameter = "mu", prior = TRUE, xlim = c(-1, 1.5), lwd = 2)

par(mar = c(4, 4, 1, 4))
plot(fit_BMA, parameter = "tau", prior = TRUE, xlim = c(0, 1), lwd = 2)

The Model-averaged estimates summarize the marginal posterior (the full mixture). For estimates conditional on the parameter being non-zero (slab only) use the conditional = TRUE argument. The marginal and conditional summaries diverge for \(\tau\), which keeps substantial posterior weight on the spike, and coincide for \(\mu\), where the spike’s posterior weight is essentially zero.

summary(fit_BMA, conditional = TRUE)

Meta-Regression with Model Averaging

Adding a moderator extends the averaging to its inclusion. With mods = ~ tailor, BMA() averages over 8 sub-models (effect \(\times\) heterogeneity \(\times\) tailor coefficient, each present or absent).

res_mod <- metafor::rma.uni(
  yi = smd, sei = se, mods = ~ tailor,
  data = dat, method = "REML"
)
res_mod

fit_BMA_mod <- BMA(
  yi = smd, sei = se, measure = "SMD",
  mods = ~ tailor,
  data = dat, seed = 1
)

summary(fit_BMA_mod)

The Components table now also lists tailor and its inclusion BF. Because we average the moderation across models in which the slope is zero, regplot() shrinks the difference towards zero when the inclusion BF is modest.

par(mar = c(4, 4, 1, 1))
regplot(fit_BMA_mod, mod = "tailor")

Alternatively, we can directly compute the estimated marginal means via the marginal_means() function and visualize the prior and posterior distributions at each moderator level.

par(mar = c(4, 4, 1, 1))
emm <- marginal_means(fit_BMA_mod)
emm

plot(emm, parameter = "tailor", prior = TRUE)

Customizing Prior Components

In Bayesian model-averaged meta-analysis, every component has a paired _null argument (i.e. prior_effect and prior_effect_null). Either side can be NULL (remove the component), a single prior distribution, or a list of prior distributions (a sub-mixture). Common customizations are tightening or widening the slab, replacing the point-zero null with a region of practical equivalence, or forcing a moderator into every model. The examples below use only_priors = TRUE to print the prior structure without fitting the models.

Tighten the effect slab with rescale_priors = 0.5:

fit_priors <- BMA(
  yi = smd, sei = se, measure = "SMD",
  rescale_priors = 0.5,
  data = dat, only_priors = TRUE
)
print_prior(fit_priors, parameter = "mu")

Replace the point-zero null on the effect with a small spread around zero (a perinull) to test against a region of practical equivalence:

fit_priors <- BMA(
  yi = smd, sei = se, measure = "SMD",
  prior_effect_null = prior("normal", list(mean = 0, sd = 0.05)),
  data = dat, only_priors = TRUE
)
print_prior(fit_priors, parameter = "mu")

Force tailor into every model by dropping its null:

fit_priors <- BMA(
  yi = smd, sei = se, measure = "SMD",
  mods = ~ tailor,
  prior_mods_null = list(tailor = NULL),
  data = dat, only_priors = TRUE
)
print_prior(fit_priors, parameter_mods = "tailor")

A tighter half-normal slab on \(\tau\):

fit_priors <- BMA(
  yi = smd, sei = se, measure = "SMD",
  prior_heterogeneity = prior(
    distribution = "normal",
    parameters   = list(mean = 0, sd = 0.25),
    truncation   = list(lower = 0)
  ),
  data = dat, only_priors = TRUE
)
print_prior(fit_priors, parameter = "tau")

Empirical informed prior distributions based on the Cochrane database are available via prior_informed_field and prior_informed_subfield; see Informed Bayesian Model-Averaged Meta-Analysis in Medicine.

Other `brma` Features

A BMA() fit is also a brma object: summary(), plot(), print_prior(), predict(), marginal_means(), funnel(), rstudent(), qqnorm(), and influence() work the same way as on a single-model fit; see Bayesian Meta-Analysis. Prior-distribution families and informed prior distributions are covered in Prior Distributions. For binary or count outcomes, BMA.glmm() does the same product-space averaging on a binomial-normal or Poisson-normal likelihood. For model averaging that also adjusts for publication bias, see Robust Bayesian Meta-Analysis.

References

Bartoš, F., & Maier, M. (2020). RoBMA: An R package for robust Bayesian meta-analyses. https://CRAN.R-project.org/package=RoBMA

Hoeting, J. A., Madigan, D., Raftery, A. E., & Volinsky, C. T. (1999). Bayesian model averaging: A tutorial. Statistical Science, 14(4), 382–401. https://doi.org/10.1214/ss/1009212519

Lodewyckx, T., Kim, W., Lee, M. D., Tuerlinckx, F., Kuppens, P., & Wagenmakers, E.-J. (2011). A tutorial on Bayes factor estimation with the product space method. Journal of Mathematical Psychology, 55(5), 331–347. https://doi.org/10.1016/j.jmp.2011.06.001