Getting started with CmdStanR

Introduction

CmdStanR (Command Stan R) is a lightweight interface to Stan for R users that provides an alternative to the traditional RStan interface. See the Comparison with RStan section later in this vignette for more details on how the two interfaces differ.

Using CmdStanR requires installing the cmdstanr R package and also CmdStan, the command line interface to Stan. First we install the R package by running the following command in R.

# we recommend running this in a fresh R session or restarting your current session
install.packages("cmdstanr", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))

We can now load the package like any other R package. We’ll also load the bayesplot and posterior packages to use later in examples.

library(cmdstanr)
library(posterior)
library(bayesplot)
color_scheme_set("brightblue")

Installing CmdStan

CmdStanR requires a working installation of version 2.35 or later of CmdStan, the command line interface to Stan. If you don’t have CmdStan installed, CmdStanR can install it for you, assuming you have a suitable C++ toolchain. The requirements are described in the CmdStan Guide:

To double check that your toolchain is set up properly you can call the check_cmdstan_toolchain() function:

check_cmdstan_toolchain()

If your toolchain is configured correctly then CmdStan can be installed by calling the install_cmdstan() function:

install_cmdstan(cores = 2)

Before CmdStanR can be used it needs to know where the CmdStan installation is located. When the package is loaded it tries to help automate this to avoid having to manually set the path every session:

  1. If the environment variable "CMDSTAN" points directly to a valid CmdStan installation at load time, that path is used for the R session. If it instead points to an existing parent directory containing versioned CmdStan installations, the installation with the largest version number is used. This is useful if your CmdStan installation is not located in the default directory that would have been used by install_cmdstan() (see #2).

  2. If no environment variable is found when loaded but any directory in the form ".cmdstan/cmdstan-[version]", for example ".cmdstan/cmdstan-2.35.0", exists in the user’s home directory (Sys.getenv("HOME"), not the current working directory) then the path to the CmdStan with the largest version number will be set as the path to CmdStan for the R session. This is the same as the default directory that install_cmdstan() uses to install the latest version of CmdStan, so if that’s how you installed CmdStan you shouldn’t need to manually set the path to CmdStan when loading CmdStanR.

If neither of these applies (or you want to subsequently change the path) you can use the set_cmdstan_path() function:

set_cmdstan_path(PATH_TO_CMDSTAN)

To check the path to the CmdStan installation and the CmdStan version number you can use cmdstan_path() and cmdstan_version():

cmdstan_path()
cmdstan_version()

Compiling a model

The cmdstan_model() function creates a new CmdStanModel object from a file containing a Stan program. Under the hood, CmdStan is called to translate a Stan program to C++ and create a compiled executable. Here we’ll use the example Stan program that comes with the CmdStan installation:

file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan")
mod <- cmdstan_model(file)

The object mod is an R6 reference object of class CmdStanModel and behaves similarly to R’s reference class objects and those in object oriented programming languages. Methods are accessed using the $ operator. This design choice allows for CmdStanR and CmdStanPy to provide a similar user experience and share many implementation details.

The Stan program can be printed using the $print() method:

mod$print()

The path to the compiled executable is returned by the $exe_file() method:

mod$exe_file()

Running MCMC

The $sample() method for CmdStanModel objects runs Stan’s default MCMC algorithm. The data argument accepts a named list of R objects (like for RStan) or a path to a data file compatible with CmdStan (JSON or R dump).

# names correspond to the data block in the Stan program
data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))

fit <- mod$sample(
  data = data_list,
  seed = 123,
  chains = 4,
  parallel_chains = 4,
  refresh = 500 # print update every 500 iters
)

There are many more arguments that can be passed to the $sample() method. For example, initial values can be specified with the init argument, including by passing a fitted model object produced by any of CmdStanR’s fitting methods. For details follow this link to its separate documentation page:

The $sample() method creates R6 CmdStanMCMC objects, which have many associated methods. Below we will demonstrate some of the most important methods. For a full list, follow this link to the CmdStanMCMC documentation:

Posterior summary statistics

Summaries from the posterior package

The $summary() method calls summarise_draws() from the posterior package. The first argument specifies the variables to summarize and any arguments after that are passed on to posterior::summarise_draws() to specify which summaries to compute, whether to use multiple cores, etc.

