LangCog RStanArm Tutorial

Preamble

The goal of this notebook is to bring together resources on the rstanarm R package, which provides an interface between the R programming environment and the Stan probabilistic programming language.

The creators of rstanarm say the goal of the package is to,

make Bayesian estimation routine for the most common regression models that applied researchers use. This will enable researchers to avoid the counter-intuitiveness of the frequentist approach to probability and statistics with only minimal changes to their existing R scripts.

If you are already fitting models using lm() or glm() form the stats package or using the mixed effects versions lmer() and glmer() from lme4, then it will be easy to translate those models to a Bayesian framework to get the Bayesian benefits. To do so, you:

• add stan_ to your models
• specify a prior distribution (or use the default, weakly informative priors)
• learn how to “inspect” the posterior distribution and report the results

A great place to get started is the CRAN rstanarm page, which has a ton of helpful vignettes. This document is my attempt to pull together pieces of information that I found helpful. In it, we will:

• Fit different kinds of linear models to a simulated dataset: MLE, LMM, and Bayesian LMM
• Learn how to extract posterior samples from a fitted model object
• See a set of tools for visualizing distributions (both priors and posteriors) in a Bayesian analysis

Set up

Clear the workspace and set some global options for the notebook.

Load packages and set plotting aesthetics.

Loading required package: Rcpp
rstanarm (Version 2.15.3, packaged: 2017-04-29 06:18:44 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())
Loading required package: Matrix
Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages --------------------------------------------------------------------------
expand(): tidyr, Matrix
filter(): dplyr, stats
lag():    dplyr, stats

Attaching package: ‘magrittr’

The following object is masked from ‘package:purrr’:

    set_names

The following object is masked from ‘package:tidyr’:

    extract

A lot of the code in this tutorial is inspired by a blogpost from Tristan Mahr where he models the effect of sleep deprivation on Reaction Time using the sleepstudy dataset. The post is great! And it’s a really nice example of how to use visualization to understand what’s going on in a data analysis model.

Here is some backgound information on the dataset:

The average reaction time per day for subjects in a sleep deprivation study. On day 0 the subjects had their normal amount of sleep. Starting that night they were restricted to 3 hours of sleep per night. The observations represent the average reaction time on a series of tests given each day to each subject.

Load the data and look at the structure.

Observations: 180
Variables: 3
$ Reaction <dbl> 249.5600, 258.7047, 250.8006, 321.4398, 356.8519, 414.6901, 382.2038, 290.1486, 4...
$ Days     <dbl> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, ...
$ Subject  <fctr> 308, 308, 308, 308, 308, 308, 308, 308, 308, 308, 309, 309, 309, 309, 309, 309, ...

Let’s get a sense of the variability due to random effects by plotting the relationship between days of sleep deprivation and Reaction Time for each participant.

We can also summarize these data in table by computing the mean and standard deviation for each day. These summary statistics will also be used as parameters in our simulated dataset.

Days	m	stdev
0	256.65	32.13
1	264.50	33.43
2	265.36	29.47
3	282.99	38.86
4	288.65	42.54
5	308.52	51.77
6	312.18	63.17
7	318.75	50.10
8	336.63	60.20
9	350.85	66.99

Simulation

To help us understand the models in the RStanArm package, let’s write a few functions to simulate these data. Note the repeated measures structure: that we have multiple measures for each participant that come from different days. Participants could vary in their baseline RTs and they could also vary in how they respond to sleep deprivation, so we want to build this into our simulation.

Next, let’s write a function to generate model parameters for each participant. Note that the user can control:

• the amount of participant-level variability using the ss_noise argument, which controls the width of the noise distribution. higher numbers result in more noise for each participant’s measurement.
• the size of the effect of sleep deprivation using the location_shift argument. higher numbers result in a smaller effect of day

Test our simulate_participant() function with no added noise and the same age effect as in the original dataset.

Now that we can simulate the data for a single participant, let’s simulate the entire experiment.

Plot the simulated data by participant with regression lines fit to each participants’ data.

Fit a complete pooling that doesn’t know that our measurements come from different participants.

