Package 'RoBMA'

Description

RoBMA provides Bayesian meta-analysis, meta-regression, multilevel meta-analysis, model averaging, and publication-bias adjustment. The main user-facing fitters are RoBMA, BMA, brma, brma.glmm, bselmodel, bPET, and bPEESE.

User guide

See Bartoš et al. (2023), Maier et al. (2023), Bartoš et al. (2022), and Bartoš et al. (2026) for the RoBMA methodology. Use vignette(package = "RoBMA") to list installed vignettes.

Author(s)

Frantisek Bartos [email protected]

References

Bartoš F, Maier M, Quintana DS, Wagenmakers E (2022). “Adjusting for publication bias in JASP and R — Selection models, PET-PEESE, and robust Bayesian meta-analysis.” Advances in Methods and Practices in Psychological Science, 5(3), 1–19. doi:10.1177/25152459221109259.

Bartoš F, Maier M, Wagenmakers E (2026). “Robust Bayesian multilevel meta-analysis: Adjusting for publication bias in the presence of dependent effect sizes.” Behavior Research Methods. doi:10.31234/osf.io/9tgp2_v2. Preprint available at https://doi.org/10.31234/osf.io/9tgp2_v2.

Bartoš F, Maier M, Wagenmakers E, Doucouliagos H, Stanley TD (2023). “Robust Bayesian meta-analysis: Model-averaging across complementary publication bias adjustment methods.” Research Synthesis Methods, 14(1), 99–116. doi:10.1002/jrsm.1594.

Maier M, Bartoš F, Wagenmakers E (2023). “Robust Bayesian Meta-Analysis: Addressing publication bias with model-averaging.” Psychological Methods, 28(1), 107–122. doi:10.1037/met0000405.

See Also

Useful links:

• https://fbartos.github.io/RoBMA/

• Report bugs at https://github.com/FBartos/RoBMA/issues

Description

Compute approximate leave-one-out cross-validation (LOO-CV) using Pareto smoothed importance sampling (PSIS) for brma model objects and store the result in the object.

Usage

## S3 method for class 'brma'
add_loo(object, unit = "estimate", r_eff = NULL, parallel = FALSE, ...)

Arguments

Details

With unit = "estimate", LOO-CV is computed with one contribution per effect-size estimate. For binomial and Poisson models, each pair of counts (ai/ci or x1i/x2i) that defines a single effect size estimate is treated as one contribution.

With unit = "cluster", LOO-CV is computed with one joint contribution per cluster. For unweighted normal models without selection this uses the analytic cluster block covariance. Selection, data-weighted normal, and GLMM models integrate the held-out cluster effect with Gauss-Hermite quadrature.

For selection models, the LOO evaluates the weighted likelihood, conditioning on the posterior omega samples.

The PSIS object is essential for model comparison via loo_compare and is automatically saved in the loo result. RoBMA stores target metadata so comparisons can reject mismatched data, unit, or conditioning-depth targets.

Important for model comparison: When comparing models via loo_compare, the selection is based on expected out-of-sample predictive performance. This evaluates how well models predict new observations, not how well they fit the observed data.

Value

The brma object with the LOO result stored in object[["loo"]][[unit]].

References

(Vehtari et al. 2017) (Vehtari et al. 2024)

See Also

loo.brma, loo, loo_compare, pareto_k_ids

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- bPET(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")

  fit <- add_loo(fit)
  loo_fit <- loo(fit)
  print(loo_fit)
  plot(loo_fit)
}

## End(Not run)

Description

Compute the marginal likelihood of a brma model using bridge sampling and store the result in the object.

Usage

## S3 method for class 'brma'
add_marglik(object, ...)

Arguments

Details

The marginal likelihood is computed using the bridgesampling package via the BayesTools::JAGS_bridgesampling wrapper. The result is stored in the object and can be extracted using bridge_sampler.brma.

Product-space model-averaging objects (BMA.norm, BMA.glmm, and RoBMA) do not expose a bridge-sampling marginal likelihood; use predictive comparison methods such as loo.brma instead.

Value

The brma object with the marginal likelihood result stored in object[["marglik"]].

See Also

bridge_sampler.brma, logml.brma, bf.brma, post_prob.brma

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- brma(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")

  fit <- add_marglik(fit)

  bridge <- bridge_sampler(fit)
  print(bridge)

  logml(fit)
}

## End(Not run)

Description

Compute the Widely Applicable Information Criterion (WAIC) for brma model objects and store the result in the object.

Usage

## S3 method for class 'brma'
add_waic(object, unit = "estimate", ...)

Arguments

Details

WAIC is an alternative to LOO-CV for estimating out-of-sample predictive accuracy. Like LOO, it evaluates expected predictive performance for new observations.

In most cases, LOO-PSIS (via add_loo) is preferred over WAIC because it provides better estimates and includes diagnostics (Pareto k values) that indicate when the approximation may be unreliable.

Value

The brma object with the WAIC result stored in object[["waic"]][[unit]].

See Also

waic.brma, add_loo, waic

Description

The data set contains correlation coefficients, sample sizes, and labels for 23 experimental studies focusing on the effect of violent video games on aggressive behavior.

Usage

Anderson2010

Format

A data.frame with 3 columns and 23 observations:

r

Correlation coefficient.

n

Sample size.

name

Study label.

Source

https://github.com/Joe-Hilgard/Anderson-meta

References

Anderson CA, Shibuya A, Ihori N, Swing EL, Bushman BJ, Sakamoto A, Rothstein HR, Saleem M (2010). “Violent video game effects on aggression, empathy, and prosocial behavior in Eastern and Western countries: A meta-analytic review.” Psychological Bulletin, 136(2), 151–173. doi:10.1037/a0018251.

Description

The data set contains raw correlation coefficients and sampling information using metafor-compatible column names (ri, vi, sei, ni, and slab), together with study-level moderators from the original data set. Three source rows have missing effect-size inputs and are omitted automatically by model-fitting functions that require complete outcomes. The original data set assessed the effect of household chaos on child executive functions (Andrews et al. 2021) and was used as an example in Bartoš et al. (2025).

Usage

Andrews2021

Format

A data.frame with 15 columns and 39 observations:

study_id

Source row identifier.

slab

Study label.

ri

Raw correlation coefficient.

vi

Sampling variance of the raw correlation coefficient.

sei

Sampling standard error of the raw correlation coefficient.

ni

Sample size.

year

Publication year.

percent_female

Percentage of female children in the sample.

age

Mean age of the children in years.

assessment_interval_months

Time between household-chaos and executive-function assessments in months.

dissertation

Whether the study was a dissertation.

percent_minority

Percentage of minority participants in the sample.

percent_low_parental_education

Percentage of parents with high school, GED, or lower education.

hc_dimension

Household-chaos dimension using the source coding.

measure

Executive-function assessment type, direct or informant.

References

Andrews K, Atkinson L, Harris M, Gonzalez A (2021). “Examining the effects of household chaos on child executive functions: A meta-analysis.” Psychological Bulletin, 147(1), 16–32. doi:10.1037/bul0000311.

Bartoš F, Maier M, Stanley TD, Wagenmakers E (2025). “Robust Bayesian meta-regression: Model-averaged moderation analysis in the presence of publication bias.” Psychological Methods. doi:10.1037/met0000737.

Description

