Models with a Single Predictor

Now we will build a model in Stan to formally estimate this relationship. Rather than creating an r code block, we want to create a stan code block. The only caveat is that we also need to add a name for our Stan model provided in the output.var argument.

Here is the model that we want to estimate in Stan:

\[ \begin{aligned} \text{mpg} &\sim \text{Normal}(\mu, \sigma) \\ \mu &= \alpha + \beta \text{weight} \end{aligned} \]

To start, we need to load cmdstanr which will allow us to interface with Stan through R.

library(cmdstanr)

## This is cmdstanr version 0.4.0

## - Online documentation and vignettes at mc-stan.org/cmdstanr

## - CmdStan path set to: /Users/nickjenkins/.cmdstanr/cmdstan-2.28.1

## - Use set_cmdstan_path() to change the path

## 
## A newer version of CmdStan is available. See ?install_cmdstan() to install it.
## To disable this check set option or environment variable CMDSTANR_NO_VER_CHECK=TRUE.

register_knitr_engine(override = FALSE) # This registers cmdstanr with knitr so that we can use
# it with R Markdown.

Now, program the model:

data {
  int n; //number of observations in the data
  vector[n] mpg; //vector of length n for the car's MPG
  vector[n] weight; //vector of length n for the car's weight
}

parameters {
  real alpha; //the intercept parameter
  real beta_w; //slope parameter for weight
  real sigma; //model variance parameter
}

model {
  //linear predictor mu
  vector[n] mu;
  
  //write the linear equation
  mu = alpha + beta_w * weight;
  
  //likelihood function
  mpg ~ normal(mu, sigma);
}

Once we finish writing the model, we need to run the code block to compile it into C++ code. This will also us to sample from the model and obtain the parameter estimates. Let’s do that now.

The next step is to prepare the data for Stan. Stan can’t use the same types of data that R can. For example, Stan requires lists, not data frames, and it cannot accept factors. We’ll use the tidybayes package to make it easier to prepare the data.

library(tidybayes)

model.data <- 
  cars.data %>%
  select(mpg, weight) %>%
  compose_data(.)

# sample from our model
linear.fit.1 <- linear.model$sample(data = model.data, refresh = 0)

## Running MCMC with 4 sequential chains...
## 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.7 seconds.

# summarize our model
print(linear.fit.1)

##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -52.26 -51.88 1.39 1.05 -54.91 -50.80 1.00     1145     1622
##    alpha   37.33  37.30 1.99 1.99  34.11  40.63 1.00     1079     1104
##    beta_w  -5.36  -5.35 0.59 0.59  -6.33  -4.41 1.00     1126     1149
##    sigma    3.20   3.14 0.46 0.42   2.57   4.00 1.00     1374     1236

Let’s run through the interpretation of this model:

• alpha For a car with a weight of zero, the expected MPG is 37.21. Obviously, a weight of zero is impossible, so we’ll want to address this in our next model.

• beta_w Comparing two cars who differ by 1000 pounds, the model predicts a difference of 5.33 miles per gallon.

• sigma The model predicts MPG within 3.18 points.

• lp__ Logarithm of the (unnormalized) posterior density. This log density can be used in various ways for model evaluation and comparison.

Ok, now that we’ve written our model, let’s make a few improvements. First, let’s center our weight variable so that we can get a more meaningful interpretation of the intercept parameter alpha. We can accomplish this by subtracting the mean from each observation. This will change the interpretation of the intercept to be the average MPG when weight is held constant at it’s average value.

cars.data <- 
  cars.data %>%
  mutate(weight_c = weight - mean(weight))

Now that we changed the name of the variable name, we also need to change our model code to incorporate this change. While we are adjusting the code, we’ll also restrict the scale parameter to be positive. This will help our model be a bit more efficient.

data {
  int n; //number of observations in the data
  vector[n] mpg; //vector of length n for the car's MPG
  vector[n] weight_c; //vector of length n for the car's weight
}

parameters {
  real alpha; //the intercept parameter
  real beta_w; //slope parameter for weight
  real<lower = 0> sigma; //variance parameter and restrict it to positive values
}

model {
  //linear predictor mu
  vector[n] mu;
  
  //write the linear equation
  mu = alpha + beta_w * weight_c;
  
  //likelihood function
  mpg ~ normal(mu, sigma);
}

Now prepare the data and re-estimate the model.

model.data <- 
  cars.data %>%
  select(mpg, weight_c) %>%
  compose_data(cars.data)

linear.fit.2 <- linear.model$sample(data = model.data, refresh = 0)

## Running MCMC with 4 sequential chains...
## 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.6 seconds.

print(linear.fit.2)

