The Hyperparameters

Okay, with the model set up, it is time to turn our attention to the hyperparameters. The bike shop has been operating for a while, so our clients do have quite a bit of prior information that they provide to us. They ask us to incorporate this prior information into our Bayesian regression model.

Here’s what we know:

1. On an average temperature day, there are typically around 5000 riders, though this average could be somewhere between 3000 and 7000.

2. For every increase in \(X_i^C\) by 1, ridership typically increases by 100 rides, though this average increase could be as low as 20 or as high as 180.

3. At any given temperature, daily ridership will tend to vary with a moderate standard deviation of 1250 rides.

Now, let’s do something we haven’t done before - let’s check on this prior specification and see if it actually reflects what we want it to. We can do this by running a simulation study. A simulation study just means that we generate data ourselves under a given model and then check to see what properties that data has. We have been given very specific prior beliefs - let’s make sure our prior is reflecting those beliefs well.

Let’s see how to do this.

One piece of prior information we were given is that is that at any given temperature, daily ridership will tend to vary with a moderate standard deviation of 1250 rides.

A second piece of prior information is that on an average temperature day, there are typically around 5000 riders, though this average could be somewhere between 3000 and 7000. This is going to require a little more work, as it discusses the variance in \(\beta_0\).

That seemed like a lot of work…why did we do all of this?? Well, it’s a really important tool that we use in Bayesian statistics to not only verify that our priors are calibrated (specified) correctly, but also to explore what different relationships different choices of hyperparameters can reflect. This is important when we don’t necessarily have strong prior beliefs. Our client was very kind to us today and gave us very specific information, but this is not always the case.

In practice, we would typically simulated more than 5 data sets. I would have done at least 100, probably more, to get a good idea of what information the priors were conveying.

The MCMC

With our model now set and our prior hyperparameters specified, it’s time to actually proceed with estimating the posterior distributions of our parameters!

For our model, we need to use MCMC to estimate the posterior distributions of our parameters. For the first time since we learned MCMC, we are not going to code it ourselves!! Instead, we are going to be using a very commonly used package called stan. To use this package, you need to load a few libraries:

# Allows us to run MCMC
library(rstanarm)
# Allows us to make pretty pictures with our MCMC draws
library(bayesplot)

The function we will be using is called stan_glm. Yes, glm like Generalized Linear Models, which is great news for us because that means this package can be used for things like Bayesian Logistic Regression or Bayesian Poisson Regression. Yay for learning a versatile function that can be applied in a variety of settings!!

So, let’s get to know Stan. The skeleton of the code you will need looks like this:

rental_model <- stan_glm( Y ~ X, data = ,
                       family = gaussian,
                       prior_intercept = normal(, ),
                       prior = normal(, ), 
                       prior_aux = exponential(),
                       chains = , iter = , seed = 84735)

Let’s talk through this.

• Y: this is your \(Y\) variable, so for us rides, which represents the number of bikes that are rented for rides on a given day.
• X : this is your CENTERED and SCALED \(X\) variable. This is important!!! This function will not center for you.
• family = gaussian: Gaussian is another word for a normal distribution, so this is telling the computer that we are using OLS regression. For different regression models, this would change.

Now we are into the prior specification part.

• prior_intercept: This is the prior on \(\beta_0\). We are specifying a normal, and you need to fill in the prior mean and standard deviation (not variance).
• prior: This is the prior on \(\beta_1\). We are specifying a normal, and you need to fill in the prior mean and standard deviation (not variance).
• prior_aux: This is the prior on \(\sigma\), and you need to feed it the prior hyperparameters \(\alpha\).

The last few arguments have to do with MCMC related things.

• chains: How many chains do we want to run? For today, let’s do 3.
• iter: This is \(S\), the number of iterations of MCMC we want to run.
• seed: This is the random seed we are all going to use today. This makes Dr. Dalzell’s life easier, so please and thank you for using the seed in the skeleton code.

Now, what you have created is a list object with lots of things in it. Let’s see how we can use our results.

Assessing Convergence

The first thing we probably want to do it to create a trace plot. Using the bayesplot package, we can do this quite easily.

mcmc_trace(rental_model, size = 0.1)

It is often also helpful to look at the empirical posterior distributions that we have created using our draws from the MCMC chains. To do that, we use the following code:

mcmc_dens_overlay(rental_model)

This shows us what the posterior distributions for each parameter looks like, based on our MCMC samples!

One other thing we may want to assess is effective sample size. Luckily, we can do this, too.

neff_ratio(rental_model)

The result running the code above is a little different from what we are used to when we assess effective sample size. Instead of returning a number, it returns a ratio. To get the effective sample size we are used to, we just multiply this by \(S\). Note: This code assumes we have only fed in post burn-in samples. If you want to burn any, remove them BEFORE running the code above.

If you want to look at the autocorrelation itself, we can do that, too! Let’s take a look for \(\beta_1\).

mcmc_acf(rental_model, pars = "XC")

Using the Results

Okay, we done a lot of things. However, one thing we have not done is answer our client’s question. Let’s do that now, shall we?

At the moment, we have a lot of samples for each of the three parameters. Is there is a neat way to summarize all of this?

library(broom.mixed)
tidy(rental_model, effects = c("fixed", "aux"),
     conf.int = TRUE, conf.level = )

The function above provides a very nice summary of our posterior results. The effects tells the computer which of the parameters we want to do posterior inference on. If you remember from the skeleton code, aux means the standard deviation \(\sigma\). The fixed means the rest (\(\beta_0, \beta_1\)).

You will also notice that we can request a credible interval. Yes it says confidence - I didn’t write this code!! It is a credible interval.

Turning in your assignment

When your Markdown document is complete, do a final knit and make sure everything compiles. You will then submit your document on Canvas. You must submit a PDF or html document to Canvas. No other formats will be accepted. Make sure you look through the final document to make sure everything has knit correctly.