Regression in the Bayesian Framework

R Packages Needed

• To follow today’s lecture, please load the following packages:

library(tidyverse)
library(bayesrules)
library(bayesplot)
library(tidybayes)
library(ggpubr)
library(rstanarm)
library(janitor)
library(broom.mixed)

Introduction

• Before today, we effectively were performing one-sample tests of means, proportions, and counts.

• Now, we will focus on incorporating just a single predictor into our analysis.

  • Overall question: what is the relationship between Y (outcome) and X (predictor)?

• We now switch to thinking about analysis in terms of regression.

y_i = \beta_0 + \beta_1 x_i + \varepsilon_i

\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i

Working Example

• Capital Bikeshare is a bike sharing service in the Washington, D.C. area. To best serve its registered members, the company must understand the demand for its service. We will analyze the number of rides taken on a random sample of n days, (Y_1, Y_2, ..., Y_n).

• Because Y_i is a count variable, you might assume that ridership might need the Poisson distribution. However, past bike riding seasons have exhibited bell-shaped daily ridership with a variability in ridership that far exceeds the typical ridership.

  • (i.e., the Poisson assumption of \mu = \sigma does not hold here)

• We will instead assume that the number of rides varies normally around some typical ridership, \mu, with standard deviation, \sigma.

Y_i|\mu,\sigma \overset{\text{ind}}{\sim} N(\mu, \sigma^2)

Working Example

• In our example, Y is the number of rides and X is the temperature.

• Our specific goal will be to model the relationship between ridership and temperature:

  • Does ridership tend to increase on warmer days?
  • If so, by how much?
  • How strong is this relationship?

• It is reasonable to assume a positive relationship between temperature and number of rides.

  • As it warms up outside, folks are more likely to puruse outdoor activities, including biking.

• For learning purposes, let’s focus on the model with a single predictor (y ~ x in R).

Building the Regression Model

• Suppose we have n data pairs,

\{ (Y_1, X_1), (Y_2, X_2), ..., (Y_n, X_n) \}

• where Y_i is the number of rides and X_i is the high temperature (^oF) on day i.

• Assuming that the relationship is linear, we can model

\mu_i = \beta_0 + \beta_1 X_i,

• where \beta_0 and \beta_1 are model coefficients.

Building the Regression Model

• What do we mean by “model coefficients?”

\mu_i = \beta_0 + \beta_1 X_i,

• \beta_0 is the baseline for where our model crosses the y-axis, i.e., when X_i=0.
  • Is this meaningful when we are talking about the average ridership when it is 0^oF in DC?

Building the Regression Model

• What do we mean by “model coefficients?”

\mu_i = \beta_0 + \beta_1 X_i,

• \beta_1 is the slope, or average change, in the outcome (Y) for a one unit increase in the predictor (X).
  • Interpretation: for a [1 unit of predictor] increase in [the predictor], [the outcome] [increases or decreases] by [abs(\beta_1)].
  • In our example, suppose \beta_1=4.5. For a 1^oF increase in the temperature, the ridership increases by 4.5 riders.

Building the Regression Model

• We are now interested in the model

Y_i|\beta_0,\beta_1,\sigma\overset{\text{ind}}{\sim} N(\mu_i, \sigma^2) \text{ with } \mu_i = \beta_0 + \beta_1 X_i

• Note that \mu_i is the local mean (for a specific value of X).
  • In our data, \mu_i is the mean ridership for day i.
• The global mean (regardless of the value of X) is given by \mu.
  • In our data, \mu is the mean ridership, regardless of day.
• Under this model, \sigma is now measuring the variability from the local mean.

Building the Regression Model

• We know the assumptions for linear regression in the frequentist framework.

• In Bayesian Normal regression,

  • The observations are independent.
  • Y can be written as a linear function of X, \mu = \beta_0 + \beta_1 X.
  • For any value of X, Y varies normally around \mu with constant standard deviation \sigma.

Building the Regression Model

• First, we will assume that our prior models of \beta_0, \beta_1, and \sigma are independent.

• It is common to use the a normal prior for \beta_0 and \beta_1.

\begin{align*} \beta_0 &\sim N(m_0, s_0^2) \\ \beta_1 &\sim N(m_1, s_1^2) \end{align*}

