bma_inference: Bayesian model averaging

In bayesics: Bayesian Analyses for One- and Two-Sample Inference and Regression Methods

Bayesian model averaging

Description

Estimates and CIs from BMA

Usage

bma_inference(
  formula,
  data,
  zellner_g = nrow(data),
  CI_level = 0.95,
  ROPE,
  mcmc_draws = 10000,
  n_models = 500,
  mc_error = 0.001,
  seed = 1,
  ...
)

Arguments

formula	A formula specifying the model.
data	Data used in linear regression model
zellner_g	numeric. Positive number giving the value of "g" in Zellner's g prior.
CI_level	Level for credible interval
ROPE	vector of positive values giving ROPE boundaries for each regression coefficient. Optionally, you can not include a ROPE boundary for the intercept. If missing, defaults go to those suggested by Kruchke (2018).
mcmc_draws	Integer. Number of draws passed into bms
n_models	Integer. The number of best models for which information is stored. See bms for more details.
mc_error	The number of posterior draws will ensure that with 99% probability the bounds of the credible intervals will be within \pm mc_error\times 4s_y, that is, within 100mc_error% of the trimmed range of y.
seed	Integer. Always set your seed!!!
...	Other arguments for bms.

Details

bma_inference leverages the bms function from its eponymous R package, and then uses lm_b to obtain inference on the regression coefficients for Bayesian model averaging.

Value

A list with the following elements:

• summary Tibble with point and interval estimates

• lm_b_fits A list of lm_b fits using zellner's g prior for all the top models from bms

• hyperparameters A named list with the user-specified zellner's g value.

• posterior_draws matrix of posterior draws of the regression parameters, marginalizing out the model

Examples

# Create data
set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5],
             x4 = rnorm(N),
             x5 = rnorm(N),
             x6 = rnorm(N),
             x7 = rnorm(N),
             x8 = rnorm(N),
             x9 = rnorm(N),
             x10 = rnorm(N))
test_data$outcome = 
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )

# Fit linear model using Bayesian model averaging
fit <-
  bma_inference(outcome ~ .,
                test_data,
                user.int = FALSE)
summary(fit)
coef(fit)
credint(fit)
plot(fit)

bayesics documentation built on March 11, 2026, 5:07 p.m.