Robust Bayesian Meta-Analysis
In RoBMA: Robust Bayesian Meta-Analyses

source("vignette-cache.R", local = knitr::knit_global())

cached_fits <- vignette_cache(
  name     = "v21-robust-bayesian-meta-analysis",
  objects  = c("fit_RoBMA", "fit_RoBMA_custom"),
  packages = c("metadat")
)

knitr::opts_chunk$set(
  collapse    = TRUE,
  comment     = "#>",
  eval        = vignette_cache_eval(cached_fits),
  message     = FALSE,
  warning     = FALSE,
  fig.width   = 7,
  fig.height  = 3.5,
  dev         = "png",
  fig.retina  = 3
)
if (.Platform$OS.type == "windows") {
  knitr::opts_chunk$set(dev.args = list(type = "cairo"))
}

library("RoBMA")

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

vignette_cache_load(cached_fits)

library("RoBMA")

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

fit_RoBMA <- RoBMA(
  yi = yi, vi = vi, measure = "SMD",
  data = dat, seed = 1
)

fit_RoBMA_custom <- RoBMA(
  yi = yi, vi = vi, measure = "SMD",
  prior_bias = list(
    prior_PET(
      distribution = "Cauchy",
      parameters   = list(location = 0, scale = 1),
      truncation   = list(lower = 0, upper = Inf)
    ),
    prior_PEESE(
      distribution = "Cauchy",
      parameters   = list(location = 0, scale = 5),
      truncation   = list(lower = 0, upper = Inf)
    ),
    prior_weightfunction(
      "one-sided",
      steps   = c(0.025),
      weights = wf_cumulative(c(1, 1))
    )
  ),
  data = dat, seed = 1
)

vignette_cache_save(cached_fits)

In the Publication-Bias Adjustment vignette, we fit PET, PEESE, and a weight-function selection model on the ego-depletion dataset and obtained three different bias-corrected effect size estimates. The methods encoded different bias mechanisms: PET and PEESE corrected for the relationship between effect sizes and standard errors or sampling variances, and the selection model adjusted for publication bias operating on p-values. There is no consensus on which of the three is the right adjustment to apply. The Bayesian Model Averaging vignette highlighted how an analogous fixed-vs-random-effects question can be resolved by combining the candidate models in one ensemble and weighting them by their posterior model probabilities; the same idea applies here.

Robust Bayesian model-averaged meta-analysis (RoBMA) is the application of that idea to publication-bias adjustment [@maier2020robust; @bartos2021no]. The default RoBMA() ensemble averages over the presence vs absence of an effect, the presence vs absence of heterogeneity, and a no-bias spike against several weight-function and PET / PEESE bias models in one fit. Inclusion Bayes factors quantify evidence for each of the three components, and the model-averaged estimates propagate uncertainty across the bias mechanisms instead of committing to a single one.

In this vignette we fit the default RoBMA ensemble on the same ego-depletion data used in the Publication-Bias Adjustment vignette, walk through the summary, the marginal posterior weights of the bias components, and the interpret() shortcut. We list the available model_type presets and close with a custom ensemble that combines only the three single-model bias adjustments fit one at a time in that vignette.

Dataset

The example dataset is dat.lehmann2018 from the metadat package: 81 standardized mean differences from the ego-depletion meta-analysis with effect sizes yi and sampling variances vi already computed. For the matching single-model bias adjustments and metafor package parity calls, see the Publication-Bias Adjustment vignette.

library("RoBMA")

data("dat.lehmann2018", package = "metadat")
dat <- dat.lehmann2018
head(dat[, c("Short_Title", "Year", "yi", "vi")])

Default RoBMA Ensemble

RoBMA() with default arguments fits the RoBMA-PSMA ensemble [@bartos2021no]: spike-and-slab prior distributions on the effect and heterogeneity, and a 9-component bias mixture (no-bias spike, six weight functions, PET, PEESE). The default prior distributions on mu, tau, and the bias parameters match the single-model fits in the Publication-Bias Adjustment vignette.

fit_RoBMA <- RoBMA(
  yi = yi, vi = vi, measure = "SMD",
  data = dat, seed = 1
)

summary(fit_RoBMA)

The Components table reports inclusion Bayes factors for the three components. Model-averaged estimates combine all 36 sub-models in the ensemble (effect $\times$ heterogeneity $\times$ bias) weighted by their posterior probabilities; we report the bias parameters only for the bias-adjusted sub-models.

summary_models() breaks the mixture prior distributions of the three components down into their individual sub-models and reports their prior weight, posterior weight, and individual inclusion Bayes factor.