fit$summary()
fit$summary(variables = c("theta", "lp__"), "mean", "sd")

# use a formula to summarize arbitrary functions, e.g. Pr(theta <= 0.5)
fit$summary("theta", pr_lt_half = ~ mean(. <= 0.5))

# summarise all variables with default and additional summary measures
fit$summary(
  variables = NULL,
  posterior::default_summary_measures(),
  extra_quantiles = ~posterior::quantile2(., probs = c(0.025, 0.975))
)

CmdStan’s stansummary utility

CmdStan itself provides a stansummary utility that can be called using the $cmdstan_summary() method. This method prints summaries produced by CmdStan but, unlike $summary(), does not return them as a data frame for further use in R.

Posterior draws

Extracting draws

The $draws() method can be used to extract the posterior draws in formats provided by the posterior package. Here we demonstrate only the draws_array and draws_df formats, but the posterior package supports other useful formats as well.

# default is a 3-D draws_array object from the posterior package
# iterations x chains x variables
draws_arr <- fit$draws() # or format="array"
str(draws_arr)

# draws x variables data frame
draws_df <- fit$draws(format = "df")
str(draws_df)
print(draws_df)

To convert an existing draws object to a different format use the posterior::as_draws_*() functions.

# this should be identical to draws_df created via draws(format = "df")
draws_df_2 <- as_draws_df(draws_arr)
identical(draws_df, draws_df_2)

In general, converting to a different draws format in this way will be slower than just setting the appropriate format initially in the call to the $draws() method, but in most cases the speed difference will be minor.

The vignette Working with Posteriors has more details on posterior draws, including how to reproduce the structured output RStan users are accustomed to getting from rstan::extract().

Plotting draws

Plotting posterior distributions is as easy as passing the object returned by the $draws() method directly to plotting functions in our bayesplot package.

mcmc_hist(fit$draws("theta"))

Sampler diagnostics

Extracting diagnostic values for each iteration and chain

The $sampler_diagnostics() method extracts the values of the sampler parameters (treedepth__, divergent__, etc.) in formats supported by the posterior package. The default is as a 3-D array (iteration x chain x variable).

# this is a draws_array object from the posterior package
str(fit$sampler_diagnostics())

# this is a draws_df object from the posterior package
str(fit$sampler_diagnostics(format = "df"))

Sampler diagnostic warnings and summaries

The $diagnostic_summary() method will display any sampler diagnostic warnings and return a summary of diagnostics for each chain.

fit$diagnostic_summary()

We see the number of divergences for each of the four chains, the number of times the maximum treedepth was hit for each chain, and the E-BFMI for each chain.

In this case there were no warnings, so in order to demonstrate the warning messages we’ll use one of the CmdStanR example models that suffers from divergences.

fit_with_warning <- cmdstanr_example("schools")

After fitting there is a warning about divergences. We can also regenerate this warning message later using fit_with_warning$diagnostic_summary().

diagnostics <- fit_with_warning$diagnostic_summary()
print(diagnostics)

# number of divergences reported in warning is the sum of the per chain values
sum(diagnostics$num_divergent)

CmdStan’s diagnose utility

CmdStan itself provides a diagnose utility that can be called using the $cmdstan_diagnose() method. This method prints CmdStan’s diagnostic output but does not return it as a structured R object.

Running optimization and variational inference

CmdStanR also supports running Stan’s optimization algorithms and its algorithms for variational approximation of full Bayesian inference. These are run via the $optimize(), $laplace(), $variational(), and $pathfinder() methods, which are called in a similar way to the $sample() method demonstrated above.

Optimization

Following CmdStan’s terminology, we can find the (penalized) maximum likelihood estimate (MLE) using $optimize().

fit_mle <- mod$optimize(data = data_list, seed = 123)
fit_mle$print() # includes lp__ (log prob calculated by Stan program)
fit_mle$mle("theta")

Here’s a plot comparing the penalized MLE to the posterior distribution of theta.

mcmc_hist(fit$draws("theta")) +
  vline_at(fit_mle$mle("theta"), size = 1.5)

For optimization, by default the mode is calculated without the Jacobian adjustment for constrained variables. Including the Jacobian adjustment by setting jacobian=TRUE produces a maximum a posteriori (MAP) estimate. More precisely, the default finds a mode of the target in the constrained parameter space, whereas jacobian=TRUE finds a mode of the corresponding density in the unconstrained space. See the Maximum Likelihood Estimation section of the CmdStan User’s Guide for more details.