Call:
lm(formula = RT ~ day, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-155.24  -56.75  -14.16   16.20  591.62 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  246.441     13.464  18.304  < 2e-16 ***
day           12.088      2.522   4.793 3.22e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 102.4 on 198 degrees of freedom
Multiple R-squared:  0.104, Adjusted R-squared:  0.09944 
F-statistic: 22.97 on 1 and 198 DF,  p-value: 3.218e-06

Visualize variability by participant

Fit mixed effect model allowing for each participant to have a different intercept and slope.

Linear mixed model fit by REML ['lmerMod']
Formula: RT ~ 1 + day + (1 + day | id)
   Data: d

REML criterion at convergence: 2407.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5153 -0.5539 -0.1383  0.1581  5.7750 

Random effects:
 Groups   Name        Variance  Std.Dev.  Corr
 id       (Intercept) 0.000e+00 0.000e+00     
          day         2.610e-15 5.109e-08  NaN
 Residual             1.049e+04 1.024e+02     
Number of obs: 200, groups:  id, 20

Fixed effects:
            Estimate Std. Error t value
(Intercept)  246.441     13.464  18.304
day           12.088      2.522   4.793

Correlation of Fixed Effects:
    (Intr)
day -0.843

We can plot the model outpout for each participant alongside their data.

It looks like when we increase the level of noise for each participant, or between-participants variability, we get different inferences from lm vs. lmer(). In the complete pooling model the slope parameter is sig, but in the mixed-effects model it is not.

Fit Bayesian Mixed Effects model using RStanArm

Set the number of cores to the number of cores on your computer.

Fit the varying intercepts and slopes model.

starting worker pid=608 on localhost:11919 at 12:08:07.334
starting worker pid=616 on localhost:11919 at 12:08:07.613

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).

Gradient evaluation took 0.000734 seconds
1000 transitions using 10 leapfrog steps per transition would take 7.34 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).

Gradient evaluation took 0.000307 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.07 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 10.8088 seconds (Warm-up)
               3.52982 seconds (Sampling)
               14.3386 seconds (Total)

Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 15.1694 seconds (Warm-up)
               5.29238 seconds (Sampling)
               20.4618 seconds (Total)

We can use the default plot function in R directly on our fitted model to visualize the posterior interval over various parameters in the model.

Under the hood, the plot method is using functions from the bayesplot library. There are a lot of other visualizations of the posterior distributions that you can make using this library, and the nice thing is that all of the functions return ggplot objects!

For more information on what you can do with the bayesplot library, see the manual and this nice vignette.

This is bayesplot version 1.2.0

We can also create the same interval plot from before.

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing
scale.

Or if you want to make a density plot.

The bayesplot and default plots are nice, but if you want more flexibility, you can sample from samples from the posterior distributions over the slope and intercept parameters. This first visualization is a countour plot using code from Tristan Mahr’s lovely blogpost.

Let’s first extract the samples from the model obect.

Now make the contour plot, showing a set of plausible values for the intercept and slope parameters.

Diagnostics and plots using Shiny Stan

A few things to look at in Shiny Stan:

• traceplots of MCMC chains
• Posterior predictive (PP) checks
• distribution of observed data vs model-generated data
• distributions of test statistics
• generating tables

Inspecting Priors

It’s always a good idea to check the prior distribution over the parameters in your model to see what kinds of values are generated and confirm that our prior information produces sensible model behavior.

We can use the prior_summary() function to get information about the priors and adjustments used in the model.

Priors for model 'm_bglmer' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 5)
     **adjusted scale = 539.76

Coefficients
 ~ normal(location = 0, scale = 2)
     **adjusted scale = 74.98

Auxiliary (sigma)
 ~ half-cauchy(location = 0, scale = 5)
     **adjusted scale = 539.76