Provides an interface to the posterior package for brma objects. These functions convert the MCMC samples from a fitted brma model to various draws formats supported by the posterior package, enabling the use of posterior's rich set of diagnostics and summary functions.

Usage

as_draws(x, ...)

as_draws_array(x, ...)

as_draws_df(x, ...)

as_draws_list(x, ...)

as_draws_matrix(x, ...)

as_draws_rvars(x, ...)

## Default S3 method:
as_draws(x, ...)

## Default S3 method:
as_draws_array(x, ...)

## Default S3 method:
as_draws_df(x, ...)

## Default S3 method:
as_draws_list(x, ...)

## Default S3 method:
as_draws_matrix(x, ...)

## Default S3 method:
as_draws_rvars(x, ...)

## S3 method for class 'brma'
as_draws(x, ...)

## S3 method for class 'brma'
as_draws_array(x, ...)

## S3 method for class 'brma'
as_draws_df(x, ...)

## S3 method for class 'brma'
as_draws_list(x, ...)

## S3 method for class 'brma'
as_draws_matrix(x, ...)

## S3 method for class 'brma'
as_draws_rvars(x, ...)

Arguments

Details

These functions are S3 generics. Their default methods forward to the corresponding posterior generics so attaching RoBMA preserves the usual posterior behavior for non-brma objects.

The following conversion functions are available:

• as_draws: converts to the default draws format

• as_draws_array: converts to a 3-D array (iteration x chain x variable)

• as_draws_df: converts to a data frame with columns for iterations, chains, and variables

• as_draws_list: converts to a list of lists

• as_draws_matrix: converts to a 2-D matrix (draw x variable)

• as_draws_rvars: converts to random variable objects

These methods require the posterior package to be installed. For brma methods, conversion is performed by first extracting the MCMC samples as a mcmc.list object and then using the corresponding posterior conversion function. brma_samples objects have separate methods documented at as_draws.brma_samples.

Value

An object of the corresponding posterior draws class.

See Also

draws, brma, as_draws.brma_samples

Description

Provides an interface to the posterior package for brma_samples objects. These functions convert the posterior samples to various draws formats supported by the posterior package.

Usage

## S3 method for class 'brma_samples'
as_draws(x, ...)

## S3 method for class 'brma_samples'
as_draws_array(x, ...)

## S3 method for class 'brma_samples'
as_draws_df(x, ...)

## S3 method for class 'brma_samples'
as_draws_list(x, ...)

## S3 method for class 'brma_samples'
as_draws_matrix(x, ...)

## S3 method for class 'brma_samples'
as_draws_rvars(x, ...)

Arguments

Details

The conversion reconstructs the MCMC chain structure from the stored nchains and niter attributes. The samples are assumed to be ordered with chains concatenated (i.e., all iterations from chain 1, then all from chain 2, etc.). Conditional RoBMA samples are intentionally stored as one flattened chain because conditioning subsets posterior rows across chains.

Value

An object of the corresponding posterior draws class.

See Also

draws, as_draws.brma

Description

Transforms an estimated brma model into a zplot object that can be summarized and plotted to assess replicability.

Usage

## S3 method for class 'brma'
as_zplot(
  object,
  significance_level = stats::qnorm(0.975),
  max_samples = 10000,
  ...
)

Arguments

Details

Zplot analysis estimates the Expected Discovery Rate (EDR), the posterior mean probability that an exact replication would be statistically significant at the supplied threshold. This provides insight into the replicability of findings in a literature.

The implementation uses extrapolation mode by default, which removes selection-model weights when evaluating the model-implied z-value density. PET/PEESE regression terms remain part of the fitted location model. Zplot diagnostics are available only for normal outcome models. GLMM objects are rejected because their raw likelihood is on a count scale while zplot diagnostics require observed effect-size z-statistics with standard errors.

The resulting object retains all original brma properties while adding zplot results, enabling both standard meta-analytic summaries and zplot diagnostics on the same object.

Value

The input object with added class "zplot_brma" and a new zplot list component containing:

estimates

a list with posterior samples for EDR and missing-study weights

data

a list with significance_level, observed z-statistics, N_significant, and N_observed

See Also

summary.zplot_brma(), plot.zplot_brma(), hist.zplot_brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- bPET(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")

  zfit <- as_zplot(fit)
  summary(zfit)
  plot(zfit)
}

## End(Not run)

Description

Converts a brma_samples object to a plain matrix, removing all brma_samples-specific attributes.

Usage

## S3 method for class 'brma_samples'
as.matrix(x, ...)

Arguments

Value

A plain matrix of posterior samples.

Description

The data set contains Cohen's d effect sizes, standard errors, and labels for 9 experimental studies of precognition from the infamous Bem (2011) as analyzed in his later meta-analysis (Bem et al. 2011).

Usage

Bem2011

Format

A data.frame with 3 columns and 9 observations:

d

Cohen's d effect size.

se

Standard error of d.

study

Study label.

References

Bem DJ (2011). “Feeling the future: Experimental evidence for anomalous retroactive influences on cognition and affect.” Journal of Personality and Social Psychology, 100(3), 407–425. doi:10.1037/a0021524.

Bem DJ, Utts J, Johnson WO (2011). “Must psychologists change the way they analyze their data?” Journal of Personality and Social Psychology, 101(4), 716–719. doi:10.1037/a0024777.

Description

Compute the Bayes factor comparing two brma models.

Usage

## S3 method for class 'brma'
bf(x1, x2, log = FALSE, ...)

## S3 method for class 'brma'
bayes_factor(x1, x2, log = FALSE, ...)

Arguments

Details

Computes the Bayes factor in favor of the model x1 over the model x2. The marginal likelihoods must first be computed using add_marglik. Both models must be fitted to the same outcome target/data, including outcome type and, when present, weights and cluster identifiers.

Value

A list of class "bf_default" with components:

• bf: (scalar) value of the Bayes factor in favor of x1 over x2.

• log: Boolean indicating whether bf corresponds to the log Bayes factor.

See Also

add_marglik, bridge_sampler.brma, post_prob.brma, logml.brma

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit1 <- brma(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")
  fit2 <- brma(
    yi           = yi,
    vi           = vi,
    data         = dat.lehmann2018,
    measure      = "SMD",
    prior_effect = FALSE
  )

  fit1 <- add_marglik(fit1)
  fit2 <- add_marglik(fit2)

  bf(fit1, fit2)
}

## End(Not run)

Description

Computes the estimated true effects (theta) from a fitted model. These correspond to Best Linear Unbiased Predictions (BLUPs) or empirical Bayes estimates.

Usage

blup(object, ...)

Arguments

Value

Method-specific return value, typically a summary table or posterior samples of BLUP or empirical-Bayes true-effect summaries.

See Also

predict.brma()

Description

Computes the estimated true effects (theta) for a fitted brma object. These correspond to Best Linear Unbiased Predictions (BLUPs) or empirical Bayes estimates.

Usage

## S3 method for class 'brma'
blup(
  object,
  bias_adjusted = FALSE,
  output_measure = NULL,
  transform = NULL,
  probs = c(0.025, 0.975),
  ...
)

Arguments

Details

This function is a convenience wrapper around predict.brma(..., type = "effect", newdata = NULL).

For unweighted two-level normal models, true effects are computed using empirical Bayes shrinkage:

$\theta_i = \lambda_i \cdot y_i + (1 - \lambda_i) \cdot \mu_i$

