Bayesian Regression Models with RStanARM

Madison R Users Group

Github repository
@tjmahr

• How I got into Bayesian statistics
• Some intuition-building about Bayes theorem
• Tour of RStanARM
• Where to learn more about Bayesian statistics

A programming language for probablistic stats.

• Write out a description of the data and model using a mathy syntax.
• The model is compiled into an executable that does the sampling.
• Pystan and RStan are interfaces that help you send/receive data to/from Stan models.

# R version
lm(log(earnings) ~ height + male)

data {
  int<lower=0> N;
  vector[N] earn;
  vector[N] height;
  vector[N] male;
}
transformed data {
  // log transformation
  vector[N] log_earn;
  log_earn = log(earn);
}
parameters {
  vector[3] beta;
  real<lower=0> sigma;
}
model {
  log_earn ~ normal(beta[1] + beta[2] * height + beta[3] * male, sigma);
}
generated quantities {
  // optional
}

• Program is a description of the data and a model.
• Language wants a more formal information about the data (type, length, bounds), presumably for more efficient sampling.
• The tilde operator ~ is a sampling statement.

RStan Applied Regression Modeling

• Batteries-included, precompiled versions of common regression models.
• glm -> stan_glm, glmer -> stan_glmer.
• Write your regular code. Add stan_ to the front, and add a prior.
• CRAN page
  • Very good! They have lots of detailed vignettes!
• Proper successor to the arm package.

# Some toy data
davis <- car::Davis %>% filter(100 < height) %>% as_data_frame

davis
#> # A tibble: 199 × 5
#>       sex weight height repwt repht
#>    <fctr>  <int>  <int> <int> <int>
#> 1       M     77    182    77   180
#> 2       F     58    161    51   159
#> 3       F     53    161    54   158
#> 4       M     68    177    70   175
#> 5       F     59    157    59   155
#> 6       M     76    170    76   165
#> 7       M     76    167    77   165
#> 8       M     69    186    73   180
#> 9       M     71    178    71   175
#> 10      M     65    171    64   170
#> # ... with 189 more rows

model <- lm(weight ~ height * sex, davis)
summary(model) %>% print(digits = 2)
#> 
#> Call:
#> lm(formula = weight ~ height * sex, data = davis)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#>  -20.9   -4.8   -0.9    4.5   41.1 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   -45.71      22.19    -2.1     0.04 *  
#> height          0.62       0.13     4.6    7e-06 ***
#> sexM          -55.62      32.54    -1.7     0.09 .  
#> height:sexM     0.37       0.19     2.0     0.05 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 8 on 195 degrees of freedom
#> Multiple R-squared:  0.64,   Adjusted R-squared:  0.64 
#> F-statistic: 1.2e+02 on 3 and 195 DF,  p-value: <2e-16

library(rstanarm)
#> Loading required package: Rcpp
#> rstanarm (Version 2.12.1, packaged: 2016-09-12 13:08:24 UTC)
#> - Do not expect the default priors to remain the same in future rstanarm versions.
#> Thus, R scripts should specify priors explicitly, even if they are just the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores())

• So… hard-code the priors.
• I only set the mc.cores option when I have a big model that going to take more than a minute to fit.

stan_model
#> stan_glm(formula = weight ~ height * sex, family = gaussian, 
#>     data = davis, prior = normal(location = 0, scale = 5), prior_intercept = normal(0, 
#>         10))
#> 
#> Estimates:
#>             Median MAD_SD
#> (Intercept) -48.1   21.5 
#> height        0.6    0.1 
#> sexM        -49.9   29.8 
#> height:sexM   0.3    0.2 
#> sigma         8.0    0.4 
#> 
#> Sample avg. posterior predictive 
#> distribution of y (X = xbar):
#>          Median MAD_SD
#> mean_PPD 65.3    0.8  
#> 
#> Observations: 199  Number of unconstrained parameters: 5

• Split into estimation and diagnostic information
• mean_PPD is the predicted value for a completely average observation

Coerce to a data-frame. Columns are parameters. One row per posterior sample.

