Bayesian Linear Regression Using brms

Introduction

This example fits a Bayesian linear regression model using the brms package.

The model is

\[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i, \]

where

\[ \epsilon_i \sim N(0,\sigma^2). \]

Load Package and Simulate Data

library(brms)

set.seed(123)

n <- 31

dat <- data.frame(
  x = seq(-15, 15, length.out = n)
)

dat$y <- 5 * dat$x + rnorm(
  n,
  mean = 0,
  sd = 10
)

Plot the Simulated Data

plot(
  dat$x,
  dat$y,
  pch = 19,
  xlab = "x",
  ylab = "y",
  main = "Simulated Data"
)

Define Prior Distributions

The prior distributions are

\[ \beta_1 \sim N(5,2^2), \]

\[ \beta_0 \sim N(0,5^2), \]

and

\[ \sigma \sim \text{Exponential}(0.1). \]

priors <- c(
  brms::set_prior(
    "normal(5, 2)",
    class = "b"
  ),
  brms::set_prior(
    "normal(0, 5)",
    class = "Intercept"
  ),
  brms::set_prior(
    "exponential(0.1)",
    class = "sigma"
  )
)

Fit the Bayesian Linear Regression Model

fit <- brms::brm(
  y ~ x,
  data = dat,
  family = gaussian(),
  prior = priors,
  chains = 4,
  iter = 4000,
  warmup = 2000,
  cores = 4,
  seed = 123,
  refresh = 0
)

## Compiling Stan program...

## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.4.1 from
## https://cran.r-project.org/bin/windows/Rtools/.

## Trying to compile a simple C file

## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.4.1 from
## https://cran.r-project.org/bin/windows/Rtools/.

## Start sampling

Model Summary

summary(fit)

##  Family: gaussian 
##   Links: mu = identity 
## Formula: y ~ x 
##    Data: dat (Number of observations: 31) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.31      1.68    -3.57     3.00 1.00     7413     5971
## x             4.81      0.20     4.41     5.20 1.00     7000     5076
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     9.94      1.35     7.71    12.97 1.00     7176     5653
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

MCMC Diagnostic Plots

# Close all active graphics devices
while (dev.cur() > 1) {
  dev.off()
}

# Draw brms diagnostic plots
plot(fit)

Interpretation

The posterior estimate of b_x represents the slope of the relationship between x and y.

The posterior estimate of Intercept represents the expected value of y when (x=0).

The posterior estimate of sigma represents the residual standard deviation.