Bayesian Modeling in R

Overview

Comprehensive Bayesian statistical modeling using Stan-based packages (brms, rstanarm), covering prior specification, posterior analysis, model comparison, and Bayesian workflow best practices.

brms: Bayesian Regression Models

Basic Models

library(brms)

# Linear regression
fit <- brm(
  formula = y ~ x1 + x2,
  data = df,
  family = gaussian(),
  seed = 123
)

# Logistic regression
fit_logit <- brm(
  y ~ x1 + x2,
  data = df,
  family = bernoulli(link = "logit")
)

# Poisson regression
fit_pois <- brm(
  count ~ x1 + x2 + offset(log(exposure)),
  data = df,
  family = poisson()
)

Prior Specification

# View default priors
get_prior(y ~ x1 + x2, data = df, family = gaussian())

# Set custom priors
custom_priors <- c(
  prior(normal(0, 10), class = "Intercept"),
  prior(normal(0, 2), class = "b"),  # All regression coefficients
  prior(normal(0, 1), class = "b", coef = "x1"),  # Specific coefficient
  prior(exponential(1), class = "sigma")  # Error SD
)

fit <- brm(
  y ~ x1 + x2,
  data = df,
  family = gaussian(),
  prior = custom_priors,
  seed = 123
)

Prior Predictive Checks

# Sample from prior only
fit_prior <- brm(
  y ~ x1 + x2,
  data = df,
  family = gaussian(),
  prior = custom_priors,
  sample_prior = "only",  # Prior predictive
  seed = 123
)

# Visualize prior predictions
pp_check(fit_prior, type = "dens_overlay", ndraws = 100)

Mixed Effects Models

# Random intercepts
fit_mixed <- brm(
  y ~ x1 + x2 + (1 | group),
  data = df,
  family = gaussian()
)

# Random slopes
fit_mixed <- brm(
  y ~ x1 + x2 + (1 + x1 | group),
  data = df,
  family = gaussian()
)

# Crossed random effects
fit_mixed <- brm(
  y ~ x1 + (1 | subject) + (1 | item),
  data = df,
  family = gaussian()
)

Control Parameters

fit <- brm(
  y ~ x1 + x2,
  data = df,
  family = gaussian(),
  chains = 4,
  iter = 4000,
  warmup = 2000,
  cores = 4,
  seed = 123,
  control = list(
    adapt_delta = 0.95,  # Higher for problematic posteriors
    max_treedepth = 15
  )
)

rstanarm: Applied Regression

library(rstanarm)

# Linear regression
fit <- stan_glm(
  y ~ x1 + x2,
  data = df,
  family = gaussian(),
  prior = normal(0, 2.5),
  prior_intercept = normal(0, 10),
  seed = 123
)

# Mixed effects
fit_mixed <- stan_lmer(
  y ~ x1 + x2 + (1 | group),
  data = df,
  seed = 123
)

# Generalized linear mixed
fit_glmer <- stan_glmer(
  y ~ x1 + x2 + (1 | group),
  data = df,
  family = binomial(),
  seed = 123
)

Posterior Analysis

Summary and Inference

# Model summary
summary(fit)

# Posterior draws
posterior <- as_draws_df(fit)

# Posterior summary
posterior_summary(fit)

# Fixed effects
fixef(fit)

# Random effects
ranef(fit)

# Credible intervals
posterior_interval(fit, prob = 0.95)

Hypothesis Testing

# Probability statements
hypothesis(fit, "x1 > 0")
hypothesis(fit, "x1 > x2")
hypothesis(fit, "x1 + x2 > 0")

# Multiple hypotheses
hypothesis(fit, c("x1 > 0", "x2 > 0", "x1 > x2"))

Posterior Predictive Checks

library(bayesplot)

# Density overlay
pp_check(fit, type = "dens_overlay", ndraws = 50)

# Histogram
pp_check(fit, type = "hist", ndraws = 8)

# Error scatter
pp_check(fit, type = "error_scatter_avg")

# Intervals
pp_check(fit, type = "intervals")

# Stat comparison
pp_check(fit, type = "stat", stat = "mean")
pp_check(fit, type = "stat_2d", stat = c("mean", "sd"))

MCMC Diagnostics

library(bayesplot)

# Trace plots
mcmc_trace(fit)

# Rhat
rhat(fit)
mcmc_rhat(rhat(fit))