fit_map <- mod$optimize(
  data = data_list,
  jacobian = TRUE,
  seed = 123
)

Laplace Approximation

The $laplace() method produces a sample from a normal approximation centered at the mode of a distribution in the unconstrained space. Following CmdStan’s terminology, if the mode is a MAP estimate, the draws provide an estimate of the mean and standard deviation of the posterior distribution. If the mode is an MLE, the draws provide an estimate of the standard error of the likelihood. Whether the mode is called MAP or MLE depends on the value of the jacobian argument when running optimization. The Laplace jacobian setting must match the setting used for optimization so that both use the same target density. See the Laplace Sampling chapter of the CmdStan User’s Guide for more details.

Here we pass in the fit_map object from above as the mode argument. If mode is omitted then optimization will be run internally before taking draws from the normal approximation.

fit_laplace <- mod$laplace(
    mode = fit_map,
    draws = 4000,
    data = data_list,
    seed = 123,
    refresh = 1000
  )
fit_laplace$print("theta")
mcmc_hist(fit_laplace$draws("theta"), binwidth = 0.025)

Variational (ADVI)

We can run Stan’s experimental Automatic Differentiation Variational Inference (ADVI) using the $variational() method. For details on the ADVI algorithm see the CmdStan User’s Guide.

fit_vb <- mod$variational(
  data = data_list,
  seed = 123,
  draws = 4000
)
fit_vb$print("theta")
mcmc_hist(fit_vb$draws("theta"), binwidth = 0.025)

Variational (Pathfinder)

Stan version 2.33 introduced the Pathfinder variational inference method. For details on how Pathfinder works see the section in the CmdStan User’s Guide. Pathfinder is run using the $pathfinder() method.

fit_pf <- mod$pathfinder(
  data = data_list,
  seed = 123,
  draws = 4000
)
fit_pf$print("theta")

Let’s extract the draws, make the same plot we made after running the other algorithms, and compare them all. In this simple example the distributions are quite similar, but this will not always be the case for more challenging problems.

mcmc_hist(fit_pf$draws("theta"), binwidth = 0.025) +
  ggplot2::labs(subtitle = "Approximate posterior from pathfinder") +
  ggplot2::xlim(0, 1)
mcmc_hist(fit_vb$draws("theta"), binwidth = 0.025) +
  ggplot2::labs(subtitle = "Approximate posterior from variational") +
  ggplot2::xlim(0, 1)
mcmc_hist(fit_laplace$draws("theta"), binwidth = 0.025) +
  ggplot2::labs(subtitle = "Approximate posterior from Laplace") +
  ggplot2::xlim(0, 1)
mcmc_hist(fit$draws("theta"), binwidth = 0.025) +
  ggplot2::labs(subtitle = "Posterior from MCMC") +
  ggplot2::xlim(0, 1)

For more details on the $optimize(), $laplace(), $variational(), and $pathfinder() methods, follow these links to their documentation pages.

Running standalone generated quantities

The $generate_quantities() method runs a Stan program’s generated quantities block using parameter values from a previous fit, without rerunning inference. It accepts fitted model objects produced by any of the fitting methods described above. See the method documentation for details and a complete example.

Saving fitted model objects

The $save_object() method provided by CmdStanR is the most convenient way to save a fitted model object to disk and ensure that all of the contents are available when reading the object back into R. By default, fit$save_object() will use the RDS format to save the object. The saved object can then be read back into R using readRDS().

fit$save_object(file = "fit.RDS")

fit2 <- readRDS("fit.RDS")

But if your model object is large, then fit$save_object() can take a long time if saving in the default RDS format. For large objects, we recommend using the much faster qs2 format. The saved object can then be read back into R using qs2::qs_read().

fit$save_object(file = "fit.qs2", format = "qs2")

fit2 <- qs2::qs_read("fit.qs2")

Comparison with RStan

Additional resources

There are additional vignettes available that discuss many other aspects of using CmdStanR. These can be found online at the CmdStanR website:

To ask a question please post on the Stan forums:

To report a bug, suggest a feature (including additions to these vignettes), or to start contributing to CmdStanR development (new contributors welcome!) please open an issue on GitHub: