This vignette accompanies the manuscript Robust Bayesian Multilevel Meta-Analysis: Adjusting for Publication Bias in the Presence of Dependent Effect Sizes preprinted at PsyArXiv (Bartoš et al., 2026).
This vignette reproduces the first example from the manuscript. For multilevel meta-regression with moderators, see the Multilevel Robust Bayesian Model-Averaged Meta-Regression vignette.
Multilevel Robust Bayesian Meta-Analysis (RoBMA) extends the standard RoBMA framework to datasets containing multiple effect sizes from the same studies. We demonstrate the method using data from Johnides et al. (2025), who meta-analyzed 412 effect sizes from 128 studies investigating secondary benefits of family-based treatments for childhood disorders.
Multilevel RoBMA is appropriate for datasets that contain multiple
effect sizes from the same studies. The multilevel approach explicitly
models the clustering structure through the cluster
argument, which:
We load the package and examine the dataset structure.
The dataset contains effect sizes (d), standard errors
(se), and study identifiers (study).
Multiple rows share the same study name, indicating that these effect sizes come from the same study. For example, the first five effect sizes are from “Price et al. (2012)” and might be more similar to each other than to effect sizes from other studies.
To reproduce the RoBMA 3.6.1 version results, which used a different prior distribution formulation and scaling settings, we need to specify the prior distributions manually (see the Prior Distributions vignette for more details on specifying prior distributions).
We specify the effect sizes as yi and standard errors as
sei argument. Setting measure = "SMD"
specifies that we are using standardized mean differences as the
effect-size input, which is important for the prior distribution
specification (here the effect size and heterogeneity prior
distributions are specified explicitly since they differ from the
default specification in the 3.6.1 version; the publication-bias
specification remains the same). The cluster argument
specifies which effect sizes belong together.
fit <- RoBMA(
yi = d, sei = se, measure = "SMD", cluster = study,
prior_effect = prior_effect_361,
prior_heterogeneity = prior_heterogeneity_361,
adapt = 5000, burnin = 5000, sample = 10000,
parallel = TRUE, seed = 1, data = Johnides2025
)Note: This model takes 10-15 minutes with parallel processing enabled. Without parallel processing, expect over an hour.
The summary() output contains two main sections.
This table shows Bayes factors testing the presence of effect, heterogeneity, and publication bias. Each row displays the prior probability, posterior probability after observing the data, and the inclusion Bayes factor comparing models with versus without that component.
In our example, the inclusion BF for the effect is 0.927, indicating
weak evidence against an effect. The inclusion BF for heterogeneity and
publication bias are reported as Inf, i.e., beyond the
numerical precision of the estimation algorithm (> 10⁶), indicating
extreme evidence for both components.
We also notice a warning about the effective sample size (ESS) being below the set target. The difference is, however, not substantial, and we can safely ignore it.
This section provides meta-analytic estimates averaged across all
models, weighted by their posterior probabilities. The average effect
size is mu, the overall heterogeneity is tau,
and the heterogeneity allocation parameter is rho. The
heterogeneity allocation parameter rho indicates the
proportion of heterogeneity allocated to the cluster-level component:
rho = 0 means all heterogeneity is within clusters,
rho = 1 means all heterogeneity is between clusters, and
rho = 0.5 means equal split. The remaining parameters
summarize the weight function and PET/PEESE regression for publication
bias.
In our example, we find a small effect size (d = 0.050, 95%
CI [0.000, 0.173]), substantial heterogeneity (τ = 0.40, 95% CI [0.348,
0.458]), and nearly balanced heterogeneity allocation (ρ = 0.461, 95% CI
[0.306, 0.603]). The large heterogeneity combined with the small pooled
effect suggests substantial variation in true effects. This distribution
of true effect sizes can be obtained from the
summary_heterogeneity() function:
The wide prediction interval (d = -0.752 to 0.845) quantifies the degree of this heterogeneity: some studies may show benefits while others show harm.
To understand which publication bias mechanisms the data support, we examine the posterior distribution across model types:
Essentially all of the posterior probability is allocated to selection models (weight functions). The publication-bias adjustment therefore reflects selective reporting rather than small-study effects (PET/PEESE regression); specifically, one-sided selection operating on marginally significant p-values and on the direction of the effect.
The weight function shows how publication probability varies across p-value ranges:
The plot displays one-sided p-value cutoffs (x-axis) against relative publication probability (y-axis), averaged across weight function specifications. A flat line at 1.0 indicates no publication bias; values below 1.0 indicate suppression.
In our example, we notice a small drop in the relative publication probability at the cutoff corresponding to p = 0.05 (one-sided, i.e., 0.10 two-sided: selection on marginally significant p-values) and a much sharper drop corresponding to the direction of the effect (p = 0.50).
To understand the importance of accounting for nested structure, we compare our results with a standard single-level RoBMA that ignores dependencies among effect sizes from the same study.
fit_simple <- RoBMA(
yi = d, sei = se, measure = "SMD",
prior_effect = prior_effect_361,
prior_heterogeneity = prior_heterogeneity_361,
adapt = 5000, burnin = 5000, sample = 10000,
parallel = TRUE, seed = 1, data = Johnides2025
)The single-level model differs by:
cluster: it treats all effect sizes as
independent,For the Johnides et al. data, the single-level RoBMA finds strong evidence for the absence of an effect, while the multilevel model provides only weak evidence against it. Properly accounting for data structure leads to more conservative and appropriate inferences.