• Then, it is the default to use the exponential for \sigma; both are restricted to positive values.

\sigma \sim \text{Exp}(l)

• Note that E[\sigma] = 1/l and \text{sd}[\sigma] = 1/l.

Building the Regression Model

• Thus, our regression model has the following formulation:

\begin{align*} &\text{data}: & Y_i|\beta_0, \beta_1, \sigma & \overset{\text{ind}}{\sim} N(\mu_i, \sigma^2) \text{ with } \mu_i = \beta_0 + \beta_1 X_i \\ &\text{priors:} &\beta_0 & \sim N(m_0, s_0^2) \\ & & \beta_1 &\sim N(m_1, s_1^2) \\ & & & \sigma \sim \text{Exp}(l) \end{align*}

Tuning the Prior Models

• Based on the past bikeshare analyses, we have the following centered prior understandings. What should our priors be?

1. On an average temperature day for DC (65^oF-70^oF), there are typically around 5000 riders, but this could vary between 3000 and 7000 riders.

2. For every one degree increase in temperature, ridership typically increases by 100 rides, but this could vary between 20 and 180 rides.

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

Tuning the Prior Models

1. On an average temperature day for DC (65^oF-70^oF), there are typically around 5000 riders, but this could vary between 3000 and 7000 riders.

\beta_{0\text{c}} \sim N(5000, 1000^2)

2. For every one degree increase in temperature, ridership typically increases by 100 rides, but this could vary between 20 and 180 rides.

\beta_{1\text{c}} \sim N(100, 40^2)

Tuning the Prior Models

3. At any given temperature, the daily ridership will tend to vary with a moderate standard deviation of 1250 rides.
   • Recall, E[\sigma] = 1/l = 1250, so l = 1/1250 = 0.0008.

\sigma \sim \text{Exp}(0.0008)

• Therefore, we have,

\begin{align*} \beta_{0\text{c}} &\sim N(5000, 1000^2) \\ \beta_{1\text{c}} &\sim N(100, 40^2) \\ \sigma &\sim \text{Exp}(0.0008) \end{align*}

Tuning the Prior Models

\beta_{0\text{c}} \sim N(5000, 1000^2) \ \ \ \beta_{1\text{c}} \sim N(100, 40^2) \ \ \ \sigma \sim \text{Exp}(0.0008)

Tuning the Prior Models

\beta_{0\text{c}} \sim N(5000, 1000^2) \ \ \ \beta_{1\text{c}} \sim N(100, 40^2) \ \ \ \sigma \sim \text{Exp}(0.0008)

Simulating the Posterior

• Let’s now update what we know using the bikes data.

• Looking at the basic relationship between the number of rides vs. the temperature (as it feels outside),

Simulating the Posterior

• Note: I am skipping the derivation / the true math behind how to find the posterior in this situation. (It involves multiple integrals!)

• We will use the stan_glm() function from rstanarm – it contains pre-defined Bayesian regression models.

  • stan_glm() also applies to the wider family of GzLM (i.e., logistic & Poisson/negbin).

bike_model <- stan_glm(rides ~ temp_feel, # data model
                       data = bikes, # dataset
                       family = gaussian, # distribution to apply
                       prior_intercept = normal(5000, 1000), # b0_c
                       prior = normal(100, 40), # b1_c
                       prior_aux = exponential(0.0008), # sigma
                       chains = 4, # 4 chains
                       iter = 5000*2, # 10000 iterations - throw out first 5000
                       seed = 84735) # starting place in RNG

Simulating the Posterior

bike_model <- stan_glm(rides ~ temp_feel, # data model
                       data = bikes, # dataset
                       family = gaussian, # distribution to apply
                       prior_intercept = normal(5000, 1000), # b0_c
                       prior = normal(100, 40), # b1_c
                       prior_aux = exponential(0.0008), # sigma
                       chains = 4, # 4 chains
                       iter = 5000*2, # 10000 iterations - throw out first 5000
                       seed = 84735) # starting place in RNG

• The first three arguments specify the structure of our data.
  • We want to state our model (rides ~ temp_feel),
  • specify the dataset (data = bikes), and
  • specify the family for the outcome (family = gaussian).