##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -51.06 -50.75 1.24 1.02 -53.53 -49.70 1.00     1681     2326
##    alpha   20.09  20.10 0.57 0.56  19.15  21.01 1.00     3359     2424
##    beta_w  -5.35  -5.36 0.57 0.56  -6.27  -4.40 1.00     3278     2625
##    sigma    3.18   3.13 0.43 0.40   2.56   3.95 1.00     3232     2957

Let’s interpret this model:

• alpha When a vehicle’s weight is held at its average value, the expected MPG is 20.09.

• beta_w This estimate has the same interpretation as before.

• sigma This estimate has the same interpretation as before.

• lp__ This estimate has a slightly lower value (in absolute value) than it did in the previous model indicating that this model performs slightly better.

Models with Multiple Predictors

To add more predictors, we just need to adjust out model code. Let’s add in the vehicle’s cylinders and horsepower. We’ll also center these variables.

data {
  int n; //number of observations in the data
  vector[n] mpg; //vector of length n for the car's MPG
  vector[n] weight_c; //vector of length n for the car's weight
  vector[n] cylinders_c; ////vector of length n for the car's cylinders
  vector[n] hp_c; //vector of length n for the car's horsepower
}

parameters {
  real alpha; //the intercept parameter
  real beta_w; //slope parameter for weight
  real beta_cyl; //slope parameter for cylinder
  real beta_hp; //slope parameter for horsepower
  real<lower = 0> sigma; //variance parameter and restrict it to positive values
}

model {
  //linear predictor mu
  vector[n] mu;
  
  //write the linear equation
  mu = alpha + beta_w * weight_c + beta_cyl * cylinders_c + beta_hp * hp_c;
  
  //likelihood function
  mpg ~ normal(mu, sigma);
}

Prepare the data and sample from the model.

model.data <- 
  cars.data %>%
  mutate(cylinders_c = cylinders - mean(cylinders),
         hp_c = hp - mean(hp)) %>%
  select(mpg, weight_c, cylinders_c, hp_c) %>%
  compose_data(.)

linear.fit.3 <- linear.model$sample(data = model.data, refresh = 0)

## Running MCMC with 4 sequential chains...
## 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.9 seconds.

print(linear.fit.3)

##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##  lp__     -45.28 -44.95 1.73 1.62 -48.59 -43.12 1.00     1428     2426
##  alpha     20.09  20.09 0.47 0.45  19.33  20.84 1.00     3536     2507
##  beta_w    -3.16  -3.16 0.80 0.81  -4.49  -1.85 1.00     3482     2566
##  beta_cyl  -0.96  -0.96 0.61 0.57  -1.98   0.04 1.00     2700     2218
##  beta_hp   -0.02  -0.02 0.01 0.01  -0.04   0.00 1.00     2892     2425
##  sigma      2.62   2.58 0.38 0.36   2.07   3.33 1.00     2546     2622

After adjusting for a car’s cylinders and horsepower two cars that differ by 1000 pounds, the model predicts a difference of 3.18 miles per gallon. Notice that lp__ is now even lower, suggesting that this latest model is a better fit.

Model Convergence

Now that we’ve built a decent model, we need to see how well it actual preforms. First, we’ll want to check that our chains have converged and are producing reliable point estimates. We can do this with a traceplot.

library(bayesplot)

## This is bayesplot version 1.8.1

## - Online documentation and vignettes at mc-stan.org/bayesplot

## - bayesplot theme set to bayesplot::theme_default()

##    * Does _not_ affect other ggplot2 plots

##    * See ?bayesplot_theme_set for details on theme setting

fit.draws <- linear.fit.3$draws() # extract the posterior draws
mcmc_trace(fit.draws)

The fuzzy caterpillar appearance indicates that the chains are mixing well and have converged to a common distribution. We can also assess the Rhat values for each parameter. As a rule of thumb, Rhat values less than 1.05 indicate good convergence. The bayesplot package makes these calculations easy.

rhats <- rhat(linear.fit.3)
mcmc_rhat(rhats)

Effective Sample Size

The effective sample size estimates the number of independent draws from the posterior distribution of a given estimate. This metric is important because Markov chains can have autocorrelation which will lead to biased parameter estimates. With the bayesplot package we can visualize the ratio of the effective sample size to the total number of samples - the larger the ratio the better. The rule of thumb here is to worry about ratios less than 0.1.

eff.ratio <- neff_ratio(linear.fit.3)
eff.ratio

##     alpha    beta_w  beta_cyl   beta_hp     sigma 
## 0.8839743 0.8695231 0.6686701 0.7003035 0.6309044

mcmc_neff(eff.ratio)

We can also check the autocorrelation in the chains with bayesplot. To use the mcmc_acf function, we’ll need to extract the posterior draws from the model.

mcmc_acf(fit.draws)

Here, we are looking to see how quickly the autocorrelation drops to zero.

Posterior Predictive Checks

One of the most powerful tools of Bayesian inference is to conduct posterior predictive checks. This check is designed to see how well our model can generate data that matches observed data. If we built a good model, it should be able to generate new observations that very closely resemble the observed data.

