Working with posteriors

library(cmdstanr)
library(posterior)
library(ggplot2)
theme_set(bayesplot::theme_default())

Summary statistics

We can easily customize the summary statistics reported by $summary() and $print().

fit <- cmdstanr_example("schools_ncp", method = "sample")
fit$summary()

By default, all variables are summarized with the following functions:

posterior::default_summary_measures()

To change the variables summarized, use the variables argument:

fit$summary(variables = c("mu", "tau", "theta"))

We can also change which functions are used:

fit$summary(variables = c("mu", "tau"), mean, sd)

To summarize all variables with non-default functions, it is necessary to explicitly set the variables argument, either to NULL or the full vector of variable names.

fit$summary(variables = NULL, "mean", "median")

Summary functions can be specified by character string, function, or using a formula (or anything else supported by rlang::as_function()). If these arguments are named, those names will be used in the tibble output. If the summary results are named they will take precedence.

my_sd <- function(x) c(My_SD = sd(x))
fit$summary(
  c("mu", "tau"), 
  MEAN = mean, 
  "median",
  my_sd,
  ~quantile(.x, probs = c(0.1, 0.9)),
  Minimum = function(x) min(x)
)

Arguments to all summary functions can also be specified with .args.

fit$summary(c("mu", "tau"), quantile, .args = list(probs = c(0.025, .05, .95, .975)))

Each summary function is applied separately to each variable and receives a matrix whose rows are saved iterations and whose columns are chains.

fit$summary(variables = "theta", dim, colMeans)

For this reason users may have unexpected results if they use stats::var() directly, as it will return a covariance matrix. An alternative is the distributional::variance() function, which can also be accessed via posterior::variance().

fit$summary(c("mu", "tau"), posterior::variance, ~var(as.vector(.x)))

Summary functions need not return numeric values when used with $summary(). The $print() method requires numeric summary columns because it rounds them to the requested number of digits.

strict_pos <- function(x) if (all(x > 0)) "yes" else "no"
fit$summary(variables = c("mu", "tau", "theta"), "Strictly Positive" = strict_pos)
# fit$print(variables = NULL, "Strictly Positive" = strict_pos)

For more information, see posterior::summarise_draws(), which is called internally by $summary().

Extracting posterior draws/samples

The $draws() method extracts draws in formats provided by the posterior package. The Getting started with CmdStanR vignette introduces the most commonly used formats and how to convert between them.

fit$draws("mu")
fit$draws("theta")
fit$draws(c("mu", "theta[1]"), format = "df")

For MCMC fits, inc_warmup = TRUE includes warmup draws, but only if save_warmup = TRUE was specified when fitting the model.

For more ways to manipulate draws, see the posterior package vignettes and documentation.

Structured draws similar to rstan::extract()

The posterior package provides two useful ways to work with variables while preserving their original dimensions.

posterior::extract_list_of_variable_arrays() returns a named list containing one array per variable. Setting with_chains = FALSE combines the chains, giving the same general structure as the list returned by rstan::extract():

draw_arrays <- extract_list_of_variable_arrays(
  fit$draws(),
  variables = c("mu", "theta"),
  with_chains = FALSE
)
str(draw_arrays)
dim(draw_arrays$theta)

The first dimension of each array indexes draws, and any remaining dimensions match the dimensions of the corresponding Stan variable.

Alternatively, the posterior package’s rvar format represents each variable as a multidimensional random variable, with its posterior draws handled behind the scenes:

draws_rvars <- as_draws_rvars(
  fit$draws(c("mu", "theta"))
)
theta_rvar <- draws_rvars$theta

# Compute the difference for every draw using natural vector indexing
# The posterior draws are handled automatically
theta_difference <- theta_rvar[1] - theta_rvar[2]
theta_difference

hist(
  draws_of(theta_difference),
  main = "Difference between theta[1] and theta[2]",
  xlab = "theta[1] - theta[2]"
)

# Direct access to the underlying draws is also available with posterior::draws_of
theta_array <- draws_of(theta_rvar)
dim(theta_array)

The object theta_rvar behaves like the vector declared in the Stan program. theta_array provides direct access to its underlying draws, with the first dimension indexing draws. See the rvar vignette for details.

Plotting the draws of a vector

Because theta_array has draws in the first dimension and the vector index (the eight schools) in the second, we can reshape it into a long data frame and overlay the individual draws.

theta_plot <- draw_arrays$theta

theta_df <- data.frame(
  .draw = rep(seq_len(nrow(theta_plot)), times = ncol(theta_plot)),
  school = rep(seq_len(ncol(theta_plot)), each = nrow(theta_plot)),
  theta = c(theta_plot)
)

ggplot(theta_df, aes(school, theta, group = .draw)) +
  geom_line(alpha = 0.01)

The reshaping above uses only base R. Tidyverse users can produce the same plot directly from the draws data frame (format = "df") with tidyr::pivot_longer(), extracting the vector index from variable names like theta[1]:

fit$draws("theta", format = "df") |>
  tidyr::pivot_longer(
    cols = dplyr::starts_with("theta"),
    names_to = "school",
    names_transform = readr::parse_number,
    values_to = "theta"
  ) |>
  ggplot(aes(school, theta, group = .draw)) +
  geom_line(alpha = 0.01)

Here school is simply the index into the theta vector. In many models the vector index corresponds to a meaningful covariate, for example the time points of a time series. In that case you can replace school with the associated covariate values to plot each draw as a function of that covariate.