Simulating the Posterior

bike_model <- stan_glm(rides ~ temp_feel, # data model
                       data = bikes, # dataset
                       family = gaussian, # distribution to apply
                       prior_intercept = normal(5000, 1000), # b0_c
                       prior = normal(100, 40), # b1_c
                       prior_aux = exponential(0.0008), # sigma
                       chains = 4, # 4 chains
                       iter = 5000*2, # 10000 iterations - throw out first 5000
                       seed = 84735) # starting place in RNG

• Then we define our priors:
  • prior_intercept \to \beta_{0\text{c}},
  • prior \to \beta_{1\text{c}}, and
  • prior_aux \to \sigma.

Simulating the Posterior

bike_model <- stan_glm(rides ~ temp_feel, # data model
                       data = bikes, # dataset
                       family = gaussian, # distribution to apply
                       prior_intercept = normal(5000, 1000), # b0_c
                       prior = normal(100, 40), # b1_c
                       prior_aux = exponential(0.0008), # sigma
                       chains = 4, # 4 chains
                       iter = 5000*2, # 10000 iterations - throw out first 5000
                       seed = 84735) # starting place in RNG

• Finally, we specify Markov chain information:
  • Number of chains (chains),
  • number of iterations (iter) for each seed, and
  • the starting place in the RNG (seed).

Simulating the Posterior

• Wait, how does this work when we are looking at three model parameters?

• We will have three vectors: one for each model parameter.

\begin{align*} <\beta_0^{(1)}, & \ \beta_0^{(2)}, ..., \beta_0^{(5000)}> \\ <\beta_1^{(1)}, & \ \beta_1^{(2)}, ..., \beta_1^{(5000)} > \\ <\sigma^{(1)}, & \ \sigma^{(2)}, ..., \sigma^{(5000)} > \\ \end{align*}

Simulation Diagnostics

• Run the simulation code (from a previous slide).

• Then, run diagnostics:

neff_ratio(bike_model)
rhat(bike_model)
mcmc_trace(bike_model, size = 0.1)
mcmc_dens_overlay(bike_model)

Simulation Diagnostics

• Diagnostics:

neff_ratio(bike_model)

(Intercept)   temp_feel       sigma 
    0.98500     0.98525     0.96260 

rhat(bike_model)

(Intercept)   temp_feel       sigma 
  1.0001302   1.0001308   0.9998656 

• Quick diagnostics indicate that the resulting chains are trustworthy.

• The effective sample size ratios are slightly above 1 and the R-hat values are very close to 1.

Simulation Diagnostics

• Diagnostics:

• We can see that the chains are stable, mixing quickly, and behaving much like an independent sample.

Simulation Diagnostics

• Diagnostics:

• The density plot lets us visualize and examine the posterior models for each of our regression parameters, \beta_0, \beta_1, \sigma.

Interpreting the Posterior

• Okay… what does this mean, though?

tidy(bike_model, effects = c("fixed"), 
     conf.int = TRUE, 
     conf.level = 0.95)

• Thus, the posterior median relationship is

y = -2189.94 + 82.16x

• For a 1 degree increase in temperature, we expect ridership to increase by about 82 rides, with 95% credible interval (72.0, 92.0).

Interpreting the Posterior

• We can look at alternatives by drawing from the simulated data in bike_model.
  • The add_fitted_draws() function is from the tidybayes package.

bikes %>%
  add_fitted_draws(bike_model, n = 50) %>%
  ggplot(aes(x = temp_feel, y = rides)) +
  geom_line(aes(y = .value, group = .draw), alpha = 0.25, color = "blue") + 
  geom_point(data = bikes, size = .5) +
  theme_minimal()

Interpreting the Posterior

• We can look at alternatives by drawing from the simulated data in bike_model.

• We see that the plausible models are not super variable – we are more confident about the relationship we’re observing.

Interpreting the Posterior

• How does this compare against the frequentist version of regression?

Working Example

• We will examine Australian weather from the weather_WU data in the bayesrules package.
  • This data contains 100 days of weather data for each of two Australian cities: Uluru and Wollongong.

Working Example