Covariance
 ~ decov(reg. = 2, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details

rstanarm comes with a built-in function for visualizing “belief change” after seeing the data.

Drawing from prior...

starting worker pid=642 on localhost:11919 at 12:08:42.835
starting worker pid=650 on localhost:11919 at 12:08:43.079

Gradient evaluation took 0.000354 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.54 seconds.
Adjust your expectations accordingly!



Gradient evaluation took 0.000375 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.75 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 5.16904 seconds (Warm-up)
               3.61552 seconds (Sampling)
               8.78456 seconds (Total)


 Elapsed Time: 5.3337 seconds (Warm-up)
               3.65644 seconds (Sampling)
               8.99013 seconds (Total)

But we can also sample from the prior directly using: prior_PD = TRUE.

starting worker pid=665 on localhost:11919 at 12:08:59.739
starting worker pid=673 on localhost:11919 at 12:08:59.955

Model Info:

 function:  stan_glmer
 family:    gaussian [identity]
 formula:   RT ~ day + (day | id)
 algorithm: sampling
 priors:    see help('prior_summary')
 sample:    2000 (posterior sample size)
 num obs:   200
 groups:    id (20)

Estimates:
              mean   sd     10%    50%    90% 
(Intercept)   -4.4  627.9 -796.5   -5.8  811.6
day           -1.5   75.5 -100.8   -3.1   96.9

Diagnostics:
            mcse Rhat n_eff
(Intercept) 14.0  1.0 2000 
day          1.7  1.0 2000 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

We can see that our initial prior of normal(0, 2) for the effect of sleep deprivation is centered around zero and thinks that an increase of +/- 76.17 milliseconds per day is reasonable. This seems sensible to me, but if you had more information about what you expect the range of effects to be, you could build this into the prior.

Let’s simulate the same model with a “wider” prior.

starting worker pid=688 on localhost:11919 at 12:09:16.065
starting worker pid=696 on localhost:11919 at 12:09:16.327

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).

Gradient evaluation took 0.000314 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.14 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).

Gradient evaluation took 0.000277 seconds
1000 transitions using 10 leapfrog steps per transition would take 2.77 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 5.28373 seconds (Warm-up)
               3.50108 seconds (Sampling)
               8.78481 seconds (Total)

Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 5.0033 seconds (Warm-up)
               3.3543 seconds (Sampling)
               8.35759 seconds (Total)


Model Info:

 function:  stan_glmer
 family:    gaussian [identity]
 formula:   RT ~ day + (day | id)
 algorithm: sampling
 priors:    see help('prior_summary')
 sample:    2000 (posterior sample size)
 num obs:   200
 groups:    id (20)

Estimates:
              mean    sd      10%     50%     90%  
(Intercept)   -16.6  1691.4 -2129.1   -12.9  2094.3
day             2.8   361.1  -458.6     2.7   464.6

Diagnostics:
            mcse Rhat n_eff
(Intercept) 37.8  1.0 2000 
day          8.1  1.0 2000 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Now we can see that the model considers values +/- 150 ms per day to be reasonable (within 1 sd), which is probably too broad.

A note on priors from the developers of rstanarm (priors vignette):

With very few exceptions, the default priors in rstanarm —the priors used if the arguments in the tables above are untouched— are not flat priors. Rather, the defaults are intended to be weakly informative. That is, they are designed to provide moderate regularization and help stabilize computation. For many (if not most) applications the defaults will perform well, but this is not guaranteed (there are no default priors that make sense for every possible model specification).

Because the scaling is based on the scales of the predictors (and possibly the outcome) these are technically data-dependent priors. However, since these priors are quite wide (and in most cases rather conservative), the amount of information used is weak and mainly takes into account the order of magnitude of the variables. This enables rstanarm to offer defaults that are reasonable for many models.

To disable automatic rescaling simply set the autoscale argument to to FALSE. For example:

starting worker pid=713 on localhost:11919 at 12:09:31.836
starting worker pid=721 on localhost:11919 at 12:09:32.074

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).

Gradient evaluation took 0.0003 seconds
1000 transitions using 10 leapfrog steps per transition would take 3 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).

Gradient evaluation took 0.000424 seconds
1000 transitions using 10 leapfrog steps per transition would take 4.24 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 3.30799 seconds (Warm-up)
               3.13123 seconds (Sampling)
               6.43922 seconds (Total)

Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 3.61395 seconds (Warm-up)
               3.20448 seconds (Sampling)
               6.81843 seconds (Total)

Priors for model 'test_no_autoscale' 
------
Intercept (after predictors centered)
 ~ student_t(df = 4, location = 0, scale = 10)

Coefficients
 ~ normal(location = 0, scale = 5)

Auxiliary (sigma)
 ~ exponential(rate = 0.1)

Covariance
 ~ decov(reg. = 2, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details

But the rstanarm developers point out that:

Disabling prior scale adjustments is usually unnecessary but is useful for when more informative prior information is available. There is an example of specifying an informative prior later in this vignette.