lm_b: Bayesian Linear Models

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

Bayesian Linear Models

Description

lm_b is used to fit linear models. It can be used to carry out regression, single stratum analysis of variance and analysis of covariance (although aov_b may provide a more convenient interface for ANOVA.)

Usage

lm_b(
  formula,
  data,
  weights,
  prior = c("zellner", "conjugate", "improper")[1],
  zellner_g,
  prior_beta_mean,
  prior_beta_precision,
  prior_var_shape,
  prior_var_rate,
  ROPE,
  CI_level = 0.95
)

Arguments

formula	A formula specifying the model.
data	A data frame in which the variables specified in the formula will be found. If missing, the variables are searched for in the standard way. However, it is strongly recommended that you use this argument so that other generics for bayesics objects work correctly.
weights	an optional vector of weights to be used in the fitting process. Should be NULL or a numeric vector. If non-NULL, it is assumed that the variance of y_i can be written as Var(y_i) = \sigma^2/w_i. While the estimands remain the same, the estimation is done by performing unweighted lm_b on W^{\frac{1}{2}}y and W^{\frac{1}{2}}X, where W is the diagonal matrix of weights. Note that this then affects the zellner prior.
prior	character. One of "zellner", "conjugate", or "improper", giving the type of prior used on the regression coefficients.
zellner_g	numeric. Positive number giving the value of "g" in Zellner's g prior. Ignored unless prior = "zellner".
prior_beta_mean	numeric vector of same length as regression coefficients (denoted p). Unless otherwise specified, automatically set to rep(0,p). Ignored unless prior = "conjugate".
prior_beta_precision	pxp matrix giving a priori precision matrix to be scaled by the residual precision.
prior_var_shape	numeric. Twice the shape parameter for the inverse gamma prior on the residual variance(s). I.e., \sigma^2\sim IG(prior_var_shape/2,prior_var_rate/2).
prior_var_rate	numeric. Twice the rate parameter for the inverse gamma prior on the residual variance(s). I.e., \sigma^2\sim IG(prior_var_shape/2,prior_var_rate/2).
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).
CI_level	numeric. Credible interval level.

Details

MODEL:

The likelihood is given by

y_i \overset{ind}{\sim} N(x_i'\beta,\sigma^2).

The prior is given by

\beta|\sigma^2 \sim N\left( \mu, \sigma^2 V^{-1} \right) \\ \sigma^2 \sim \Gamma^{-1}(a/2,b/2).

• For Zellner's g prior, \frac{1}{g}X'X.

• The default for the conjugate prior is based on arguments from standardized regression. The default V is set according to the following: "a priori, we are 95% certain that a standard deviation increase in X will not lead to more than a 5 standard deviation in the mean of y." This leads to

  V^{1/2} = diag(v^{1/2},s_{x_1},\ldots,s_{x_p}) / 2.5,

  v^{-1} := \sum \frac{\bar x_j^2}{s_j^2}

  where s_y is the standard deviation of y, and s_{x_j} is the standard deviation of the j^{th} covariate.

• Unless prior_var_shape AND prior_var_rate are provided, the inverse gamma prior on the residual variance will place 50% prior probability that the correlation between the response and the fitted values is between 0.1 and 0.9.

• If prior = "improper", then the prior is

  \pi(\beta,\sigma^2) \propto \frac{1}{\sigma^2}.

ROPE:

If missing, the ROPE bounds will be given under the principle of "half of a small effect size." Using Cohen's D of 0.2 as a small effect size, the ROPE is built under the principle that moving the full range of X (i.e., \pm 2 s_x) will not move the mean of y by more than the overall mean of y minus 0.1s_y to the overall mean of y plus 0.1s_y. The result is a ROPE equal to |\beta_j| < 0.05s_y/s_j. If the covariate is binary, then this is simply |\beta_j| < 0.2s_y.

Value

lm_b() returns an object of class "lm_b", which behaves as a list with the following elements:

• summary - tibble giving results for regression coefficients

• posterior_parameters - list giving the posterior parameters

• hyperparameters - list giving the user input (or default) hyperparameters used

• fitted - posterior mean of the individuals' means

• residuals - posterior mean of the residuals

• formula, data - input by user

Examples

# Generate some data
set.seed(2025)
N = 500
test_data = 
  data.frame(x1 = rnorm(N),
             x2 = rnorm(N),
             x3 = letters[1:5])
test_data$outcome = 
  rnorm(N,-1 + test_data$x1 + 2 * (test_data$x3 %in% c("d","e")) )

# Find good hyperparameters for the residual variance
s2_hyperparameters = 
  find_invgamma_parms(lower_R2 = 0.05,
                      upper_R2 = 0.95,
                      probability = 0.8,
                      response_variance = var(test_data$outcome))
## Check (should equal ~ 0.8)
extraDistr::pinvgamma((1.0 - 0.05) * var(test_data$outcome),
                      0.5 * s2_hyperparameters[1],
                      0.5 * s2_hyperparameters[2]) -
  extraDistr::pinvgamma((1.0 - 0.95) * var(test_data$outcome),
                        0.5 * s2_hyperparameters[1],
                        0.5 * s2_hyperparameters[2])

# Fit the linear model using a conjugate prior
fit1 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data,
       prior = "conj",
       prior_var_shape = s2_hyperparameters["shape"],
       prior_var_rate = s2_hyperparameters["rate"])
fit1
summary(fit1,
        CI_level = 0.99)
plot(fit1)
coef(fit1)
credint(fit1,
        CI_level = 0.9)
bayes_factors(fit1)
bayes_factors(fit1,
              by = "var")
AIC(fit1)
BIC(fit1)
DIC(fit1)
WAIC(fit1)
vcov(fit1)
preds = predict(fit1)

# Try other priors
fit2 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data) # Default is prior = "zellner"
fit3 <-
  lm_b(outcome ~ x1 + x2 + x3,
       data = test_data,
       prior = "improper")

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