---
title: "Robust Bayesian Model-Averaged Meta-Regression"
author: "František Bartoš"
date: "11th of December 2024 (updated: 29th of April 2026)"
output:
  rmarkdown::html_vignette:
    self_contained: yes
bibliography: ../inst/REFERENCES.bib
csl: ../inst/apa.csl
link-citations: true
vignette: >
  %\VignetteIndexEntry{Robust Bayesian Model-Averaged Meta-Regression}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown_notangle}
---

```{r child = "_vignette-nowrap.md", echo = FALSE, eval = TRUE}
```

```{r setup, include = FALSE}
source("vignette-cache.R", local = knitr::knit_global())

cached_fits <- vignette_cache(
  name     = "v31-robma-metaregression",
  objects  = c("fit_BMA", "fit_RoBMA"),
  packages = c("metafor", "emmeans")
)

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"))
}
```

```{r load-models, include = FALSE}
library("RoBMA")

data("Andrews2021", package = "RoBMA")

Andrews2021_z <- metafor::escalc(
  ri = ri, ni = ni,
  measure = "ZCOR", data = Andrews2021
)

vignette_cache_load(cached_fits)
```

```{r precompute-models, include = FALSE, eval = FALSE}
library("RoBMA")

data("Andrews2021", package = "RoBMA")

Andrews2021_z <- metafor::escalc(
  ri = ri, ni = ni,
  measure = "ZCOR", data = Andrews2021
)

prior_effect_361 <- prior("normal", list(mean = 0, sd = 0.5))
prior_heterogeneity_361 <- prior("invgamma", list(shape = 1, scale = 0.075))

fit_BMA <- BMA(
  yi = yi, vi = vi, measure = "ZCOR",
  mods = ~ measure + age,
  prior_effect        = prior_effect_361,
  prior_heterogeneity = prior_heterogeneity_361,
  data = Andrews2021_z, parallel = TRUE, seed = 1
)

fit_RoBMA <- RoBMA(
  yi = yi, vi = vi, measure = "ZCOR",
  mods = ~ measure + age,
  prior_effect        = prior_effect_361,
  prior_heterogeneity = prior_heterogeneity_361,
  data = Andrews2021_z, parallel = TRUE, seed = 1
)

vignette_cache_save(cached_fits)
```