where $\lambda_i = \tau^2 / (\tau^2 + se_i^2)$. With likelihood weights, $se_i^2$ is replaced by the weighted sampling variance $se_i^2 / w_i$.

For GLMM models (binomial, Poisson), the estimate-level random effects are extracted directly from the posterior samples.

For multilevel (3-level) normal models, cluster-level effects are estimated jointly within cluster blocks and estimate-level effects are then shrunk conditional on those cluster effects.

Value

A brma_samples object containing posterior draws of BLUP or empirical-Bayes true-effect summaries with one column per estimate. For existing normal data, these are conditional BLUP means, not simulated latent-effect draws. When printed, displays a summary table. Use summary() to obtain the summary table directly. The samples can be converted to posterior draws formats using as_draws().

See Also

predict.brma(), pooled_effect(), pooled_heterogeneity(), true_effects()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- brma(
    yi      = yi,
    vi      = vi,
    data    = dat.lehmann2018,
    measure = "SMD",
    seed    = 1,
    silent  = TRUE
  )

  blup(fit)
}

## End(Not run)

Description

Fits Bayesian model-averaged meta-analytic models without publication-bias adjustment. BMA() is an alias for BMA.norm().

Usage

BMA(
  yi,
  vi,
  sei,
  weights,
  ni,
  mods,
  scale,
  cluster,
  data,
  slab,
  subset,
  measure,
  prior_effect,
  prior_heterogeneity,
  prior_mods,
  prior_scale,
  prior_heterogeneity_allocation,
  prior_effect_null,
  prior_heterogeneity_null,
  prior_mods_null,
  prior_scale_null,
  prior_heterogeneity_allocation_null,
  standardize_continuous_predictors = TRUE,
  set_contrast_factor_predictors = "meandif",
  prior_unit_information_sd,
  rescale_priors = 1,
  prior_informed_field,
  prior_informed_subfield,
  sample = 5000,
  burnin = 2000,
  adapt = 500,
  chains = 3,
  thin = 1,
  parallel = FALSE,
  autofit = FALSE,
  autofit_control = set_autofit_control(),
  convergence_checks = set_convergence_checks(),
  seed = NULL,
  silent = TRUE,
  ...
)

Arguments

Details

BMA.norm (and its alias BMA) provides a simplified interface for Bayesian model-averaged meta-analysis when publication bias adjustment is not needed. It uses the same product-space mixture-prior machinery as RoBMA() but omits selection models and PET-PEESE bias adjustment.

For publication bias adjusted meta-analysis, use RoBMA directly.

Value

A fitted object of class c("BMA.norm", "RoBMA", "brma"). The object contains checked data, checked mixture priors, the JAGS fit, cached summary, and cached coefficients.

See Also

RoBMA() brma() summary.brma() plot.brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")

  fit <- BMA(
    yi      = yi,
    vi      = vi,
    mods    = ~ Preregistered,
    data    = dat.lehmann2018,
    measure = "SMD",
    seed    = 1,
    silent  = TRUE
  )

  summary(fit)
}

## End(Not run)

Description

Function for fitting Bayesian model-averaged meta-analytic models directly to binary or count data using a generalized linear mixed model (GLMM) framework. Unlike RoBMA, this function does not adjust for publication bias, as weight function and regression-based bias adjustment methods are not available for GLMM outcomes.

Usage

BMA.glmm(
  ai,
  bi,
  ci,
  di,
  n1i,
  n2i,
  x1i,
  x2i,
  t1i,
  t2i,
  weights,
  mods,
  scale,
  cluster,
  data,
  slab,
  subset,
  measure = "OR",
  prior_effect,
  prior_heterogeneity,
  prior_mods,
  prior_scale,
  prior_heterogeneity_allocation,
  prior_baserate,
  prior_lograte,
  prior_effect_null,
  prior_heterogeneity_null,
  prior_mods_null,
  prior_scale_null,
  prior_heterogeneity_allocation_null,
  standardize_continuous_predictors = TRUE,
  set_contrast_factor_predictors = "treatment",
  prior_unit_information_sd,
  rescale_priors = 1,
  prior_informed_field,
  prior_informed_subfield,
  sample = 5000,
  burnin = 2000,
  adapt = 500,
  chains = 3,
  thin = 1,
  parallel = FALSE,
  autofit = FALSE,
  autofit_control = set_autofit_control(),
  convergence_checks = set_convergence_checks(),
  seed = NULL,
  silent = TRUE,
  ...
)

Arguments

Details

BMA.glmm combines the data input style of brma.glmm with the mixture prior specification of RoBMA for Bayesian model-averaging.

Model for odds ratios (measure = "OR") uses a binomial-normal model as described in Jackson et al. (2018).

Model for incidence rate ratios (measure = "IRR") uses a Poisson-normal model as described in Bagos and Nikolopoulos (2009).

When weights are supplied, they are treated as likelihood weights on the paired two-arm study contribution.

Value

A fitted object of class c("BMA.glmm", "RoBMA", "brma.glmm", "brma"). The object contains checked data, checked mixture priors, the JAGS fit, cached summary, and cached coefficients.

See Also

brma.glmm() RoBMA() summary.brma() plot.brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.bcg, package = "metadat")

  fit <- BMA.glmm(
    ai      = tpos,
    bi      = tneg,
    ci      = cpos,
    di      = cneg,
    data    = dat.bcg,
    measure = "OR",
    seed    = 1,
    silent  = TRUE
  )

  summary(fit)
}

## End(Not run)

Description

Function for fitting random-effects, meta-regression, multilevel, and location-scale meta-analytic PEESE models.

Usage

bPEESE(
  yi,
  vi,
  sei,
  weights,
  ni,
  mods,
  scale,
  cluster,
  data,
  slab,
  subset,
  measure,
  prior_effect,
  prior_heterogeneity,
  prior_mods,
  prior_scale,
  prior_heterogeneity_allocation,
  prior_bias,
  standardize_continuous_predictors = TRUE,
  set_contrast_factor_predictors = "treatment",
  prior_unit_information_sd,
  rescale_priors = 1,
  prior_informed_field,
  prior_informed_subfield,
  effect_direction = "detect",
  sample = 5000,
  burnin = 2000,
  adapt = 500,
  chains = 3,
  thin = 1,
  parallel = FALSE,
  autofit = FALSE,
  autofit_control = set_autofit_control(),
  convergence_checks = set_convergence_checks(),
  seed = NULL,
  silent,
  ...
)

Arguments

Value

A fitted object of class c("bPEESE", "brma") containing a single PEESE publication-bias model fit.

See Also

publication_bias_prior_specification, RoBMA(), bPET(), bselmodel(), summary.brma(), funnel.brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")

  fit <- bPEESE(
    yi      = yi,
    vi      = vi,
    data    = dat.lehmann2018,
    measure = "SMD",
    seed    = 1,
    silent  = TRUE
  )

  summary(fit)
  funnel(fit)
}

## End(Not run)

Description

Function for fitting random-effects, meta-regression, multilevel, and location-scale meta-analytic PET models.

Usage

