Robust Bayesian Meta-Analysis (RoBMA)

Introduction

Any single publication-bias correction method embeds strong assumptions about the nature of the selection mechanism. PET-PEESE assumes selection is correlated with standard error. Weight-function models assume selection operates through p-values. Choosing between them requires subjective judgment, and reporting only the best-fitting model ignores model uncertainty.

Robust Bayesian Meta-Analysis (RoBMA; Maier, Bartoš & Wagenmakers, 2023) addresses this by averaging over a family of meta-analytic models rather than committing to one. RoBMA produces a single model-averaged posterior for the true effect that reflects uncertainty about both the heterogeneity structure and the publication-bias mechanism.

bayesma implements RoBMA via robma().

Model family

RoBMA fits a grid of models formed by crossing:

Heterogeneity: common-effect vs. random-effects
Publication bias: none vs. PET-PEESE vs. selection weight function

This yields a default grid of $2 \times 3 = 6$ component models. Each model $M_k$ is fitted separately.

Model	Heterogeneity	Bias correction
$M_1$	Common effect	None
$M_2$	Random effects	None
$M_3$	Common effect	PET-PEESE
$M_4$	Random effects	PET-PEESE
$M_5$	Common effect	Selection weight
$M_6$	Random effects	Selection weight

Custom model grids can be specified via the bias_models and heterogeneity arguments.

Model averaging

Each component model is assigned a prior probability $p(M_k)$ (equal by default: $1/6$ ). After fitting, Bayes factors are used to compute posterior model probabilities:

$p(M_k \mid \mathbf{y}) \propto p(\mathbf{y} \mid M_k) \cdot p(M_k)$

where $p(\mathbf{y} \mid M_k)$ is the marginal likelihood of model $k$ , estimated via bridge sampling.

The model-averaged posterior for any quantity $\psi$ (e.g., $\mu$ ) is:

$p(\psi \mid \mathbf{y}) = \sum_{k=1}^{K} p(M_k \mid \mathbf{y}) \cdot p(\psi \mid \mathbf{y}, M_k)$

Bayes factors for heterogeneity and bias

RoBMA computes two summary Bayes factors by collapsing over the model grid:

Bayes factor for heterogeneity (random-effects vs. common-effect):

$\text{BF}_\tau = \frac{P(\tau > 0 \mid \mathbf{y})}{P(\tau = 0 \mid \mathbf{y})}$

Bayes factor for publication bias:

$\text{BF}_b = \frac{P(\text{bias} \mid \mathbf{y})}{P(\text{no bias} \mid \mathbf{y})}$

These are the primary inferential summaries: $\text{BF}_\tau$ quantifies the evidence that heterogeneity exists; $\text{BF}_b$ quantifies the evidence that publication bias is present.

Priors

Component model priors:

$\mu \sim \mathcal{N}(0,\; 1), \qquad \tau \sim \text{Half-Cauchy}(0,\; 0.5)$

PET-PEESE coefficient:

$\beta_\text{PET} \sim \mathcal{N}(0,\; 1)$

Weight-function selection weights:

$\omega_j \sim \text{Dirichlet}(\mathbf{1})$

Priors can be modified via robma_default_priors() and passed to robma().

Fitting RoBMA

#| eval: false
fit_robma <- robma(data)

robma_output(fit_robma)
robma_table(fit_robma)

robma_output() prints the model-averaged summary. robma_table() returns a gt table with posterior model probabilities, the Bayes factors for heterogeneity and bias, and the model-averaged posterior for $\mu$ .

Interpreting RoBMA output

Model-averaged posterior for $\mu$

The model-averaged posterior for $\mu$ integrates over all component models, weighted by their posterior model probabilities. This is the primary bias-corrected estimate.

Posterior model probabilities

robma_table() reports the posterior model probability for each component model $M_k$ . Large posterior probability on models without a bias component ( $M_1$ , $M_2$ ) indicates little evidence for publication bias.

Bayes factors

$\text{BF}$	Interpretation
$< 1$	Evidence against the presence of heterogeneity / bias
$1$ – $3$	Anecdotal
$3$ – $10$	Moderate
$10$ – $30$	Strong
$> 30$	Very strong

Sensitivity analysis

Bridge sampling estimates have Monte Carlo error. robma_sensitivity() reruns the bridge sampler with different random seeds and reports the variance of the log-marginal-likelihood estimates:

#| eval: false
fit_robma <- robma(data, iter_sampling = 5000)  # more samples → more stable BF
robma_sensitivity(fit_robma)

Posterior model probabilities that change substantially across seeds indicate that the MCMC chains need more samples before bridge sampling is reliable.

Custom model grids

#| eval: false
fit_custom <- robma(
  data,
  bias_models     = c("none", "selection_weight"),  # exclude PET-PEESE
  heterogeneity   = c("random"),                    # random-effects only
  priors_effect   = normal(0, 0.5),
  null_range      = c(-0.1, 0.1)
)

The null_range argument specifies a region of practical equivalence around zero. Model probabilities are then computed for models that predict the effect is in vs. outside this range.

Limitations

Bridge sampling requires well-mixed MCMC chains. Use iter_sampling ≥ 4000 and check convergence diagnostics before interpreting Bayes factors.
RoBMA is computationally intensive — fitting six component models with bridge sampling takes several minutes per dataset. Use robma(data, parallel = TRUE) to fit component models in parallel.
The default model grid assumes that the bias direction is positive (published effects are inflated upward). For outcomes where bias deflates effects, specify appropriate priors on the weight functions.

References

Maier M, Bartoš F, Wagenmakers EJ (2023). Robust Bayesian meta-analysis: Addressing publication bias with model-averaging. Psychological Methods, 28(1), 107–122.

Bartoš F, Maier M, Wagenmakers EJ, Doucouliagos H, Stanley TD (2023). Robust Bayesian meta-analysis: Model-averaging across formulations of publication bias. Psychological Methods.