• Let’s keep only the variables on afternoon temperatures (temp3pm) and a subset of possible predictors that we’d have access to in the morning:

weather_WU <- weather_WU %>% 
  select(location, windspeed9am, humidity9am, pressure9am, temp9am, temp3pm)
head(weather_WU)

Working Example

• Examining temperature,

Working Example

• Let’s model the 3 pm temperature as a function of the 9 am temperature on a given day i.
  • Outcome: Y_i = 3 pm temp
  • Predictor: X_{i1} = 9 am temp
• Then, we can model it using the Bayesian normal regression model,

\begin{align*} Y_i | \beta_0, \beta_1, \sigma &\overset{\text{ind}}{\sim} N(\mu_i, \sigma^2), \text{ with } \mu_i = \beta_0 + \beta_1 X_{i1} \\ \beta_{0c} &\sim N(25, 5^2) \\ \beta_1 &\sim N(0,3.1^2) \\ \sigma & \sim \text{Exp}(0.13) \end{align*}

• Note that we are using the centered intercept as 0 degree mornings are rare in Australia.

Working Example

• Simulating this model,

weather_model_1 <- stan_glm(temp3pm ~ temp9am, 
                            data = weather_WU, 
                            family = gaussian,
                            prior_intercept = normal(25, 5, autoscale = TRUE),
                            prior = normal(0, 2.5, autoscale = TRUE),
                            prior_aux = exponential(1, autoscale = TRUE),
                            chains = 4, 
                            iter = 5000*2, 
                            seed = 84735)

Working Example

• Note that we asked stan_glm() to autoscale our priors. What did it change them to?

prior_summary(weather_model_1)

Priors for model 'weather_model_1' 
------
Intercept (after predictors centered)
  Specified prior:
    ~ normal(location = 25, scale = 5)
  Adjusted prior:
    ~ normal(location = 25, scale = 38)

Coefficients
  Specified prior:
    ~ normal(location = 0, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 0, scale = 3.1)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.13)
------
See help('prior_summary.stanreg') for more details

Working Example

Working Example

Working Example

posterior_interval(weather_model_1, prob = 0.95)

                 2.5%    97.5%
(Intercept) 2.2619956 6.149827
temp9am     0.9455295 1.137758
sigma       3.7552460 4.568298

• 95% credible interval for \beta_1: (0.95, 1.14)

• 95% credible interval for \sigma: (3.76, 4.57)

Working Example

Categorical Predictors

• What if we look at the data by location?

ggplot(weather_WU, aes(x = temp3pm, fill = location)) + 
  geom_density(alpha = 0.5) + 
  theme_minimal()

• We should probably look at location as a predictor…

Categorical Predictors

• Let’s let X_{i2} be an indicator for the location,

X_{i2} = \begin{cases} 1 & \text{Wollongong} \\ 0 & \text{otherwise (i.e., Uluru).} \end{cases}

• We are treating “not-Wollongong” as our reference group – in this case, it is Uluru.

\begin{array}{rl} \text{data:} & Y_i \mid \beta_0, \beta_1, \sigma \overset{\text{ind}}{\sim} N(\mu_i, \sigma^2) \quad \text{with} \quad \mu_i = \beta_0 + \beta_1 X_{i2} \\ \text{priors:} & \beta_{0c} \sim N(25, 5^2) \\ & \beta_1 \sim N(0, 38^2) \\ & \sigma \sim \text{Exp}(0.13). \end{array}

Categorical Predictors

• Let’s think about our model.

y = \beta_0 + \beta_1 x_2

• What do our coefficients mean?
  • \beta_0 is the typical 3 pm temperature in Uluru (x_2=0).
  • \beta_1 is the typical difference in 3 pm temperature in Wollongong (x_2=1) as compared to Uluru (x_2=0).
  • \sigma represents the standard deviation in 3 pm temperatures in Wollongong and Uluru.

Categorical Predictors

• Let’s think about our model.

y = \beta_0 + \beta_1 x_2

• \beta_1 is the typical difference in 3 pm temperature in Wollongong (x_2=1) as compared to Uluru (x_2=0).

• When we use a binary predictor, this results in two models, effectively.

  • One when x_1=0 (Uluru) and one when x_1=1 (Wollongong).