bPET(
  yi,
  vi,
  sei,
  weights,
  ni,
  mods,
  scale,
  cluster,
  data,
  slab,
  subset,
  measure,
  prior_effect,
  prior_heterogeneity,
  prior_mods,
  prior_scale,
  prior_heterogeneity_allocation,
  prior_bias,
  standardize_continuous_predictors = TRUE,
  set_contrast_factor_predictors = "treatment",
  prior_unit_information_sd,
  rescale_priors = 1,
  prior_informed_field,
  prior_informed_subfield,
  effect_direction = "detect",
  sample = 5000,
  burnin = 2000,
  adapt = 500,
  chains = 3,
  thin = 1,
  parallel = FALSE,
  autofit = FALSE,
  autofit_control = set_autofit_control(),
  convergence_checks = set_convergence_checks(),
  seed = NULL,
  silent,
  ...
)

Arguments

Value

A fitted object of class c("bPET", "brma") containing a single PET publication-bias model fit.

See Also

publication_bias_prior_specification, RoBMA(), bPEESE(), bselmodel(), summary.brma(), funnel.brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")

  fit <- bPET(
    yi      = yi,
    vi      = vi,
    data    = dat.lehmann2018,
    measure = "SMD",
    seed    = 1,
    silent  = TRUE
  )

  summary(fit)
  funnel(fit)
}

## End(Not run)

Description

Extract the marginal likelihood bridge sampling object from a brma model. The marginal likelihood must first be computed using add_marglik.

Usage

## S3 method for class 'brma'
bridge_sampler(samples, ...)

Arguments

Details

This function extracts the bridge sampling object that was previously computed and stored using add_marglik. If the marginal likelihood has not been computed, an error is thrown. Product-space model-averaging objects (BMA.norm, BMA.glmm, and RoBMA) do not expose bridge-sampling marginal likelihoods.

The returned object can be used for Bayesian model comparison via bf and post_prob.

Value

An object of class "bridge" as returned by bridge_sampler.

See Also

add_marglik, bridge_sampler, logml.brma, bf.brma, post_prob.brma

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- brma(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")

  fit <- add_marglik(fit)

  bridge <- bridge_sampler(fit)
  print(bridge)
}

## End(Not run)

Description

Function for fitting normal-likelihood/effect-size random-effects, meta-regression, multilevel, and location-scale meta-analytic models. brma.norm() is an alias for brma(); raw-count GLMM models use brma.glmm().

Usage

brma(
  yi,
  vi,
  sei,
  weights,
  ni,
  mods,
  scale,
  cluster,
  data,
  slab,
  subset,
  measure,
  prior_effect,
  prior_heterogeneity,
  prior_mods,
  prior_scale,
  prior_heterogeneity_allocation,
  standardize_continuous_predictors = TRUE,
  set_contrast_factor_predictors = "treatment",
  prior_unit_information_sd,
  rescale_priors = 1,
  prior_informed_field,
  prior_informed_subfield,
  sample = 5000,
  burnin = 2000,
  adapt = 500,
  chains = 3,
  thin = 1,
  parallel = FALSE,
  autofit = FALSE,
  autofit_control = set_autofit_control(),
  convergence_checks = set_convergence_checks(),
  seed = NULL,
  silent,
  ...
)

Arguments

Details

Prior distributions

Prior distributions must be specified for all model parameters. This typically includes the pooled effect $\mu$ and between-study heterogeneity $\tau$. In the case of meta-regression, the pooled effect $\mu$ corresponds to the intercept, and additional prior distributions for the regression coefficients are required. In the case of a location-scale model, the between-study heterogeneity corresponds to the intercept of the scale regression, and additional prior distributions for the scale regression coefficients are required.

There are several ways to specify the prior distributions:

1. via a standardized effect size measure with known unit information standard deviation,

2. by estimating unit information standard deviation using sample sizes ni,

3. by manually setting prior_unit_information_sd,

4. by specifying informed empirical prior distributions via prior_informed_field and prior_informed_subfield,

5. or via fully custom specification using the prior_effect, prior_heterogeneity, prior_mods, prior_scale, and prior_heterogeneity_allocation arguments.

In all cases, the prior behavior can be further modified by the rescale_priors, standardize_continuous_predictors, and set_contrast_factor_predictors arguments.

(1) Specifying prior distributions via standardized effect size measures with known

unit information standard deviation This is the easiest way to specify prior distributions. The width of prior distributions is based on a fraction of the known unit information standard deviation (UISD) (Röver et al. 2021). The default prior distributions for the parameters are set as follows:

The heterogeneity moderation parameters are multiplicative, as such they are independent of UISD.

The default fraction of the UISD can be changed via RoBMA.options() using one of the following arguments: "default_UISD.effect", "default_UISD.heterogeneity", "default_UISD.mods", "default_UISD.scale".

The known UISD for standardized effect size measures (measure) are set as follows:

See Chapter 2.4 in Spiegelhalter et al. (2004) and Chapter 1 in Grieve (2022).

(2) Estimating unit information standard deviation using sample sizes

When effect sizes are on a non-standardized scale (measure = "GEN") or use a standardized effect size without known UISD, the UISD can be estimated from sample sizes (ni) and standard errors (sei) following Equation 6 in Röver et al. (2021). The estimated UISD is then used to scale the default prior distributions as described in section (1).

Note that the known UISD for standardized effect size measures (section (1)) is used if available, even when ni is provided.

(3) Manually setting unit information standard deviation

Alternatively, the UISD can be specified directly via the prior_unit_information_sd argument. This is useful when the appropriate scale for prior distributions is known a priori or when multiple analyses are to be performed on subsets of the same data (re-estimating UISD on different subsets of the data can lead to slightly different prior distributions for different subsets; see estimate_unit_information_sd()). The specified prior_unit_information_sd is then used to scale the default prior distributions as described in section (1).

Note that the manually specified prior_unit_information_sd takes precedence over the estimated UISD from ni (section (2)) and the known UISD from measure (section (1)). It cannot be combined with prior_informed_field.

(4) Specifying informed empirical prior distributions

Informed prior distributions can be specified via the prior_informed_field and prior_informed_subfield arguments. Currently, only prior_informed_field = "medicine" with subfields defined in BayesTools::prior_informed_medicine_names is supported, which uses empirically derived prior distributions from medical meta-analyses as described in Bartoš et al. (2021) and Bartoš et al. (2023).

When prior_informed_field = "medicine", the default prior_informed_subfield is "Cochrane" (i.e., using the whole CDSR database as a reference). The informed prior distributions are available for the following effect size measures:

Note that informed prior distributions are only available for the effect size ($\mu$) and heterogeneity ($\tau$) parameters. For effect moderation (prior_mods), the informed effect prior is scaled by a factor of RoBMA.get_option("default_informed_priors.mods"). For heterogeneity moderation (prior_scale), a normal prior with standard deviation specified by RoBMA.get_option("default_informed_priors.scale") is used.

(5) Fully custom prior distributions

Prior distributions can be fully customized by directly specifying the prior_effect, prior_heterogeneity, prior_mods, prior_scale, and prior_heterogeneity_allocation arguments. These should be prior distribution objects created via BayesTools::prior() or related functions (e.g., prior_factor()).

Rescaling prior distributions

The rescale_priors argument allows rescaling supported prior distributions by a multiplicative factor. For example, rescale_priors = 2 doubles the standard deviations/scales of normal, Cauchy, t, and inverse-gamma prior distributions, making them more diffuse. Point and none priors are unchanged.

Handling of continuous and factor predictors

When standardize_continuous_predictors = TRUE, continuous moderator and scale predictors are internally standardized before fitting. Reported summaries are transformed back to the original predictor scale by default; use standardized_coefficients = TRUE in summary() or related methods to inspect coefficients on the standardized scale.