samples <- stan_model %>% as.data.frame %>% tbl_df
samples
#> # A tibble: 4,000 × 5
#>    `(Intercept)`    height       sexM `height:sexM`    sigma
#>            <dbl>     <dbl>      <dbl>         <dbl>    <dbl>
#> 1      -56.71142 0.6895427 -68.550433    0.44384885 7.537152
#> 2      -56.70521 0.6909141 -69.411803    0.44392578 7.784114
#> 3      -56.34580 0.6865143 -69.497981    0.44949406 7.795977
#> 4      -30.71059 0.5202990 -58.429948    0.40582612 7.991750
#> 5      -34.68467 0.5452743 -64.090664    0.42866842 7.784008
#> 6      -87.38731 0.8865953   2.711235    0.02298263 8.324237
#> 7      -98.27211 0.9459078  15.229156   -0.05056481 8.665427
#> 8      -58.71536 0.6977810 -13.289139    0.13100168 7.894166
#> 9      -63.19315 0.7205729 -22.734855    0.17883580 8.501026
#> 10     -71.94144 0.7791263 -52.550236    0.33906170 8.051117
#> # ... with 3,990 more rows

Basic plot() method shows 80% and 95% intervals.

plot(stan_model)

# Match a subset of the parameters
plot(stan_model, regex_pars = "height")

This thing is pretty awesome.

launch_shinystan(stan_model)

[Demo outside of slides.]

• Diagnostics
  • Mixing of chains
  • R hat, n eff
  • PPcheck
• Parameters plot
  • Different parameters selected
  • KDE
• Explore
  • Inspect one parameter in multiview
  • Bivariate view

Get a data-frame with the parameters from each sample. We now have 4000 plausble regression lines.

df_model <- stan_model %>% as_data_frame()
df_model
#> # A tibble: 4,000 × 5
#>    `(Intercept)`    height       sexM `height:sexM`    sigma
#>            <dbl>     <dbl>      <dbl>         <dbl>    <dbl>
#> 1      -56.71142 0.6895427 -68.550433    0.44384885 7.537152
#> 2      -56.70521 0.6909141 -69.411803    0.44392578 7.784114
#> 3      -56.34580 0.6865143 -69.497981    0.44949406 7.795977
#> 4      -30.71059 0.5202990 -58.429948    0.40582612 7.991750
#> 5      -34.68467 0.5452743 -64.090664    0.42866842 7.784008
#> 6      -87.38731 0.8865953   2.711235    0.02298263 8.324237
#> 7      -98.27211 0.9459078  15.229156   -0.05056481 8.665427
#> 8      -58.71536 0.6977810 -13.289139    0.13100168 7.894166
#> 9      -63.19315 0.7205729 -22.734855    0.17883580 8.501026
#> 10     -71.94144 0.7791263 -52.550236    0.33906170 8.051117
#> # ... with 3,990 more rows

Plot lines with the median parameter values and a random sample of lines.

fits <- sample_n(df_model2, 200)
medians <- df_model2 %>% summarise_each(funs = funs(median))

p2 <- ggplot(davis) +
  aes(x = height, y = weight, color = sex) +
  geom_abline(aes(color = "F", intercept = F_Intercept,
                  slope = F_Slope), data = fits, alpha = .075) +
  geom_abline(aes(color = "M", intercept = M_Intercept,
                  slope = M_Slope), data = fits, alpha = .075) +
  geom_abline(aes(color = "F", intercept = F_Intercept,
                  slope = F_Slope), data = medians, size = 1.25) +
  geom_abline(aes(color = "M", intercept = M_Intercept,
                  slope = M_Slope), data = medians, size = 1.25) +
  geom_point() +
  theme(legend.position = c(0, 1), legend.justification = c(0, 1))

• Depicts uncertainty by drawing line of best fit plus many other plausible regression lines.
• Visualization works best when there is only one group and we don't mind having regression lines stretch past the range of the data.

Get a sample of height values within each group's range.

# # do some stuff
# ...

new_data
#> # A tibble: 160 × 3
#>    Observation   sex   height
#>          <chr> <chr>    <dbl>
#> 1            1     F 148.0000
#> 2            2     F 148.3797
#> 3            3     F 148.7595
#> 4            4     F 149.1392
#> 5            5     F 149.5190
#> 6            6     F 149.8987
#> 7            7     F 150.2785
#> 8            8     F 150.6582
#> 9            9     F 151.0380
#> 10          10     F 151.4177
#> # ... with 150 more rows