\begin{align*} y = \beta_0 + \beta_1 x_2 \to \text{(U)} \ y &= \beta_0 \\ \text{(W)} \ y&= \beta_0+\beta_1 \end{align*}

Categorical Predictors

• Let’s simulate our posterior using weakly informative priors,

weather_model_2 <- stan_glm(temp3pm ~ location,
                            data = weather_WU,
                            family = gaussian,
                            prior_intercept = normal(25, 5, autoscale = TRUE),
                            prior = normal(0, 2.5, autoscale = TRUE), 
                            prior_aux = exponential(1, autoscale = TRUE),
                            chains = 4, 
                            iter = 5000*2, 
                            seed = 84735)

Categorical Predictors

mcmc_trace(weather_model_2, size = 0.1) + theme_minimal()

Categorical Predictors

mcmc_dens_overlay(weather_model_2) + theme_minimal()

Categorical Predictors

tidy(weather_model_2, effects = c("fixed"),
     conf.int = TRUE, conf.level = 0.95) 

• The estimated model is as follows,

\hat{y} = 29.72 - 10.31 x

• There is a 10 degree difference between the two cities.

Categorical Predictors

• We can also look at the temperatures by location,

Multiple Predictors

• What if we want to include multiple predictors?
  • Notice in the code, our model now has multiple predictors (temp9am and location).

weather_model_3 <- stan_glm(temp3pm ~ temp9am + location,
                            data = weather_WU, family = gaussian,
                            prior_intercept = normal(25, 5, autoscale = TRUE),
                            prior = normal(0, 2.5, autoscale = TRUE),
                            prior_aux = exponential(1, autoscale = TRUE),
                            chains = 4, 
                            iter = 5000*2, 
                            seed = 84735)

Multiple Predictors

• The simulation results in 20,000 posterior plausible relationships between temperature and location.

Multiple Predictors

tidy(weather_model_3, effects = c("fixed"),
     conf.int = TRUE, conf.level = 0.95) 

\hat{y} = 11.33 + 0.86 x_1 - 7.06 x_2

• As the 9 am temperature increases by 1 degree, the expected 3 pm temperature increases by 0.86 degrees (95% CI 0.80, 0.92).

• The 3 pm temperature is about 7 degrees lower in Wollongong than in Uluru (95% CI -7.75, -6.36).

Multiple Predictors

• Looking at our posterior predictive models,

Evaluating the Model

• We will now learn how to evaluate how good our regression model is.

• How fair is the model?

  • How was the data collected?
  • By whom and for what purpose?
  • How might the results of the data collection or analysis impact individuals and society?
  • What biases might be baked into the analysis?

Evaluating the Model

• We want to now learn how to evaluate how good our regression model is.

• How wrong is the model?

  • Are our assumptions reasonable?
  • Is the predictor set appropriate and complete?

• How accurate are the posterior predictive models?

  • How far are the posterior predictive models from reality?

Working Example

• Consider our example before the break, the Capital Bikeshare data.

• We simulated under the following model:

\begin{align*} Y_i | \beta_0, \beta_1, \sigma &\overset{\text{ind}}{\sim} N(\beta_0 + \beta_1X_i,\sigma^2) \\ \beta_{0\text{c}} &\sim N(5000, 1000^2) \\ \beta_{1\text{c}} &\sim N(100, 40^2) \\ \sigma &\sim \text{Exp}(0.0008) \end{align*}

Is the Model Fair?

• How was the data collected?

  • It is reasonable to assume that the bikeshare was collected electronically.
  • We do not have individual-level information, only aggregated ride data.

• By whom and for what purpose was the data collected?

  • It is reasonable for a bike sharing company to track data such as number of rides in a day.
  • The company is likely doing this to make evidence-based business decisions.

Is the Model Fair?

• How might the results of the analysis, or the data collection itself, impact individuals and society?

  • Because there is no individual-level information, this will not negatively impact customers.
  • Improving the service might help reduce car traffic in D.C. - this is good!

• Note!

  • Just because we cannot imagine societal impacts does not mean that none exist.
  • It’s critical to recognize our perspective based on our experiences.
  • How can we truly determine the impacts?
    • We would need to engage with the company and its community.