Factor predictors use the contrast specified by set_contrast_factor_predictors. The default "treatment" follows standard treatment coding. Model-averaging functions default to "meandif" so inclusion priors for factor levels are centered on deviations from the grand mean.

Value

A fitted object of class c("brma.norm", "brma"). The object contains checked data, checked priors, the JAGS fit, cached summary, and cached coefficients. If the corresponding package options are enabled, it can also contain cached LOO, WAIC, or marginal likelihood results. The advanced internal only_data = TRUE and only_priors = TRUE paths return partially constructed objects.

References

Bartoš F, Gronau QF, Timmers B, Otte WM, Ly A, Wagenmakers E (2021). “Bayesian model-averaged meta-analysis in medicine.” Statistics in Medicine, 40(30), 6743–6761. doi:10.1002/sim.9170.

Bartoš F, Otte WM, Gronau QF, Timmers B, Ly A, Wagenmakers E (2023). “Empirical prior distributions for Bayesian meta-analyses of binary and time-to-event outcomes.” doi:10.48550/arXiv.2306.11468. Preprint available at https://doi.org/10.48550/arXiv.2306.11468.

Grieve AP (2022). Hybrid frequentist/Bayesian power and Bayesian power in planning clinical trials. Chapman and Hall/CRC.

Röver C, Bender R, Dias S, Schmid CH, Schmidli H, Sturtz S, Weber S, Friede T (2021). “On weakly informative prior distributions for the heterogeneity parameter in Bayesian random-effects meta-analysis.” Research Synthesis Methods, 12(4), 448–474. doi:10.1002/jrsm.1475.

Spiegelhalter DJ, Abrams KR, Myles JP (2004). Bayesian approaches to clinical trials and health-care evaluation. John Wiley and Sons. doi:10.1002/0470092602.

See Also

RoBMA(), BMA(), brma.glmm(), summary.brma(), plot.brma(), predict.brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE) &&
    requireNamespace("metafor", quietly = TRUE)) {
  data(dat.bcg, package = "metadat")
  dat <- metafor::escalc(
    measure = "RR",
    ai      = tpos,
    bi      = tneg,
    ci      = cpos,
    di      = cneg,
    data    = dat.bcg
  )

  fit <- brma(
    yi      = yi,
    vi      = vi,
    mods    = ~ ablat + year,
    data    = dat,
    measure = "RR",
    seed    = 1,
    silent  = TRUE
  )

  summary(fit)
  predict(fit, type = "terms")
}

## End(Not run)

Description

Function for fitting random-effects, meta-regression, multilevel, and location-scale meta-analytic models directly to either binary or count data.

Usage

brma.glmm(
  ai,
  bi,
  ci,
  di,
  n1i,
  n2i,
  x1i,
  x2i,
  t1i,
  t2i,
  weights,
  mods,
  scale,
  cluster,
  data,
  slab,
  subset,
  measure = "OR",
  prior_effect,
  prior_heterogeneity,
  prior_mods,
  prior_scale,
  prior_heterogeneity_allocation,
  prior_baserate,
  prior_lograte,
  standardize_continuous_predictors = TRUE,
  set_contrast_factor_predictors = "treatment",
  prior_unit_information_sd,
  rescale_priors = 1,
  prior_informed_field,
  prior_informed_subfield,
  sample = 5000,
  burnin = 2000,
  adapt = 500,
  chains = 3,
  thin = 1,
  parallel = FALSE,
  autofit = FALSE,
  autofit_control = set_autofit_control(),
  convergence_checks = set_convergence_checks(),
  seed = NULL,
  silent,
  ...
)

Arguments

Details

Model for odds ratios (measure = "OR") corresponds to Model 4 described in Jackson et al. (2018). logit(pi[i]) is the study-specific midpoint of the two arm logits. prior_baserate defines the estimate-specific prior distribution on pi[i].

Model for incidence rate ratios (measure = "IRR") corresponds to Bagos and Nikolopoulos (2009). phi[i] is the study-specific midpoint of the two arm log incidence rates. prior_lograte defines the estimate-specific prior distribution on phi[i]. If unspecified, a unit-information prior is based on the data and used independently for each estimate.

When weights are supplied, they are treated as likelihood weights on the paired two-arm study contribution.

Value

A fitted object of class c("brma.glmm", "brma"). The object contains checked data, checked priors, the JAGS fit, cached summary, and cached coefficients. If the corresponding package options are enabled, it can also contain cached LOO, WAIC, or marginal likelihood results.

See Also

brma(), BMA.glmm(), summary.brma(), predict.brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.bcg, package = "metadat")

  fit <- brma.glmm(
    ai      = tpos,
    bi      = tneg,
    ci      = cpos,
    di      = cneg,
    mods    = ~ alloc,
    data    = dat.bcg,
    measure = "OR",
    seed    = 1,
    silent  = TRUE
  )

  summary(fit)
}

## End(Not run)

Description

Function for fitting random-effects, meta-regression, multilevel, and location-scale meta-analytic selection models.

Usage

bselmodel(
  yi,
  vi,
  sei,
  weights,
  ni,
  mods,
  scale,
  cluster,
  data,
  slab,
  subset,
  measure,
  prior_effect,
  prior_heterogeneity,
  prior_mods,
  prior_scale,
  prior_heterogeneity_allocation,
  prior_bias,
  standardize_continuous_predictors = TRUE,
  set_contrast_factor_predictors = "treatment",
  prior_unit_information_sd,
  rescale_priors = 1,
  prior_informed_field,
  prior_informed_subfield,
  effect_direction = "detect",
  steps,
  sample = 5000,
  burnin = 2000,
  adapt = 500,
  chains = 3,
  thin = 1,
  parallel = FALSE,
  autofit = FALSE,
  autofit_control = set_autofit_control(),
  convergence_checks = set_convergence_checks(),
  seed = NULL,
  silent,
  ...
)

Arguments

Details

bselmodel() is a normal/effect-size selection-model constructor. Custom prior_bias can be a weightfunction prior or a supported BayesTools selection-kernel prior; p-hacking kernels are not supported in active RoBMA.

Value

A fitted object of class c("bselmodel", "brma") containing a single Bayesian selection model fit.

See Also

publication_bias_prior_specification, RoBMA(), bPET(), bPEESE(), summary.brma(), funnel.brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")

  fit <- bselmodel(
    yi      = yi,
    vi      = vi,
    data    = dat.lehmann2018,
    measure = "SMD",
    steps   = 0.025,
    seed    = 1,
    silent  = TRUE
  )

  summary(fit)
  funnel(fit)
}

## End(Not run)

Description

Check Pareto k diagnostics for a brma model object and warn if any values are high.

Usage

## S3 method for class 'brma'
check_loo(object, unit = "estimate", ...)

Arguments

Details

LOO must first be computed with object <- add_loo(object, unit = unit). The method warns when any Pareto $k$ diagnostic is greater than 0.7.

Value

NULL (throws warning if diagnostics are unreliable).

Description

Extract model coefficients (posterior means) from a fitted brma object.

Usage

## S3 method for class 'brma'
coef(object, ...)

Arguments

Value

A named numeric vector of posterior mean coefficients.

See Also

summary.brma()

Description

BayesTools provides several contrast matrix functions for Bayesian factor analysis. These functions create different types of contrast matrices that can be used with factor variables in Bayesian models.

