Homework1

Bayesian Computation in Stan and R

Load/Setup STRAN

library(rstan) # observe startup messages

## Loading required package: StanHeaders

## Loading required package: ggplot2

## rstan (Version 2.21.2, 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)

options(mc.cores = parallel::detectCores()) # Use all local cores (does this work for all r packages?)
rstan_options(auto_write = TRUE) # saves stan to hard drive
Sys.setenv(TZ = "America/Toronto") # this deals with a time zone bug issue that was causes an error when running the model

Example 1: Eight Schools

This is an example in Section 5.5 of Gelman et al (2003), which studied coaching effects from eight schools. For simplicity, we call this example “eight schools.”

Write STAN file

We start by writing a Stan program for the model in a text file. If you are using RStudio version 1.2.x or greater, click on File -> New File -> Stan File . Otherwise, open your favorite text editor. Either way, paste in the following and save your work to a file called schools.stan in R’s working directory (which can be seen by executing getwd()) You can also embed plots, for example:

// saved as schools.stan
data {
  int<lower=0> J;         // number of schools 
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates 
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
  target += normal_lpdf(eta | 0, 1);       // prior log-density
  target += normal_lpdf(y | theta, sigma); // log-likelihood
}

Be sure that your Stan programs ends in a blank line without any characters including spaces and comments.

In this Stan program, we let theta be a transformation of mu, eta, and tau instead of declaring theta in the parameters block, which allows the sampler will run more efficiently .

Prepare data

We can prepare the data (which typically is a named list) in R with:

And we can get a fit with the following R command. Note that the argument to file = should point to where the file is on your file system unless you have put it in the working directory of R in which case the below will work.

fit <- stan(file = 'schools.stan', data = schools_dat)

## Warning: There were 9 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

##unexpected WARNINGS: Warning message: In utils::install.packages(“openssl”, repos = “https://cran.rstudio.com/”) : installation of package ‘openssl’ had non-zero exit status

Output

print(fit)

## Inference for Stan model: schools.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##            mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## mu         7.82    0.17 5.20  -2.11   4.48   7.71  10.81  17.90   938    1
## tau        6.29    0.15 5.41   0.15   2.22   4.93   8.84  20.39  1285    1
## eta[1]     0.38    0.01 0.95  -1.56  -0.25   0.40   1.02   2.15  3988    1
## eta[2]     0.00    0.01 0.91  -1.78  -0.60  -0.02   0.60   1.79  3844    1
## eta[3]    -0.18    0.01 0.90  -1.92  -0.79  -0.18   0.42   1.63  3694    1
## eta[4]    -0.02    0.02 0.87  -1.72  -0.61  -0.01   0.57   1.69  3247    1
## eta[5]    -0.32    0.02 0.87  -1.97  -0.89  -0.35   0.22   1.48  3341    1
## eta[6]    -0.19    0.02 0.91  -1.95  -0.79  -0.21   0.41   1.58  3247    1
## eta[7]     0.34    0.01 0.93  -1.61  -0.26   0.37   0.95   2.13  3832    1
## eta[8]     0.07    0.02 0.95  -1.82  -0.55   0.07   0.71   1.94  3109    1
## theta[1]  11.09    0.19 8.33  -2.30   5.59   9.83  15.20  31.27  1973    1
## theta[2]   7.84    0.09 6.23  -4.50   3.98   7.71  11.59  20.49  4680    1
## theta[3]   6.19    0.13 7.54 -11.68   2.30   6.76  10.71  20.05  3129    1
## theta[4]   7.57    0.11 6.43  -5.17   3.76   7.51  11.40  20.19  3623    1
## theta[5]   5.13    0.10 6.17  -8.43   1.56   5.59   9.32  15.88  3749    1
## theta[6]   6.21    0.11 6.66  -8.34   2.43   6.44  10.44  18.47  3820    1
## theta[7]  10.45    0.13 6.93  -1.70   5.84   9.86  14.48  25.89  2843    1
## theta[8]   8.41    0.25 8.19  -7.16   3.79   8.05  12.35  26.17  1048    1
## lp__     -39.71    0.08 2.68 -45.52 -41.42 -39.44 -37.77 -35.15  1262    1
## 
## Samples were drawn using NUTS(diag_e) at Sat Sep 25 06:48:00 2021.
## For each parameter, 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).

Plot Paramteres

plot(fit)

## 'pars' not specified. Showing first 10 parameters by default.

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

pairs(fit, pars = c("mu", "tau", "lp__"))

la <- extract(fit, permuted = TRUE) # return a list of arrays 
mu <- la$mu 

### return an array of three dimensions: iterations, chains, parameters 
a <- extract(fit, permuted = FALSE) 

### use S3 functions on stanfit objects
a2 <- as.array(fit)
m <- as.matrix(fit)
d <- as.data.frame(fit)