Is the Model Fair?

• What biases might be baked into this analysis?

  • The data is only for one company.
    • Other bikeshare companies may have different usage patterns.
  • The data is only for one city.
    • Other cities may have different usage patterns for bikeshare.
    • Other cities will have different weather patterns.
  • How generalizable is our model?

Is the Model Fair?

• How do I approach thinking about bias in analysis?

  • Who is and is not represented in the data?
  • What assumptions were made about the data?
  • What assumptions are made about the truth/the population?
  • Is the analysis method appropriate for the question posed?
    • e.g., did they use logistic regression for a binary outcome?
  • Were there any “reaching” assumptions made in order to complete the analysis?
    • e.g., assuming data is MCAR when it is likely MAR or MNAR.
    • e.g., assuming independent observations when there is likely correlation.
  • Are the results stated appropriately? (vs. overstated)
    • Are broad, sweeping generalizations made from a narrow dataset?

How Wrong is the Model?

How Wrong is the Model?

• Recall the data model,

Y_i | \beta_0, \beta_1, \sigma \overset{\text{ind}}{\sim} N(\beta_0 + \beta_1X_i,\sigma^2)

• What are the assumptions on our model?

  • The observed data, Y_i, is independent of other observed data, Y_j for i \neq j.
  • Y can be written as a linear function of X, \mu=\beta_0 + \beta_1 X.
  • At any X, Y varys normally around \mu with consistent variability, \sigma.

How Wrong is the Model?

• The observed data, Y_i, is independent of other observed data, Y_j for i \neq j.

  • We can technically test for this, but it is not often employed by statisticians.
  • Instead, we rely on knowledge of the data collection process.
  • If the data collection process suggests dependence, we should consider a different model.
    • e.g., a modeling approach that will account for the correlation between observations.

How Wrong is the Model?

• The observed data, Y_i, is independent of other observed data, Y_j for i \neq j.

• In this example, the data are daily counts of bike rides.

  • There is inherently some correlation – the ridership today is likely correlated to the ridership yesterday.
  • However, this can be explained by the time of year and features associated with the time of year, like temperature.
  • We should ask ourselves if it’s reasonable to assume that the temperature cancels out the time correlation.

• “Reasonable” does not mean “perfect,” but instead, “good enough to proceed.”

How Wrong is the Model?

• Returning to our bike data, we can look at the scatterplot of temperature and number of rides.

How Wrong is the Model?

• We can overlay a line of best fit using geom_smooth(method = "lm").

How Wrong is the Model?

• Or we can overlay a locally estimated scatterplot smoothing (LOESS) line using geom_smooth().

How Wrong is the Model?

• How different are the two lines?

How Wrong is the Model?

• When does “visualization” suggest that a linear model is not appropriate?

  • Does the LOESS line deviates substantially from the linear line?
  • Here, we see that the LOESS line is fairly close to the linear line.
  • Thus, it is reasonable to assume a linear relationship between X and Y.

• Are there “issues” with this approach?

  • When we have more than one predictor, yes.
  • We are only visualizing one predictor at a time.
  • We will need to use other tools to assess linearity in multiple regression.

How Wrong is the Model?

• To examine linearity in multiple regression, we will perform a posterior predictive check.

• If the combined model assumptions are reasonable, then our posterior model should be able to simulate ridership data that’s similar to the original 500 rides observations.

• For each of the 20,000 posterior plausible parameter sets, \{\beta_0, \beta_1, \sigma\}, we can predict ridership data from the observed temperature data.
  • This gives us 20,000 unique datasets of predicted ridership, each with n=500.

How Wrong is the Model?

• Recall our analysis,

bike_model <- stan_glm(rides ~ temp_feel, 
                       data = bikes, 
                       family = gaussian,
                       prior_intercept = normal(5000, 1000, autoscale = TRUE), 
                       prior = normal(100, 40, autoscale = TRUE), 
                       prior_aux = exponential(0.0008, autoscale = TRUE), # sigma
                       chains = 4, 
                       iter = 5000*2, 
                       seed = 84735) 

How Wrong is the Model?

• Recall our analysis,

tidy(bike_model, conf.int = TRUE) 