# Effective sample size
neff_ratio(fit)
mcmc_neff(neff_ratio(fit))

# Pairs plot (divergences)
mcmc_pairs(fit, pars = c("b_x1", "b_x2", "sigma"))

# Energy
mcmc_nuts_energy(nuts_params(fit))

Model Comparison

LOO Cross-Validation

# LOO-CV
loo_fit1 <- loo(fit1)
loo_fit2 <- loo(fit2)

# Compare models
loo_compare(loo_fit1, loo_fit2)

# Pareto k diagnostics
plot(loo_fit1)

WAIC

# WAIC
waic_fit1 <- waic(fit1)
waic_fit2 <- waic(fit2)

# Compare
loo_compare(waic_fit1, waic_fit2)

Bayes Factors

library(bridgesampling)

# Compute marginal likelihood
bridge_fit1 <- bridge_sampler(fit1)
bridge_fit2 <- bridge_sampler(fit2)

# Bayes factor
bayes_factor(bridge_fit1, bridge_fit2)

Model Stacking

library(loo)

# Model weights based on LOO
model_weights <- loo_model_weights(
  list(fit1, fit2, fit3),
  method = "stacking"
)

Predictions

Posterior Predictions

# Expected values (fitted)
fitted(fit, newdata = new_data)

# Predictions with uncertainty
predict(fit, newdata = new_data)

# Full posterior predictive draws
pp_draws <- posterior_predict(fit, newdata = new_data)

Marginal Effects

# Conditional effects
conditional_effects(fit)

# Specific effects
conditional_effects(fit, effects = "x1")

# Interaction effects
conditional_effects(fit, effects = "x1:x2")

# Plot with data
plot(conditional_effects(fit, effects = "x1"), points = TRUE)

Marginal Means

library(emmeans)

# Estimated marginal means
emmeans(fit, ~ treatment)

# Contrasts
emmeans(fit, pairwise ~ treatment)

Advanced Topics

Non-Linear Models

# Non-linear formula
fit_nl <- brm(
  bf(y ~ a * exp(-b * x), a ~ 1, b ~ 1, nl = TRUE),
  data = df,
  prior = c(
    prior(normal(10, 5), nlpar = "a"),
    prior(normal(0.5, 0.2), nlpar = "b")
  ),
  family = gaussian()
)

Distributional Models

# Model mean and variance
fit_dist <- brm(
  bf(y ~ x1 + x2, sigma ~ x1),  # Heteroscedasticity
  data = df,
  family = gaussian()
)

Ordinal Regression

fit_ordinal <- brm(
  rating ~ x1 + x2,
  data = df,
  family = cumulative("logit")
)

Zero-Inflated Models

fit_zi <- brm(
  count ~ x1 + x2,
  data = df,
  family = zero_inflated_poisson()
)

# With predictors for zero-inflation
fit_zi <- brm(
  bf(count ~ x1 + x2, zi ~ x3),
  data = df,
  family = zero_inflated_poisson()
)

Survival Models

# Cox model (brms)
fit_surv <- brm(
  time | cens(censored) ~ x1 + x2,
  data = df,
  family = cox()
)

# Parametric survival
fit_weibull <- brm(
  time | cens(censored) ~ x1 + x2,
  data = df,
  family = weibull()
)

Bayesian Workflow

Complete Workflow Example

library(brms)
library(bayesplot)

# 1. Prior predictive check
fit_prior <- brm(
  y ~ x1 + x2,
  data = df,
  family = gaussian(),
  prior = c(
    prior(normal(0, 10), class = "Intercept"),
    prior(normal(0, 2), class = "b"),
    prior(exponential(1), class = "sigma")
  ),
  sample_prior = "only",
  seed = 123
)
pp_check(fit_prior, ndraws = 50)

# 2. Fit model
fit <- update(fit_prior, sample_prior = "no")

# 3. Check convergence
summary(fit)
mcmc_trace(fit)

# 4. Posterior predictive check
pp_check(fit, ndraws = 50)

# 5. Model comparison
fit_alt <- brm(y ~ x1, data = df, family = gaussian())
loo_compare(loo(fit), loo(fit_alt))

# 6. Inference
fixef(fit)
hypothesis(fit, "x1 > 0")

# 7. Predictions
conditional_effects(fit)

Key Packages Summary