Usage

contr.orthonormal(n, contrasts = TRUE)

contr.meandif(n, contrasts = TRUE)

contr.independent(n, contrasts = TRUE)

Arguments

Details

The package includes the following contrast functions:

contr.orthonormal

Return a matrix of orthonormal contrasts. Code is based on stanova::contr.bayes and corresponding to description by Rouder et al. (2012). Returns a matrix with n rows and k columns, with k = n - 1 if contrasts = TRUE and k = n if contrasts = FALSE.

contr.meandif

Return a matrix of mean difference contrasts. This is an adjustment to the contr.orthonormal that ascertains that the prior distributions on difference between the grand mean and factor level are identical independent of the number of factor levels (which does not hold for the orthonormal contrast). Furthermore, the contrast is re-scaled so the specified prior distribution exactly corresponds to the prior distribution on difference between each factor level and the grand mean – this is approximately twice the scale of contr.orthonormal. Returns a matrix with n rows and k columns, with k = n - 1 if contrasts = TRUE and k = n if contrasts = FALSE.

contr.independent

Return a matrix of independent contrasts – a level for each term. Returns a matrix with n rows and k columns, with k = n if contrasts = TRUE and k = n if contrasts = FALSE.

References

Rouder JN, Morey RD, Speckman PL, Province JM (2012). “Default Bayes factors for ANOVA designs.” Journal of Mathematical Psychology, 56(5), 356–374. doi:10.1016/j.jmp.2012.08.001.

Examples

# Orthonormal contrasts
contr.orthonormal(c(1, 2))
contr.orthonormal(c(1, 2, 3))

# Mean difference contrasts
contr.meandif(c(1, 2))
contr.meandif(c(1, 2, 3))

# Independent contrasts
contr.independent(c(1, 2))
contr.independent(c(1, 2, 3))

Description

Computes Cook's distance for a fitted brma object. Cook's distance measures the aggregate influence of each observation on the model coefficients.

Usage

## S3 method for class 'brma'
cooks.distance(model, ...)

Arguments

Details

Cook's distance is computed as a PSIS leave-one-out deletion diagnostic. For each observation $i$, normalized PSIS weights estimate the fitted values under the leave-one-out posterior. The distance is the posterior Mahalanobis distance between the full-data and leave-one-out fitted-value vectors:

$D_i = \frac{\Delta_i' V_\mu^+ \Delta_i}{P}$

where $\Delta_i = \hat{\mu} - \hat{\mu}_{(-i)}$, $V_\mu^+$ is the generalized inverse of the full-posterior fitted-value covariance, and $P$ is the rank of the fixed-effect model matrix.

Value

A numeric vector of Cook's distance values, one for each observation.

See Also

influence.brma, dffits.brma, hatvalues.brma

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- bPET(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")
  fit <- add_loo(fit)

  cooks.distance(fit)
}

## End(Not run)

Description

Computes COVRATIO for a fitted brma object. COVRATIO measures the change in the determinant of the covariance matrix of the estimates when observation $i$ is removed.

Usage

## S3 method for class 'brma'
covratio(model, type = "mods", ...)

Arguments

Details

COVRATIO is computed using importance sampling weights to approximate the leave-one-out covariance matrices without refitting the model. Estimate-unit LOO must first be computed with model <- add_loo(model, unit = "estimate"), unless internal weights are supplied.

$COVRATIO_i = \frac{\det(Cov(\beta)_{-i})}{\det(Cov(\beta))}$

Values > 1 indicate that the observation improves precision (decreases variance), while values < 1 indicate that the observation decreases precision (increases variance). If any included parameter has zero posterior variance, or if a full or LOO covariance determinant is zero or non-finite, COVRATIO is undefined. In that case, values are reported as NaN with a printed note when available.

Value

A named numeric vector of COVRATIO values, one for each observation.

See Also

influence.brma, dffits.brma, cooks.distance.brma

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- bPET(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")
  fit <- add_loo(fit)

  covratio(fit)
}

## End(Not run)

Description

Shared data-input arguments used by the RoBMA fitting functions.

Normal models use approximate effect-size estimates supplied through yi with either vi or sei. GLMM models use the raw two-arm count arguments for binomial (measure = "OR") or Poisson (measure = "IRR") outcomes.

Arguments

Description

Computes DFBETAS (Difference in BETAS, standardized) for a fitted brma object. DFBETAS measures the influence of each observation on the estimated model coefficients. Positive values indicate that deleting the observation yields a smaller estimate, negative values indicate that deleting the observation yields a larger estimate.

Usage

## S3 method for class 'brma'
dfbetas(
  model,
  type = "mods",
  standardized_coefficients = FALSE,
  transform_factors = TRUE,
  return_loo_estimates = FALSE,
  ...
)

Arguments

Details

This function computes DFBETAS values using the Leave-One-Out (LOO) approximation based on Pareto Smoothed Importance Sampling (PSIS) weights. Ideally, DFBETAS is defined as:

$DFBETAS_{ij} = \frac{\hat{\beta}_j - \hat{\beta}_{j(-i)}}{SE(\hat{\beta}_{j(-i)})}$

where $\hat{\beta}_j$ is the estimate of the $j$-th coefficient using the full data, $\hat{\beta}_{j(-i)}$ is the estimate when observation $i$ is omitted, and $SE(\hat{\beta}_{j(-i)})$ is the standard error of the coefficient when observation $i$ is omitted.

In the Bayesian context using LOO approximation:

• $\hat{\beta}_{j(-i)}$ is estimated as the importance sampling weighted mean of the posterior samples, using PSIS weights $w_{is}$.

• $SE(\hat{\beta}_{j(-i)})$ is estimated as the importance sampling weighted standard deviation of the posterior samples.

This approximation allows computing influence statistics without refitting the model $K$ times, making it computationally efficient. For type = "bias", fixed identification parameters (e.g., the reference $\omega = 1$ interval) can have zero LOO posterior standard deviation. These parameters are retained in the output, but their DFBETAS values are reported as NaN because the standardized diagnostic is undefined.

Note: This function requires that LOO-CV has been computed for the model using add_loo.

Value

If return_loo_estimates = FALSE, a data frame with $K$ rows (observations) and $P$ columns (parameters), containing DFBETAS values. If return_loo_estimates = TRUE, returns the corresponding leave-one-out coefficient estimates. Row names correspond to study labels (if available) or indices.

See Also

add_loo, loo_weights.brma

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- bPET(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")
  fit <- add_loo(fit)

  inf <- dfbetas(fit)
  plot(inf[, 1], type = "h")
}

## End(Not run)

Description

Computes DFFITS (Difference in FITS, standardized) for a fitted brma object. DFFITS measures how much the fitted value for observation $i$ changes if observation $i$ is removed, standardized by the estimated standard error of the fit.

Usage

## S3 method for class 'brma'
dffits(model, ...)

Arguments

Details

DFFITS values are computed as a PSIS leave-one-out deletion diagnostic. For each observation $i$, the leave-one-out posterior mean fitted value at that observation is estimated with normalized PSIS weights and compared to the full-posterior fitted value:

$DFFITS_i = \frac{\hat{\mu}_i - \hat{\mu}_{i(-i)}}{SD_{(-i)}(\mu_i)}$

This targets deletion influence on fitted values directly. It does not use LOO-PIT residuals, which are predictive outlier diagnostics rather than fitted-value deletion diagnostics.

Estimate-unit LOO must first be computed with model <- add_loo(model, unit = "estimate"). If the leave-one-out posterior SD of a fitted value is near zero, the corresponding DFFITS value is returned as NA.

Value

A named numeric vector of DFFITS values, one for each observation.

See Also

influence.brma, cooks.distance.brma, hatvalues.brma

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- bPET(yi = yi, vi = vi, data = dat.lehmann2018, measure = "SMD")
  fit <- add_loo(fit)

  dffits(fit)
}