In order to perform posterior predictive checks, we will need to add in some code to our model. Specifically we need to calculate replications of our outcome variable. We can do this using the generated quantities section.

data {
  int n; //number of observations in the data
  vector[n] mpg; //vector of length n for the car's MPG
  vector[n] weight_c; //vector of length n for the car's weight
  vector[n] cylinders_c; ////vector of length n for the car's cylinders
  vector[n] hp_c; //vector of length n for the car's horsepower
}

parameters {
  real alpha; //the intercept parameter
  real beta_w; //slope parameter for weight
  real beta_cyl; //slope parameter for cylinder
  real beta_hp; //slope parameter for horsepower
  real<lower = 0> sigma; //variance parameter and restrict it to positive values
}

model {
  //linear predictor mu
  vector[n] mu;
  
  //write the linear equation
  mu = alpha + beta_w * weight_c + beta_cyl * cylinders_c + beta_hp * hp_c;
  
  //likelihood function
  mpg ~ normal(mu, sigma);
}

generated quantities {
  //replications for the posterior predictive distribution
  real y_rep[n] = normal_rng(alpha + beta_w * weight_c + beta_cyl * 
    cylinders_c + beta_hp * hp_c, sigma);
}

In the code block above, normal_rng is the Stan function to generate observations from a normal distribution. So, y_rep generates new data points from a normal distribution using the linear model we built mu and a variance sigma. Now let’s re-estimate the model:

linear.fit.3 <- linear.model$sample(data = model.data, refresh = 0)

## Running MCMC with 4 sequential chains...
## 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.9 seconds.

print(linear.fit.3)

##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##  lp__     -45.27 -44.93 1.77 1.62 -48.58 -43.13 1.01     1340     1990
##  alpha     20.08  20.07 0.49 0.48  19.25  20.87 1.00     3018     2417
##  beta_w    -3.15  -3.16 0.77 0.76  -4.41  -1.88 1.00     2949     2618
##  beta_cyl  -0.95  -0.96 0.58 0.57  -1.89   0.00 1.00     2458     1928
##  beta_hp   -0.02  -0.02 0.01 0.01  -0.04   0.00 1.00     2763     2393
##  sigma      2.64   2.60 0.38 0.36   2.10   3.33 1.00     2681     2335
##  y_rep[1]  22.85  22.85 2.73 2.65  18.31  27.40 1.00     4020     4054
##  y_rep[2]  21.97  22.02 2.74 2.72  17.50  26.46 1.00     3903     3851
##  y_rep[3]  25.88  25.86 2.76 2.68  21.44  30.42 1.00     3879     3712
##  y_rep[4]  20.89  20.92 2.74 2.65  16.31  25.28 1.00     3798     3738
## 
##  # showing 10 of 38 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

In our model output, we now have a replicated y value for every row of data. We can use these values to plot the replicated data against the observed data.

y <- cars.data$mpg

# convert the cmdstanr fit to an rstan object
library(rstan)

## Loading required package: StanHeaders

## rstan (Version 2.21.3, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

## 
## Attaching package: 'rstan'

## The following object is masked from 'package:tidyr':
## 
##     extract

stanfit <- read_stan_csv(linear.fit.3$output_files())

# extract the fitted values
y.rep <- extract(stanfit)[["y_rep"]]

ppc_dens_overlay(y = cars.data$mpg, yrep = y.rep[1:100, ])

The closer the replicated values (yrep) get to the observed values (y) the more accurate the model. Here it looks like we could probably do a bit better, though the lose fit is likely due to the small sample size (which adds more uncertainty).

Coeficient Plot

stan_plot(linear.fit.4, 
          pars = c("alpha", "beta_w", "beta_cyl", "beta_hp", "sigma"))

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

Fitted Regression Line

Here we can plot the fitted regression line

# Fitted Line
ggplot(data = cars.data, aes(x = weight_c, y = mpg)) +
  geom_point() +
  stat_function(fun = function(x) mean(posterior$alpha) + mean(posterior$beta_w) * x)

# Fitted Line with Uncertainty ------------------------------------------------
fit.plot <- 
  ggplot(data = cars.data, aes(x = weight_c, y = mpg)) +
  geom_point()

# select a random sample of 100 draws from the posterior distribution
sims <-
  posterior %>%
  mutate(n = row_number()) %>%
  sample_n(size = 100)

# add these draws to the plot
lines <- 
  purrr::map(1:100, function(i) stat_function(fun = function(x) sims[i, 1] + sims[i, 2] * x, 
             size = 0.08, color = "gray"))
fit.plot <- fit.plot + lines

# add the mean line to the plot
fit.plot <- 
  fit.plot +
  stat_function(fun = function(x) mean(posterior$alpha) + mean(posterior$beta_w) * x)

fit.plot