With the normal predict() function on a classical regression model, we use the model parameters to get the expected mean (\( \mu_i \)) for each row in newdata.

With posterior_linpred(), we do the same thing, but for each of the 4000 posterior samples.

linpreds <- posterior_linpred(stan_model, newdata = new_data)
dim(linpreds)
#> [1] 4000  160

p3 <- ggplot(davis) +
  aes(x = height, y = weight, color = sex, group = sex) +
  geom_point() +
  geom_ribbon(aes(ymin = ymin, ymax = ymax, y = NULL, color = NULL),
              data = df_predictions, fill = "grey60", alpha = .4) +
  geom_line(aes(y = median), data = df_predictions, size = 1.25) +
  theme(legend.position = c(0, 1), legend.justification = c(0, 1))

RStanARM does not like posterior_linpred(). Documentation says:

See also: posterior_predict to draw from the posterior predictive distribution of the outcome, which is almost always preferable.

Why? Because we have samples of the error term sigma, and we're not using them!

Same plan as before, but use posterior_predict() to get predictions for observations (\( y_i \)), instead of estimates of the mean (\( \mu_i \)).

# ?posterior_predict
preds <- posterior_predict(stan_model, newdata = new_data)
dim(preds)
#> [1] 4000  160

p4 <- ggplot(davis) +
  aes(x = height, y = weight, color = sex, group = sex) +
  geom_point() +
  geom_ribbon(aes(ymin = ymin, ymax = ymax, y = NULL, color = NULL),
              data = df_predictions, fill = "grey60", alpha = .2) +
  geom_line(aes(y = median),
              data = df_predictions, size = 1.25) +
  theme(legend.position = c(0, 1), legend.justification = c(0, 1))

Do we need the sex-height interaction?

Do we need the sex predictor at all?

Let's fit alternative models.

stan_model_no_inter <- update(stan_model, weight ~ sex + height)
stan_model_no_group <- update(stan_model, weight ~ height)

For classical regression models, we could compare models with AIC and BIC.

model_no_inter <- update(model, weight ~ sex + height)
model_no_group <- update(model, weight ~ height)
model_list <- list(
  no_group = model_no_group,
  no_inter = model_no_inter,
  inter = model)

model_list %>% lapply(AIC) %>% unlist
#> no_group no_inter    inter 
#> 1421.572 1401.590 1399.692

model_list %>% lapply(BIC) %>% unlist
#> no_group no_inter    inter 
#> 1431.452 1414.764 1416.158

Or with anova().

anova(model_no_group, model_no_inter, model)
#> Analysis of Variance Table
#> 
#> Model 1: weight ~ height
#> Model 2: weight ~ sex + height
#> Model 3: weight ~ height * sex
#>   Res.Df   RSS Df Sum of Sq       F    Pr(>F)    
#> 1    197 14312                                   
#> 2    196 12815  1   1496.72 23.2250 2.892e-06 ***
#> 3    195 12567  1    248.64  3.8582   0.05092 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

For RStanARM, we compare the models by using approximation leave-one-out cross-validation. The details are outside of the scope of this talk.

loo1 <- loo(stan_model_no_group)
loo2 <- loo(stan_model_no_inter)
loo3 <- loo(stan_model)

• In our models, the prior normal(0, 5) gets rescaled to match each predictor variable.
• This rescaling is not well documented.
• Some options:
  • scale() our measurements so that they are scale free (with mean = 0 and SD = 1).
  • Sample from the prior to see what kinds of values are generated by our prior and confirm that our prior information yields sensible model behavior.

RStanARM will compare samples from the posterior and the prior.

posterior_vs_prior(stan_model)

We can also sample from the prior directly.

stan_model_prior <- stan_glm(
  weight ~ height * sex,
  data = davis,
  family = gaussian,
  prior = normal(0, 5),
  prior_intercept = normal(0, 10),
  # this line is new:
  prior_PD = TRUE
)

• How I got into Bayesian statistics
• Some intuition-building about Bayes theorem
• Tour of RStanARM
• Where to learn more about Bayesian statistics