## End(Not run)

Description

Estimates the unit information standard deviation (UISD) from sample sizes and standard errors. The UISD is used to scale weakly informative prior distributions for meta-analytic parameters.

Usage

estimate_unit_information_sd(sei, ni)

Arguments

Details

The unit information standard deviation is computed following Equation 6 in Röver et al. (2021):

$\text{UISD} = \sqrt{\frac{\sum n_i}{\sum \text{se}_i^{-2}}}$

where $n_i$ is the sample size and $\text{se}_i$ is the standard error for each study.

This function is useful when you want to compute the UISD once and pass it to multiple analyses via the prior_unit_information_sd argument in brma() or related functions. This ensures consistent prior scaling across analyses performed on different subsets of the same data.

Value

Returns a single numeric value representing the estimated unit information standard deviation.

References

Röver C, Bender R, Dias S, Schmid CH, Schmidli H, Sturtz S, Weber S, Friede T (2021). “On weakly informative prior distributions for the heterogeneity parameter in Bayesian random-effects meta-analysis.” Research Synthesis Methods, 12(4), 448–474. doi:10.1002/jrsm.1475.

Examples

# Example with simulated data
sei <- c(0.2, 0.3, 0.25, 0.15)
ni  <- c(50, 30, 40, 80)
estimate_unit_information_sd(sei = sei, ni = ni)

Description

Extract in-sample fitted values from a fitted brma object.

Usage

## S3 method for class 'brma'
fitted(
  object,
  unit = "estimate",
  conditioning_depth = "marginal",
  component = "location",
  bias_adjusted = FALSE,
  output_measure = NULL,
  transform = NULL,
  conditional = FALSE,
  ...
)

Arguments

Details

This method is a compact adapter around predict.brma. It summarizes posterior prediction draws by their column means and returns a base numeric vector, matching the usual fitted contract. Use predict() directly when posterior draws or intervals are needed.

The default conditioning_depth = "marginal" corresponds to predict(object, type = "terms") and matches the usual fitted-value convention for meta-regression. For normal models, conditioning_depth = "estimate" corresponds to BLUP means for the observed estimates.

For component = "all", conditioning_depth, output_measure, and transform apply only to the location component. The scale component always returns fitted $\tau_i$ values.

Value

A numeric vector of fitted values, or a named list with location and scale components when component = "all".

See Also

predict.brma(), residuals.brma(), blup.brma()

Description

The brma family of functions uses the following arguments to specify the MCMC sampling and fitting settings.

Arguments

See Also

brma, set_autofit_control, set_convergence_checks

Description

funnel.brma creates a funnel plot for a fitted brma object. For intercept-only models without scale regression, the default outcome mode displays observed effect sizes against the fitted sampling distribution. For models with location or scale moderators, the default residual mode displays residuals against a standard-error funnel.

Usage

## S3 method for class 'brma'
funnel(
  x,
  residual,
  type = "LOO-PIT",
  unit = "estimate",
  conditioning_depth = "marginal",
  sampling_heterogeneity = TRUE,
  sampling_bias = TRUE,
  max_samples = 10000,
  plot_type = "base",
  ...
)

Arguments

Details

The funnel plot has two modes. If residual is not specified, the mode is chosen automatically from the fitted model: intercept-only models without scale regression use outcome mode, whereas models with location or scale moderators use residual mode.

Outcome mode (intercept-only models without scale regression): Displays observed effect sizes on the x-axis and standard errors on the y-axis. The reference line follows the center of the fitted sampling distribution. When sampling_bias = FALSE, this center is the pooled effect; when PET/PEESE bias adjustment is incorporated, the center line can vary with the standard error. The funnel region represents the central 95\ region of the sampling distribution, optionally incorporating heterogeneity and publication bias.

Residual mode (models with moderators or scale regression): Displays residuals on the x-axis and standard errors on the y-axis. The funnel region represents the central 95\ $N(0, \mathrm{SE}^2)$. With type = "LOO-PIT", the plotted residuals and standard errors are the raw-scale LOO predictive companions returned by rstudent.brma; the PIT-normalized z values are used by qqnorm.brma and influence diagnostics. With type = "rstandard", the plotted values are internally standardized residual companions from rstandard.brma. With type = "outcome", these are raw outcome residuals. Under a correctly specified model, most points should fall within this region.

The type argument controls how residuals are computed in residual mode. See residuals.brma for details on each type. The sampling_heterogeneity and sampling_bias arguments are ignored in residual mode.

For GLMM models, observed effect sizes are computed from the raw frequency data using formulas equivalent to metafor::escalc. Residual-mode GLMM funnels use approximate effect-size-scale residual/PIT companions, not exact PIT diagnostics for the raw count likelihood.

Value

If as_data = TRUE, funnel.brma returns a list with the data used for plotting, including the plotted points, funnel polygons, plotting limits, labels, and reference line. Otherwise, it returns NULL invisibly if plot_type = "base" or a ggplot object if plot_type = "ggplot".

See Also

residuals.brma(), rstandard.brma(), rstudent.brma(), predict.brma()

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE) &&
    requireNamespace("metafor", quietly = TRUE)) {
  data(dat.bcg, package = "metadat")
  dat <- metafor::escalc(
    measure = "RR",
    ai      = tpos,
    bi      = tneg,
    ci      = cpos,
    di      = cneg,
    data    = dat.bcg
  )

  fit <- brma(yi = yi, vi = vi, data = dat, measure = "RR")
  funnel(fit)
  funnel(fit, pch = 19, col = "blue", bg = "lightblue")

  fit_reg <- brma(
    yi      = yi,
    vi      = vi,
    mods    = ~ ablat + year,
    data    = dat,
    measure = "RR"
  )
  fit_reg <- add_loo(fit_reg)
  funnel(fit_reg)
  funnel(fit_reg, type = "outcome")

  funnel_data <- funnel(fit, as_data = TRUE)
  funnel(fit, plot_type = "ggplot")
}

## End(Not run)

Description

Computes hat values (leverages) from a fitted brma object. Returns posterior mean leverages, one for each observation.

Usage

## S3 method for class 'brma'
hatvalues(model, ...)

Arguments

Details

Hat values (leverages) measure the influence of each observation on the fitted values. In a Bayesian meta-analysis, the random effects variance $\tau^2$ is uncertain, so the hat matrix depends on the posterior samples of $\tau^2$.

This function computes the diagonal elements of the hat matrix:

$h_i = (X (X^T W X)^{-1} X^T W)_{ii}$

where $W$ is the weight matrix inverse to the marginal variance matrix.

The hat matrix is computed for each posterior draw and then averaged over draws, matching the vector output shape used by metafor.

This method is available only for normal outcome models without weightfunction selection.

Value