• This gives the model,

\hat{y} = -2177.9 + 81.9x \ \ \ \to \ \ \ \hat{\text{rides}} = -2177.9 + 81.9\text{temp}

How Wrong is the Model?

• How do we perform the posterior predictive check?

• The pp_check() function plots a number of the corresponding simulated datasets.

pp_check(bike_model, 
         nreps=50, 
         seed = 84735) + 
  xlab("Number of Rides") + 
  theme_minimal()

How Wrong is the Model?

• The pp_check() function plots a number of the corresponding simulated datasets.

How Accurate are Posterior Predictive Models?

• Although we are interested in predicting the future, we can begin by evaluating how well it predicts the data that was used to build the model.

• We can use the posterior_predict() function to generate posterior predictive summaries to either our data or new data.

predictions <- posterior_predict(bike_model, newdata=bikes, seed = 84735)
dim(predictions)

[1] 20000   500

• We have 20,000 posterior predictive datasets, each with 500 predicted ride counts.

How Accurate are Posterior Predictive Models?

• We can also use the ppc_intervals() function to visualize the posterior predictive intervals for each observation.

ppc_intervals(bikes$rides, 
              yrep = predictions, 
              x = bikes$temp_feel, 
              prob = 0.5, 
              prob_outer = 0.95) + 
  theme_minimal()

How Accurate are Posterior Predictive Models?

• We can also use the ppc_intervals() function to visualize the posterior predictive intervals for each observation.

How Accurate are Posterior Predictive Models?

• Let Y_1, Y_2, ..., Y_n denote n observed outcomes. Each Y_i has a corresponding posterior predictive model with mean Y_i' and standard deviation \text{sd}_i

• Median Absolute Error (MAE): the typical difference between the observed value and posterior predictive mean.

\text{MAE} = \text{median} | Y_i - Y_i'|

• Scaled Median Absolute Error (MAE scaled): the typical difference between the observed value and posterior predictive mean.

\text{MAE scaled} = \text{median} \frac{| Y_i - Y_i'|}{\text{sd}_i}

How Accurate are Posterior Predictive Models?

• Let Y_1, Y_2, ..., Y_n denote n observed outcomes. Each Y_i has a corresponding posterior predictive model with mean Y_i' and standard deviation \text{sd}_i

• Within 50: Proportion of observed values that fall within their 50% posterior prediction interval.

• Within 95: Proportion of observed values that fall within their 95% posterior prediction interval.

• We get these four values out of the prediction_summary() function,

set.seed(84735)
prediction_summary(bike_model, data = bikes)

How Accurate are Posterior Predictive Models?

• Let’s now interpret these values.

set.seed(84735)
prediction_summary(bike_model, data = bikes)

• The MAE of 989.7 means that the typical difference between the observed ride count and the posterior predictive mean is about 990 rides.
  • The scaled MAE of 0.77 means that the typical difference between the observed ride count and the posterior predictive mean is about 0.77 standard deviations.
• Only 43.8% of test observations fall within their respective 50% prediction interval, while 95.8% fall within their respective 95% prediction interval.

How Accurate are Posterior Predictive Models?

• Looking at these side by side,

How Accurate are Posterior Predictive Models?

• What is our final conclusion? Does our Bayesian model produce accurate predictions?

• Unfortunately, the answer is subjective and dependent on context.

  • Is an MAE of 990 rides acceptable for the business?
  • Is having only 43.8% of observations within the 50% prediction interval acceptable?

• Per usual, we are just the statisticians. The business must decide what is acceptable.

How Accurate are Posterior Predictive Models?

• How can we improve the model’s posterior predictive accuracy?

  • Add additional predictors or terms (e.g., interactions) to the model.
  • Consider different distributions for the outcome variable.
    • e.g., Poisson, negative binomial, binomial, etc.
      
  • Collect more data.
    • This is always the cheeky answer.

Wrap Up

• Today we covered Bayesian general linear regression and basic evaluation of regression models.
  • Is the model fair?
  • How wrong is the model?,
  • How accurate are the posterior predictive models?
• We as statisticians provide the numbers, but our collaborator decides what is acceptable in terms of “good enough” for the model’s predictive accuracy.