bayesics (version 2.0.2)

bma_inference: 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,
  ...
)

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

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.

Examples

# \donttest{
# 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)
# }