A numeric vector of length K, one posterior mean leverage per observation.

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE) &&
    requireNamespace("metafor", quietly = TRUE)) {
  data(dat.bcg, package = "metadat")
  dat <- metafor::escalc(
    measure = "RR",
    ai      = tpos,
    bi      = tneg,
    ci      = cpos,
    di      = cneg,
    data    = dat.bcg
  )

  fit <- brma(yi = yi, vi = vi, data = dat, measure = "RR")
  hatvalues(fit)
}

## End(Not run)

Description

The data set contains effect sizes (percent increase in earnings), standard errors, study identifiers, sample sizes, and the type of customer contact (no, some, or direct). The meta-analysis examined the relationship between perceived beauty and professional success (Havránková et al. 2025).

Usage

Havrankova2025

Format

A data.frame with 5 columns and 1159 observations:

y

Effect size: percent increase in earnings.

se

Standard error of y.

facing_customer

Type of customer contact.

study_id

Study identifier.

N

Sample size.

References

Havránková Z, Havránek T, Bortnikova K, Bartoš F (2025). “Meta-analysis of field studies on beauty and professional success.” Preprint available at https://meta-analysis.cz/beauty/beauty.pdf.

Description

Plots a histogram of observed z-values from the meta-analysis.

Usage

## S3 method for class 'zplot_brma'
hist(
  x,
  plot_type = "base",
  from = -6,
  to = 6,
  by = 0.5,
  length.out = NULL,
  add = FALSE,
  plot_thresholds = TRUE,
  dots_thresholds = NULL,
  dots_hist = NULL,
  dots_all = NULL,
  ...
)

Arguments

Details

Z-statistics are computed as effect size divided by standard error (yi / sei). Histogram bins are adjusted to align with significance thresholds when selection model priors are present.

Value

NULL invisibly for base graphics, or a ggplot2 object.

See Also

plot.zplot_brma(), lines.zplot_brma()

Description

The data set contains Cohen's d effect sizes, variances, and study characteristics including outcome type, feedback level, social comparison type, number of sessions, sample type, sample size, and country. The meta-analysis examined social comparison as a behavior change technique across the behavioral sciences (Hoppen et al. 2025).

Usage

Hoppen2025

Format

A data.frame with 9 columns and 37 observations:

d

Cohen's d effect size.

v

Sampling variance of d.

outcome

Outcome type.

feedback_level

Feedback level.

social_comparison_type

Social-comparison type.

sessions

Number of sessions.

sample_type

Sample type.

sample_size

Sample size.

country

Country.

References

Hoppen TH, Cuno RM, Nelson J, Lemmel F, Schlechter P, Morina N (2025). “Meta-analysis of randomized controlled trials examining social comparison as a behaviour change technique across the behavioural sciences.” Nature Human Behaviour, 9(8), 1595–1612. doi:10.1038/s41562-025-02209-2.

Description

Computes DFFITS, Cook's distance, COVRATIO, and other influence diagnostics for a fitted brma object.

Usage

## S3 method for class 'brma'
influence(model, ...)

Arguments

Value

An object of class "infl.brma", which corresponds to the structure of metafor::influence objects. It is a list containing:

Undefined determinant- or variance-standardized diagnostics are reported as NaN and printed with an explanatory note.

See Also

dffits.brma, cooks.distance.brma, covratio.brma, dfbetas.brma, hatvalues.brma

Examples

## Not run: 
if (requireNamespace("metadat", quietly = TRUE)) {
  data(dat.lehmann2018, package = "metadat")
  fit <- brma(
    yi      = yi,
    vi      = vi,
    data    = dat.lehmann2018,
    measure = "SMD",
    seed    = 1,
    silent  = TRUE
  )
  fit <- add_loo(fit)

  inf <- influence(fit)
  print(inf)
}

## End(Not run)

Description

Creates a concise textual interpretation of fitted RoBMA brma objects.

Usage

interpret(object, ...)

## Default S3 method:
interpret(object, ...)

## S3 method for class 'brma'
interpret(
  object,
  output_measure = NULL,
  transform = NULL,
  conditional = FALSE,
  scope = "core",
  probs = c(0.025, 0.975),
  central = NULL,
  priors = FALSE,
  digits = 3,
  ...
)

## S3 method for class 'interpret.brma'
print(x, ...)

Arguments

Value

A character vector with class "interpret.brma". The normalized BayesTools interpretation records are stored in the "records" attribute; the brma method also stores "scope", "conditional", and optionally "priors".

Description

The data set contains Cohen's d effect sizes, standard errors, and study labels from a meta-analysis investigating the extent to which family-based treatments for children with mental health, physical health, and developmental disorders provide secondary benefits to the children's siblings and caregivers (Johnides et al. 2025).

Usage

Johnides2025

Format

A data.frame with 3 columns and 412 observations:

study

Study label.

d

Cohen's d effect size.

se

Standard error of d.

References

Johnides BD, Borduin CM, Sheerin KM, Kuppens S (2025). “Secondary benefits of family member participation in treatments for childhood disorders: A multilevel meta-analytic review.” Psychological Bulletin, 151(1), 1–32. doi:10.1037/bul0000462.

Description

The data set contains partial correlation coefficients, standard errors, study labels, sample sizes, type of the educational outcome, intensity of the employment, gender of the student population, study location, study design, whether the study controlled for endogeneity, and whether the study controlled for motivation. The original data set including additional variables and the publication are available from the project page. (Note that some standard errors and employment intensities are missing.)

Usage

Kroupova2021

Format

A data.frame with 11 columns and 881 observations:

r

Partial correlation coefficient.

se

Standard error of r.

study

Study label.

sample_size

Sample size.

education_outcome

Type of educational outcome.

employment_intensity

Employment intensity.

students_gender

Gender composition of the student sample.

location

Study location.

design

Study design.

endogenity_control

Whether endogeneity was controlled.

motivation_control

Whether motivation was controlled.

Source

http://meta-analysis.cz/students

References

Kroupova K, Havranek T, Irsova Z (2021). “Student employment and education: A meta-analysis.” CEPR Discussion Paper. https://www.ssrn.com/abstract=3928863.

Description

Adds model-implied density lines to an existing zplot.

Usage

## S3 method for class 'zplot_brma'
lines(
  x,
  plot_type = "base",
  probs = c(0.025, 0.975),
  max_samples = 10000,
  plot_ci = TRUE,
  extrapolate = FALSE,
  from = -6,
  to = 6,
  by = 0.05,
  length.out = NULL,
  col = "black",
  as_data = FALSE,
  ...
)

Arguments

Details

When extrapolate = FALSE, the density includes all bias adjustments (PET/PEESE regression, selection weights) representing the fitted model. When extrapolate = TRUE, bias adjustments are removed to show the hypothetical distribution without publication bias. Under selection models, the curve is scaled by inverse selection probability and need not integrate to one.

Value

NULL invisibly for base graphics, ggplot2 layers for ggplot, or a data frame with columns x, y, y_lCI, y_uCI if as_data = TRUE.

See Also

plot.zplot_brma(), hist.zplot_brma()

Description

Extract the pointwise log-likelihood matrix from a brma model object. This is an S x K or S x G matrix where S is the number of posterior samples, K is the number of estimates, and G is the number of clusters. This method implements the S3 \'logLik\' generic for brma objects and returns the matrix of pointwise log-likelihoods (one column per observation, one row per sample).

Usage

## S3 method for class 'brma'
logLik(object, unit = "estimate", ...)