**This vignette accompanies the manuscript [Robust Bayesian Meta-Regression: Model-Averaged Moderation Analysis in the Presence of Publication Bias](https://doi.org/10.1037/met0000737) published at *Psychological Methods* [@bartos2023robust].**

Robust Bayesian model-averaged meta-regression (RoBMA-reg) extends the robust Bayesian model-averaged meta-analysis (RoBMA) by including covariates in the meta-analytic model. 
RoBMA-reg allows for estimating and testing the moderating effects of study-level covariates on the meta-analytic effect in a unified framework (e.g., accounting for uncertainty in the presence vs. absence of the effect, heterogeneity, and publication bias). 
This vignette illustrates how to fit a robust Bayesian model-averaged meta-regression using the `RoBMA` R package. 
We reproduce the example from @bartos2023robust, who re-analyzed a meta-analysis of the effect of household chaos on child executive functions with the mean age and assessment type covariates based on @andrews2021examining's meta-analysis. 

First, we fit a frequentist meta-regression using the `metafor` package.
Second, we explain the Bayesian meta-regression model specification, the default prior distributions for continuous and categorical moderators, and the standardized effect-size input specification.
Third, we estimate Bayesian model-averaged meta-regression (without publication bias adjustment).
Finally, we estimate the complete robust Bayesian model-averaged meta-regression.

## Data

We start by loading the `Andrews2021` dataset included in the `RoBMA` R package, which contains 39 source rows from the meta-analysis of the effect of household chaos on child executive functions. 
The dataset includes raw correlation coefficients (`ri`), sample sizes (`ni`), study labels (`slab`), and study-level moderators such as the type of executive function assessment (`measure`) and the mean age of the children (`age`). 
The analyses below use the 36 rows with complete effect-size inputs.

```{r}
library(RoBMA)
data("Andrews2021", package = "RoBMA")
head(Andrews2021)[,c("slab", "ri", "ni", "measure", "age")]
```

We use `metafor::escalc()` to compute Fisher's $z$ effect sizes directly from correlation coefficients and sample sizes which we will use for all the subsequent analyses.

```{r escalc}
Andrews2021_z <- metafor::escalc(
  ri = ri, ni = ni,
  measure = "ZCOR", data = Andrews2021
)
```

## Frequentist Meta-Regression

We start by fitting a frequentist meta-regression using the `metafor` package [@metafor].
While @andrews2021examining estimated univariate meta-regressions for each moderator, we directly proceed by analyzing both moderators simultaneously.
The model is estimated on Fisher's $z$ scale and the pooled estimates are transformed to correlation coefficients only for reporting.

```{r}
fit_rma <- metafor::rma(yi, vi, mods = ~ measure + age, data = Andrews2021_z)
fit_rma
```

The results reveal a statistically significant moderation effect of the executive function assessment type on the effect of household chaos on child executive functions (p = 0.013).
To explore the moderation effect further, we estimate the estimated marginal means for the executive function assessment type using the `emmeans` R package [@emmeans] and back-transform the estimates to correlations.

```{r}
emm_rma <- emmeans::emmeans(metafor::emmprep(fit_rma), specs = "measure")
emm_rma
metafor::transf.ztor(as.data.frame(emm_rma)[, c("emmean", "asymp.LCL", "asymp.UCL")])
```

Studies using the informant-completed questionnaires show a stronger effect of household chaos on child executive functions, r = 0.227, 95\% CI [0.160, 0.293], than the direct assessment, r = 0.111, 95\% CI [0.050, 0.172]; both types of studies show statistically significant effects.

The mean age of the children does not significantly moderate the effect (p = 0.600) with the estimated regression coefficient of b = 0.003, 95\% CI [-0.009, 0.016] on Fisher's $z$ scale.
As usual, frequentist inference only allows us to fail to reject the null hypothesis.
Here, we try to overcome this limitation with Bayesian model-averaged meta-regression.


## Bayesian Meta-Regression Specification

Before we proceed with the Bayesian model-averaged meta-regression, we provide a quick overview of the regression model specification. 
In contrast to frequentist meta-regression, we need to specify prior distributions on the regression coefficients, which encode the tested hypotheses about the presence vs. absence of the moderation (specifying different prior distributions corresponds to different hypotheses and results in different conclusions). 
Importantly, the treatment of continuous and categorical covariates differs in the Bayesian model-averaged meta-regression.

### Continuous vs. Categorical Moderators and Prior Distributions

The package automatically standardizes continuous moderators, which makes the moderator prior distributions scale-invariant and keeps the intercept prior distribution centered on the grand mean effect.
This setting can be overridden with `standardize_continuous_predictors = FALSE`.

Categorical moderators use mean-difference contrasts by default for model-averaged models (`set_contrast_factor_predictors = "meandif"`), modeling deviations of each category from the grand mean effect.
The `"meandif"` contrasts achieve label invariance and the prior distribution for the intercept corresponds to the grand mean effect.
Alternatively, `"treatment"` contrasts place prior distributions on differences from the reference level, with the intercept corresponding to the effect in that level.

### Prior Distributions

Note that the effect-size and heterogeneity default prior distributions in the current (4.0) version of the package differ from those used by `RoBMA` 3.6.1.
The current default prior distributions for Bayesian meta-analyses are calibrated to the unit-information standard deviation of the effect-size measure, specified via the `measure` argument (see the [*Prior Distributions*](v01-prior-distributions.html) vignette).

To reproduce the prior specification used in @bartos2023robust, we specify the prior distributions explicitly below.
The old effect-size and heterogeneity prior distributions, expressed on the fitted Fisher's $z$ scale, are specified explicitly.[^prior-scale]

```{r prior-specification}
prior_effect_361 <- prior("normal", list(mean = 0, sd = 0.5))
prior_heterogeneity_361 <- prior("invgamma", list(shape = 1, scale = 0.075))
```

[^prior-scale]: `RoBMA` 3.6.1 used Cohen's $d$ as the prior scale for correlation inputs while fitting Fisher's $z$ internally. The prior distributions above express the old effect-size and heterogeneity scales on Fisher's $z$ using the local approximation $z \approx d / 2$.

## Bayesian Model-Averaged Meta-Regression

We fit the Bayesian model-averaged meta-regression using `BMA()`, which omits publication-bias adjustment.
We specify the Fisher's $z$ effect sizes as `yi` and their sampling variances as `vi`, and set the `measure` argument to `"ZCOR"`, which tells `BMA()` and `RoBMA()` how to interpret the input scale.
The moderators are passed via the `mods` formula, here `~ measure + age`.
We also set `parallel = TRUE` to speed up computation and `seed = 1` for reproducibility.

```{r fit-bma, eval = FALSE}
fit_BMA <- BMA(
  yi = yi, vi = vi, measure = "ZCOR",
  mods = ~ measure + age,
  prior_effect        = prior_effect_361,
  prior_heterogeneity = prior_heterogeneity_361,
  data = Andrews2021_z, parallel = TRUE, seed = 1
)
```

`BMA()` specifies the combination of all models assuming presence vs. absence of the effect, heterogeneity, moderation by `measure`, and moderation by `age`, which corresponds to 2\*2\*2\*2=16 models.
The new version of the package uses product-space parameterization, which greatly improves estimation time by fitting all of the models simultaneously within a single MCMC.

Once the ensemble is estimated, `summary()` reports model-averaging results on the fitted Fisher's $z$ scale.
Dedicated effect-summary helpers handle transformations to the correlation scale.

```{r}
summary(fit_BMA, include_mcmc_diagnostics = FALSE)
pooled_effect(fit_BMA, output_measure = "COR")
```

The summary function produces output with multiple sections.
The first section contains the `Components summary` with the hypothesis test results for the overall effect size and heterogeneity. 
We find overwhelming evidence for both with inclusion Bayes factors (`Inclusion BF`) above 10,000 (the numerical precision of Bayes factors above 100 is usually limited by the number of posterior samples).

The second section contains the `Meta-regression components summary` with the hypothesis test results for the moderators. 
We find weak evidence for the moderation by the executive function assessment type, $\text{BF}_{\text{measure}} = 2.67$. 
Furthermore, we find moderate evidence for the null hypothesis of no moderation by mean age of the children, $\text{BF}_{\text{age}} = 0.122$ (i.e., BF for the null is 1/0.122 = 8.20). 
These findings extend the frequentist meta-regression by disentangling the absence of evidence from the evidence of absence.

The third section contains the `Model-averaged estimates` with the model-averaged estimates for mean effect and between-study heterogeneity.
`pooled_effect()` reports the pooled mean effect on the correlation scale via `output_measure = "COR"`.

The fourth section contains the `Model-averaged meta-regression estimates` with the model-averaged regression coefficient estimates. 
The main difference from the usual frequentist meta-regression output is that the categorical predictors are summarized as a difference from the grand mean for each factor level. 
Here, the `intercept` regression coefficient estimate corresponds to the grand mean effect and the `measure[dif: direct]` regression coefficient estimate of -0.044, 95\% CI [-0.105, 0.000] corresponds to the difference between the direct assessment and the grand mean. 
As such, the results suggest that the effect size in studies using direct assessment is lower in comparison to the grand mean of the studies. The `age` regression coefficient of 0.000, 95\% CI [-0.000, 0.009] corresponds to the increase in the mean effect for a one-unit increase in mean age of children.

Similarly to the frequentist meta-regression, we can use `marginal_means()` and `summary()` to obtain the marginal estimates for each factor level on the correlation scale.
```{r}
BMA_marginal <- marginal_means(fit_BMA, output_measure = "COR")
summary(BMA_marginal)
```
The estimated marginal means are similar to the frequentist results. 
Studies using the informant-completed questionnaires again show a stronger effect of household chaos on child executive functions, r = 0.210, 95\% CI [0.128, 0.289], than the direct assessment, r = 0.124, 95\% CI [0.054, 0.196]. 

The last column summarizes results from a test against a null hypothesis of marginal means equals 0. 
Here, we find extreme evidence that the effect size for studies using the informant-completed questionnaires differs from zero, $\text{BF}_{10} > 10{,}000$, and strong evidence that the effect size for studies using direct assessment differs from zero, $\text{BF}_{10} = 23.6$.
The test is performed using the change from prior to posterior distribution at 0 (i.e., the Savage-Dickey density ratio) assuming the presence of the overall effect or the presence of difference according to the tested factor. 
Because the tests use prior and posterior samples, calculating the Bayes factor can be problematic when the posterior distribution is far from the tested value. 
In such cases, warning messages are printed and $\text{BF}_{10} = \infty$ returned (like here); while the actual Bayes factor is less than infinity, it is still too large to be computed precisely given the posterior samples.

The full model-averaged posterior marginal means distribution can be visualized by the `plot()` method.
```{r fig_BMA, fig.width = 6, fig.height = 4, out.width = "75%", fig.align = "center"}
par(mar = c(4, 4, 1, 4))
plot(BMA_marginal, parameter = "measure", lwd = 2)
```

## Publication-Bias-Adjusted Model-Averaged Meta-Regression

Finally, we adjust the Bayesian model-averaged meta-regression model by fitting the robust Bayesian model-averaged meta-regression.
In contrast to the previous publication-bias-unadjusted ensemble, `RoBMA()` extends the model ensemble with a publication-bias component specified via 6 weight functions and PET-PEESE [@bartos2021no].
The estimation time further increases as the ensemble now contains 144 models.

```{r fit-robma, eval = FALSE}
fit_RoBMA <- RoBMA(
  yi = yi, vi = vi, measure = "ZCOR",
  mods = ~ measure + age,
  prior_effect        = prior_effect_361,
  prior_heterogeneity = prior_heterogeneity_361,
  data = Andrews2021_z, parallel = TRUE, seed = 1
)
```
```{r}
summary(fit_RoBMA, include_mcmc_diagnostics = FALSE)
pooled_effect(fit_RoBMA, output_measure = "COR")
```
All previously described functions for manipulating the fitted model work identically with the publication bias adjusted model. 
As such, we just briefly mention the main differences found after adjusting for publication bias.

RoBMA-reg reveals moderate evidence of publication bias $\text{BF}_{\text{pb}} = 4.26$. 
Furthermore, accounting for publication bias turns the previously found evidence for the overall effect into a weak evidence for the effect $\text{BF}_{10} = 1.25$ and notably reduces the mean effect estimate to r = 0.061, 95\% CI [-0.020, 0.189].

```{r}
RoBMA_marginal <- marginal_means(fit_RoBMA, output_measure = "COR")
summary(RoBMA_marginal)
```

The estimated marginal means now suggest that studies using the informant-completed questionnaires show a much smaller effect of household chaos on child executive functions, r = 0.128, 95\% CI [0.000, 0.116] with moderate evidence in favor of the effect, $\text{BF}_{10} = 6.33$, 
while studies using direct assessment even provide weak evidence against the effect of household chaos on child executive functions, $\text{BF}_{10} = 0.626$, with most likely effect sizes around zero, r = 0.013, 95\% CI [-0.115, 0.175].

A visual summary of the estimated marginal means highlights the considerably wider model-averaged posterior distributions of the marginal means, a consequence of accounting and adjusting for publication bias.

```{r fig_RoBMA, fig.width = 6, fig.height = 4, out.width = "75%", fig.align = "center"}
par(mar = c(4, 4, 1, 4))
plot(RoBMA_marginal, parameter = "measure", lwd = 2)
```

The Bayesian model-averaged meta-regression models are compatible with the remaining custom specification, visualization, and summary functions included in the `RoBMA` R package.
For prior customization, see the [*Prior Distributions*](v01-prior-distributions.html) vignette; for shared meta-regression plotting and diagnostics, see the [*Bayesian Meta-Analysis*](v02-bayesian-meta-analysis.html) vignette.

## References