In this vignette we present RStan, the R interface to Stan. Stan is a C++ library for Bayesian inference using the No-U-Turn sampler (a variant of Hamiltonian Monte Carlo) or frequentist inference via optimization. We illustrate the features of RStan through an example in Gelman et al. (2003).
Stan is a C++ library for Bayesian modeling and inference that primarily uses the No-U-Turn sampler (NUTS) (Hoffman and Gelman 2012) to obtain posterior simulations given a user-specified model and data. Alternatively, Stan can utilize the LBFGS optimization algorithm to maximize an objective function, such as a log-likelihood. The R package rstan provides RStan, the R interface to Stan. The rstan package allows one to conveniently fit Stan models from R (R Core Team 2014) and access the output, including posterior inferences and intermediate quantities such as evaluations of the log posterior density and its gradients.
In this vignette we provide a concise introduction to the functionality included in the rstan package. Stan’s website mc-stan.org has additional details and provides up-to-date information about how to operate both Stan and its many interfaces including RStan. See, for example, RStan Getting Started (The Stan Development Team 2014).
Stan has a modeling language, which is similar to but not identical to that of the Bayesian graphical modeling package BUGS (Lunn et al. 2000). A parser translates a model expressed in the Stan language to C++ code, whereupon it is compiled to an executable program and loaded as a Dynamic Shared Object (DSO) in R which can then be called by the user.
A C++ compiler, such as g++ or clang++, is required for
this process. For instructions on installing a C++ compiler for use with
RStan see RStan-Getting-Started.
The rstan package also depends heavily on several other R packages:
These dependencies should be automatically installed if you install the rstan package via one of the conventional mechanisms.
The following is a typical workflow for using Stan via RStan for Bayesian inference.
.stan extension, although it
can also be done using a character string within R.stanc
function.Conveniently, steps 2, 3, and 4, above, are all performed implicitly
by a single call to the stan function.
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 (yj) | Standard Error (σ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 σj is assumed known.
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, (y1, …, yJ),
and the vector of standard errors of the estimates (σ1, …, σ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 σy must all be
positive.
The parameters block declares the parameters whose
posterior distribution is sought. These are the the mean, μ, and standard deviation, τ, of the school effects, plus the
standardized school-level effects η. In this model, we let the
unstandardized school-level effects, θ, be a transformed parameter
constructed by scaling the standardized effects by τ and shifting them by μ rather than directly declaring
θ 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(⋅, ⋅) 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.
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.
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
)Error in stanc(file = file, model_code = model_code, model_name = model_name, : 0
Syntax error in 'string', line 3, column 9 to column 10, parsing error:
-------------------------------------------------
1: data {
2: int<lower=0> J; // number of schools
3: real y[J]; // estimated treatment effects
^
4: real<lower=0> sigma[J]; // s.e. of effect estimates
5: }
-------------------------------------------------
";" expected after variable declaration.
It looks like you are trying to use the old array syntax.
Please use the new syntax:
array[J] real y;The stan function wraps the following three steps:
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:
Error: object 'fit1' not foundThe 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)).
stan FunctionThe 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.
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)).
"stanfit" ClassThe 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):
Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'fit1' not foundThe 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:
Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'traceplot': object 'fit1' not foundTo assess the convergence of the Markov chains, in addition to
visually inspecting traceplots we can calculate the split R̂ statistic. Split R̂ is an updated version of the 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 R̂ for
each parameter is included as one of the columns in the output from the
summary and print methods.
Error: object 'fit1' not foundAgain, see the additional vignette on stanfit objects for more details.
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:
Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_sampler_params': object 'fit1' not foundError in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': second argument must be a list