Example
Throughout the rest of the vignette we’ll use a hierarchical meta-analysis model described in section 5.5 of Gelman et al. (2003) as a running example. A hierarchical model is used to model the effect of coaching programs on college admissions tests. The data, shown in the table below, summarize the results of experiments conducted in eight high schools, with an estimated standard error for each. These data and model are of historical interest as an example of full Bayesian inference (Rubin 1981). For short, we call this the Eight Schools examples.
| School | Estimate (\(y_j\)) | Standard Error (\(\sigma_j\)) |
|---|---|---|
| A | 28 | 15 |
| B | 8 | 10 |
| C | -3 | 16 |
| D | 7 | 11 |
| E | -1 | 9 |
| F | 1 | 11 |
| G | 18 | 10 |
| H | 12 | 18 |
We use the Eight Schools example here because it is simple but also represents a nontrivial Markov chain simulation problem in that there is dependence between the parameters of original interest in the study — the effects of coaching in each of the eight schools — and the hyperparameter representing the variation of these effects in the modeled population. Certain implementations of a Gibbs sampler or a Hamiltonian Monte Carlo sampler can be slow to converge in this example.
The statistical model of interest is specified as
\[ \begin{aligned} y_j &\sim \mathsf{Normal}(\theta_j, \sigma_j), \quad j=1,\ldots,8 \\ \theta_j &\sim \mathsf{Normal}(\mu, \tau), \quad j=1,\ldots,8 \\ p(\mu, \tau) &\propto 1, \end{aligned} \]
where each \(\sigma_j\) is assumed known.
Write a Stan Program
RStan allows a Stan program to be coded in a text file (typically with suffix .stan) or in a R character vector (of length one). We put the following code for the Eight Schools model into the file schools.stan:
data {
int<lower=0> J; // number of schools
real y[J]; // estimated treatment effects
real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
real mu;
real<lower=0> tau;
vector[J] eta;
}
transformed parameters {
vector[J] theta;
theta = mu + tau * eta;
}
model {
target += normal_lpdf(eta | 0, 1);
target += normal_lpdf(y | theta, sigma);
}The first section of the Stan program above, the data block, specifies the data that is conditioned upon in Bayes Rule: the number of schools, \(J\), the vector of estimates, \((y_1, \ldots, y_J)\), and the vector of standard errors of the estimates \((\sigma_{1}, \ldots, \sigma_{J})\). Data are declared as integer or real and can be vectors (or, more generally, arrays) if dimensions are specified. Data can also be constrained; for example, in the above model \(J\) has been restricted to be at least \(1\) and the components of \(\sigma_y\) must all be positive.
The parameters block declares the parameters whose posterior distribution is sought. These are the the mean, \(\mu\), and standard deviation, \(\tau\), of the school effects, plus the standardized school-level effects \(\eta\). In this model, we let the unstandardized school-level effects, \(\theta\), be a transformed parameter constructed by scaling the standardized effects by \(\tau\) and shifting them by \(\mu\) rather than directly declaring \(\theta\) as a parameter. By parameterizing the model this way, the sampler runs more efficiently because the resulting multivariate geometry is more amendable to Hamiltonian Monte Carlo (Neal 2011).
Finally, the model block looks similar to standard statistical notation. (Just be careful: the second argument to Stan’s normal\((\cdot,\cdot)\) distribution is the standard deviation, not the variance as is usual in statistical notation). We have written the model in vector notation, which allows Stan to make use of more efficient algorithmic differentiation (AD). It would also be possible — but less efficient — to write the model by replacing normal_lpdf(y | theta,sigma) with a loop over the \(J\) schools,
for (j in 1:J)
target += normal_lpdf(y[j] | theta[j],sigma[j]);Stan has versions of many of the most useful R functions for statistical modeling, including probability distributions, matrix operations, and various special functions. However, the names of the Stan functions may differ from their R counterparts and, more subtly, the parameterizations of probability distributions in Stan may differ from those in R for the same distribution. To mitigate this problem, the lookup function can be passed an R function or character string naming an R function, and RStan will attempt to look up the corresponding Stan function, display its arguments, and give the page number in The Stan Development Team (2016) where the function is discussed.
StanFunction
415 normal_id_glm_lpdf
418 normal_log
419 normal_lpdf
553 std_normal_lpdf
Arguments
415 (real, matrix, real, vector, T);(vector, row_vector, vector, vector, vector)
418 (real, real, T);(vector, vector, vector)
419 (real, real, T);(vector, vector, vector)
553 (T);(vector)
ReturnType
415 T;real
418 T;real
419 T;real
553 T;real[1] "no matching Stan functions"If the lookup function fails to find an R function that corresponds to a Stan function, it will treat its argument as a regular expression and attempt to find matches with the names of Stan functions.
Preparing the Data
The stan function accepts data as a named list, a character vector of object names, or an environment. Alternatively, the data argument can be omitted and R will search for objects that have the same names as those declared in the data block of the Stan program. Here is the data for the Eight Schools example:
schools_data <- list(
J = 8,
y = c(28, 8, -3, 7, -1, 1, 18, 12),
sigma = c(15, 10, 16, 11, 9, 11, 10, 18)
)It would also be possible (indeed, encouraged) to read in the data from a file rather than to directly enter the numbers in the R script.
Sample from the Posterior Distribution
Next, we can call the stan function to draw posterior samples:
library(rstan)
fit1 <- stan(
file = "schools.stan", # Stan program
data = schools_data, # named list of data
chains = 4, # number of Markov chains
warmup = 1000, # number of warmup iterations per chain
iter = 2000, # total number of iterations per chain
cores = 1, # number of cores (could use one per chain)
refresh = 0 # no progress shown
)The stan function wraps the following three steps:
- Translate a model in Stan code to C++ code
- Compile the C++ code to a dynamic shared object (DSO) and load the DSO
- Sample given some user-specified data and other settings
A single call to stan performs all three steps, but they can also be executed one by one (see the help pages for stanc, stan_model, and sampling), which can be useful for debugging. In addition, Stan saves the DSO so that when the same model is fit again (possibly with new data and settings) we can avoid recompilation. If an error happens after the model is compiled but before sampling (e.g., problems with inputs like data and initial values), we can still reuse the compiled model.
The stan function returns a stanfit object, which is an S4 object of class "stanfit". For those who are not familiar with the concept of class and S4 class in R, refer to Chambers (2008). An S4 class consists of some attributes (data) to model an object and some methods to model the behavior of the object. From a user’s perspective, once a stanfit object is created, we are mainly concerned about what methods are defined.
If no error occurs, the returned stanfit object includes the sample drawn from the posterior distribution for the model parameters and other quantities defined in the model. If there is an error (e.g. a syntax error in the Stan program), stan will either quit or return a stanfit object that contains no posterior draws.
For class "stanfit", many methods such as print and plot are defined for working with the MCMC sample. For example, the following shows a summary of the parameters from the Eight Schools model using the print method:
Inference for Stan model: anon_model.
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 10% 50% 90% n_eff Rhat
theta[1] 11.38 0.15 8.40 2.16 10.36 22.58 3340 1
theta[2] 7.78 0.10 6.42 -0.25 7.76 15.56 4493 1
theta[3] 6.11 0.15 7.67 -3.35 6.58 14.77 2793 1
theta[4] 7.51 0.10 6.62 -0.50 7.46 15.40 4384 1
theta[5] 5.02 0.12 6.47 -3.31 5.54 12.68 3115 1
theta[6] 6.20 0.10 6.66 -2.30 6.50 14.00 4409 1
theta[7] 10.63 0.12 6.78 2.57 10.07 19.70 3426 1
theta[8] 8.17 0.15 7.85 -0.71 7.94 17.48 2829 1
mu 7.66 0.12 5.27 1.23 7.58 14.10 1811 1
tau 6.58 0.15 5.63 0.97 5.29 13.80 1420 1
lp__ -39.61 0.08 2.72 -43.20 -39.33 -36.36 1127 1
Samples were drawn using NUTS(diag_e) at Mon Mar 4 18:36:43 2024.
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).The last line of this output, lp__, is the logarithm of the (unnormalized) posterior density as calculated by Stan. This log density can be used in various ways for model evaluation and comparison (see, e.g., Vehtari and Ojanen (2012)).
Arguments to the stan Function
The primary arguments for sampling (in functions stan and sampling) include data, initial values, and the options of the sampler such as chains, iter, and warmup. In particular, warmup specifies the number of iterations that are used by the NUTS sampler for the adaptation phase before sampling begins. After the warmup, the sampler turns off adaptation and continues until a total of iter iterations (including warmup) have been completed. There is no theoretical guarantee that the draws obtained during warmup are from the posterior distribution, so the warmup draws should only be used for diagnosis and not inference. The summaries for the parameters shown by the print method are calculated using only post-warmup draws.
The optional init argument can be used to specify initial values for the Markov chains. There are several ways to specify initial values, and the details can be found in the documentation of the stan function. The vast majority of the time it is adequate to allow Stan to generate its own initial values randomly. However, sometimes it is better to specify the initial values for at least a subset of the objects declared in the parameters block of a Stan program.
Stan uses a random number generator (RNG) that supports parallelism. The initialization of the RNG is determined by the arguments seed and chain_id. Even if we are running multiple chains from one call to the stan function we only need to specify one seed, which is randomly generated by R if not specified.
Data Preprocessing and Passing
The data passed to stan will go through a preprocessing procedure. The details of this preprocessing are documented in the documentation for the stan function. Here we stress a few important steps. First, RStan allows the user to pass more objects as data than what is declared in the data block (silently omitting any unnecessary objects). In general, an element in the list of data passed to Stan from R should be numeric and its dimension should match the declaration in the data block of the model. So for example, the factor type in R is not supported as a data element for RStan and must be converted to integer codes via as.integer. The Stan modeling language distinguishes between integers and doubles (type int and real in Stan modeling language, respectively). The stan function will convert some R data (which is double-precision usually) to integers if possible.
The Stan language has scalars and other types that are sets of scalars, e.g. vectors, matrices, and arrays. As R does not have true scalars, RStan treats vectors of length one as scalars. However, consider a model with a data block defined as
data {
int<lower=1> N;
real y[N];
} in which N can be \(1\) as a special case. So if we know that N is always larger than \(1\), we can use a vector of length N in R as the data input for y (for example, a vector created by y <- rnorm(N)). If we want to prevent RStan from treating the input data for y as a scalar when \(N`\) is \(1\), we need to explicitly make it an array as the following R code shows:
y <- as.array(y)Stan cannot handle missing values in data automatically, so no element of the data can contain NA values. An important step in RStan’s data preprocessing is to check missing values and issue an error if any are found. There are, however, various ways of writing Stan programs that account for missing data (see The Stan Development Team (2016)).
Methods for the "stanfit" Class
The other vignette included with the rstan package discusses stanfit objects in greater detail and gives examples of accessing the most important content contained in the objects (e.g., posterior draws, diagnostic summaries). Also, a full list of available methods can be found in the documentation for the "stanfit" class at help("stanfit", "rstan"). Here we give only a few examples.
The plot method for stanfit objects provides various graphical overviews of the output. The default plot shows posterior uncertainty intervals (by default 80% (inner) and 95% (outer)) and the posterior median for all the parameters as well as lp__ (the log of posterior density function up to an additive constant):
'pars' not specified. Showing first 10 parameters by default.ci_level: 0.8 (80% intervals)outer_level: 0.95 (95% intervals)
The optional plotfun argument can be used to select among the various available plots. See help("plot,stanfit-method").
The traceplot method is used to plot the time series of the posterior draws. If we include the warmup draws by setting inc_warmup=TRUE, the background color of the warmup area is different from the post-warmup phase:
To assess the convergence of the Markov chains, in addition to visually inspecting traceplots we can calculate the split \(\hat{R}\) statistic. Split \(\hat{R}\) is an updated version of the \(\hat{R}\) statistic proposed in Gelman and Rubin (1992) that is based on splitting each chain into two halves. See the Stan manual for more details. The estimated \(\hat{R}\) for each parameter is included as one of the columns in the output from the summary and print methods.
Inference for Stan model: anon_model.
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.66 0.12 5.27 -2.83 4.34 7.58 10.97 18.10 1811 1
tau 6.58 0.15 5.63 0.20 2.45 5.29 9.14 20.59 1420 1
Samples were drawn using NUTS(diag_e) at Mon Mar 4 18:36:43 2024.
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).Again, see the additional vignette on stanfit objects for more details.
Sampling Difficulties
The best way to visualize the output of a model is through the ShinyStan interface, which can be accessed via the shinystan R package. ShinyStan facilitates both the visualization of parameter distributions and diagnosing problems with the sampler. The documentation for the shinystan package provides instructions for using the interface with stanfit objects.
In addition to using ShinyStan, it is also possible to diagnose some sampling problems using functions in the rstan package. The get_sampler_params function returns information on parameters related the performance of the sampler:
# all chains combined
sampler_params <- get_sampler_params(fit1, inc_warmup = TRUE)
summary(do.call(rbind, sampler_params), digits = 2) accept_stat__ stepsize__ treedepth__ n_leapfrog__ divergent__
Min. :0.00 Min. : 0.033 Min. :0.0 Min. : 1 Min. :0.000
1st Qu.:0.80 1st Qu.: 0.285 1st Qu.:3.0 1st Qu.: 7 1st Qu.:0.000
Median :0.95 Median : 0.343 Median :3.0 Median : 15 Median :0.000
Mean :0.84 Mean : 0.385 Mean :3.4 Mean : 13 Mean :0.009
3rd Qu.:0.99 3rd Qu.: 0.405 3rd Qu.:4.0 3rd Qu.: 15 3rd Qu.:0.000
Max. :1.00 Max. :10.289 Max. :7.0 Max. :127 Max. :1.000
energy__
Min. :35
1st Qu.:42
Median :44
Mean :45
3rd Qu.:47
Max. :61 [[1]]
accept_stat__ stepsize__ treedepth__ n_leapfrog__ divergent__
Min. :0.00 Min. :0.063 Min. :0.0 Min. : 1 Min. :0.0000
1st Qu.:0.75 1st Qu.:0.343 1st Qu.:3.0 1st Qu.: 7 1st Qu.:0.0000
Median :0.94 Median :0.343 Median :3.0 Median :15 Median :0.0000
Mean :0.82 Mean :0.395 Mean :3.4 Mean :12 Mean :0.0085
3rd Qu.:0.99 3rd Qu.:0.383 3rd Qu.:4.0 3rd Qu.:15 3rd Qu.:0.0000
Max. :1.00 Max. :8.548 Max. :6.0 Max. :95 Max. :1.0000
energy__
Min. :36
1st Qu.:42
Median :44
Mean :45
3rd Qu.:47
Max. :61
[[2]]
accept_stat__ stepsize__ treedepth__ n_leapfrog__ divergent__
Min. :0.00 Min. :0.033 Min. :0.0 Min. : 1 Min. :0.000
1st Qu.:0.83 1st Qu.:0.278 1st Qu.:3.0 1st Qu.: 7 1st Qu.:0.000
Median :0.96 Median :0.278 Median :4.0 Median : 15 Median :0.000
Mean :0.85 Mean :0.364 Mean :3.5 Mean : 13 Mean :0.011
3rd Qu.:0.99 3rd Qu.:0.410 3rd Qu.:4.0 3rd Qu.: 15 3rd Qu.:0.000
Max. :1.00 Max. :5.402 Max. :7.0 Max. :127 Max. :1.000
energy__
Min. :35
1st Qu.:42
Median :44
Mean :44
3rd Qu.:46
Max. :61
[[3]]
accept_stat__ stepsize__ treedepth__ n_leapfrog__ divergent__
Min. :0.00 Min. : 0.045 Min. :0.0 Min. : 1 Min. :0.000
1st Qu.:0.82 1st Qu.: 0.356 1st Qu.:3.0 1st Qu.: 7 1st Qu.:0.000
Median :0.95 Median : 0.356 Median :3.0 Median :15 Median :0.000
Mean :0.84 Mean : 0.420 Mean :3.3 Mean :12 Mean :0.009
3rd Qu.:0.99 3rd Qu.: 0.424 3rd Qu.:4.0 3rd Qu.:15 3rd Qu.:0.000
Max. :1.00 Max. :10.289 Max. :6.0 Max. :63 Max. :1.000
energy__
Min. :36
1st Qu.:42
Median :44
Mean :45
3rd Qu.:47
Max. :57
[[4]]
accept_stat__ stepsize__ treedepth__ n_leapfrog__ divergent__
Min. :0.00 Min. :0.047 Min. :0.0 Min. : 1 Min. :0.0000
1st Qu.:0.82 1st Qu.:0.294 1st Qu.:3.0 1st Qu.: 7 1st Qu.:0.0000
Median :0.96 Median :0.294 Median :4.0 Median : 15 Median :0.0000
Mean :0.84 Mean :0.362 Mean :3.5 Mean : 13 Mean :0.0075
3rd Qu.:0.99 3rd Qu.:0.399 3rd Qu.:4.0 3rd Qu.: 15 3rd Qu.:0.0000
Max. :1.00 Max. :5.423 Max. :6.0 Max. :127 Max. :1.0000
energy__
Min. :35
1st Qu.:42
Median :44
Mean :45
3rd Qu.:47
Max. :59 Here we see that there are a small number of divergent transitions, which are identified by divergent__ being \(1\). Ideally, there should be no divergent transitions after the warmup phase. The best way to try to eliminate divergent transitions is by increasing the target acceptance probability, which by default is \(0.8\). In this case the mean of accept_stat__ is close to \(0.8\) for all chains, but has a very skewed distribution because the median is near \(0.95\). We could go back and call stan again and specify the optional argument control=list(adapt_delta=0.9) to try to eliminate the divergent transitions. However, sometimes when the target acceptance rate is high, the stepsize is very small and the sampler hits its limit on the number of leapfrog steps it can take per iteration. In this case, it is a non-issue because each chain has a treedepth__ of at most \(7\) and the default is \(10\). But if any treedepth__ were \(11\), then it would be wise to increase the limit by passing control=list(max_treedepth=12) (for example) to stan. See the vignette on stanfit objects for more on the structure of the object returned by get_sampler_params.
We can also make a graphical representation of (much of the) the same information using pairs. The “pairs”" plot can be used to get a sense of whether any sampling difficulties are occurring in the tails or near the mode:
Warning in par(usr): argument 1 does not name a graphical parameter
Warning in par(usr): argument 1 does not name a graphical parameter
Warning in par(usr): argument 1 does not name a graphical parameter
In the plot above, the marginal distribution of each selected parameter is included as a histogram along the diagonal. By default, draws with below-median accept_stat__ (MCMC proposal acceptance rate) are plotted below the diagonal and those with above-median accept_stat__ are plotted above the diagonal (this can be changed using the condition argument). Each off-diagonal square represents a bivariate distribution of the draws for the intersection of the row-variable and the column-variable. Ideally, the below-diagonal intersection and the above-diagonal intersection of the same two variables should have distributions that are mirror images of each other. Any yellow points would indicate transitions where the maximum treedepth__ was hit, and red points indicate a divergent transition.