summary_models(fit_RoBMA)

The posterior model weights show which sub-models the data prefer within each component's mixture. The bias side of the ensemble concentrates posterior mass on PET (posterior probability $0.77$, individual inclusion BF $\approx 23$) with a smaller share on PEESE; the weight functions and the no-bias spike all receive negligible posterior weight. The inclusion BF column is the Bayes factor for that one sub-model against the rest of its mixture, so it does not need to agree with the component-level inclusion BF reported by summary().

plot() overlays the prior and the model-averaged posterior of any model parameter. On mu the spike at zero (right-axis arrow) carries the posterior probability that the effect is null and the slab carries the rest; the grey shapes are the corresponding prior distribution.

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

interpret() distills the components and the pooled estimates into prose. With scope = "all", the bias section is included.

interpret(fit_RoBMA, scope = "all")

Preset Ensembles

model_type selects a preset bias-model ensemble. Each preset combines the no-bias spike with a different alternative set; the prior weight on "any bias" is $0.5$ in every preset.

| model_type | Alternative bias models | |:-------------|:-------------------------------------------------------------------------------------| | "PSMA" | 6 weight functions, PET, PEESE (the default RoBMA-PSMA ensemble [@bartos2021no]) | | "PP" | PET and PEESE only | | "6w" | 6 weight functions only (no PET / PEESE) | | "2w" | 2 two-sided weight functions only |

The PEESE scale is rescaled by $\text{UISD}{\text{ratio}} = \text{UISD}{\text{SMD}} / \text{UISD}_{\text{measure}}$ so the prior distribution describes a comparable bias on each effect-size scale (the ratio is $1$ for measure = "SMD"). Bias-adjustment prior distributions are not affected by rescale_priors. Use print_prior(fit, parameter = "bias") on a RoBMA(..., only_priors = TRUE) fit to inspect the alternative components and prior weights inside any preset:

fit_priors <- RoBMA(
  yi = yi, vi = vi, measure = "SMD",
  model_type = "PP",
  data = dat, only_priors = TRUE
)
print_prior(fit_priors, parameter = "bias")

Custom Bias Ensembles

Beyond the presets, we can specify the bias model space directly via prior_bias. This is useful when the publication process under study has a specific shape (for example, two-sided selection only or a single weight-function step at a non-default cutoff), or for sensitivity analyses against the preset choice.

The example below builds a 3-model bias ensemble out of the three single-model bias adjustments from the Publication-Bias Adjustment vignette (default bPET(), default bPEESE(), and default bselmodel()), so the ensemble averages over the same three bias mechanisms that vignette fit one at a time, plus the no-bias spike.

fit_RoBMA_custom <- RoBMA(
  yi = yi, vi = vi, measure = "SMD",
  prior_bias = list(
    prior_PET(
      distribution = "Cauchy",
      parameters   = list(location = 0, scale = 1),
      truncation   = list(lower = 0, upper = Inf)
    ),
    prior_PEESE(
      distribution = "Cauchy",
      parameters   = list(location = 0, scale = 5),
      truncation   = list(lower = 0, upper = Inf)
    ),
    prior_weightfunction(
      "one-sided",
      steps   = c(0.025),
      weights = wf_cumulative(c(1, 1))
    )
  ),
  data = dat, seed = 1
)

interpret(fit_RoBMA_custom, scope = c("components", "estimates"))

The two ensembles agree on all three component inclusion conclusions: strong evidence for heterogeneity, moderate evidence against the effect, and strong evidence in favour of publication bias. The model-averaged effect on mu is essentially zero in both fits with overlapping credible intervals. The agreement reflects the bias-mixture composition reported by summary_models() above: the default ensemble already concentrates posterior weight on PET, so dropping the weight functions that received negligible weight does not move the model-averaged inference.

The prior probability of any bias is $0.75$ in the custom fit and $0.5$ in the default, because the three user-supplied alternative prior distributions all default to prior_weights = 1 against the no-bias spike's weight of $1$. The component inclusion Bayes factor is invariant to this asymmetry; pass prior_weights explicitly on each prior in prior_bias to match a different prior odds split.

Other Inference Helpers

RoBMA() returns an object of class c("RoBMA", "brma"), so all summary(), plot(), predict(), funnel(), marginal_means(), rstudent(), qqnorm(), and influence() methods documented in the Bayesian Meta-Analysis vignette work the same way. For single-model bias adjustment without averaging see the Publication-Bias Adjustment vignette; for the model-averaged workflow without publication-bias adjustment see Bayesian Model Averaging; the prior distributions on mu, tau, and the moderators are documented in Prior Distributions.