Using tidy data with Bayesian models

Philosophy

There are a few core ideas that run through the tidybayes API that should (hopefully) make it easy to use:

1. Tidy data does not always mean all parameter names as values. In contrast to the ggmcmc library (which translates model results into a data frame with a Parameter and value column), the spread_draws function in tidybayes produces data frames where the columns are named after parameters and (in some cases) indices of those parameters, as automatically as possible and using a syntax as close to the same way you would refer to those variables in the model’s language as possible. A similar function to ggmcmc’s approach is also provided in gather_draws, since sometimes you do want variable names as values in a column. The goal is for tidybayes to do the tedious work of figuring out how to make a data frame look the way you need it to, including turning parameters with indices like "b[1,2]" and the like into tidy data for you.

2. Fit into the tidyverse. tidybayes methods fit into a workflow familiar to users of the tidyverse (dplyr, tidyr, ggplot2, etc), which means fitting into the pipe (%>%) workflow, using and respecting grouped data frames (thus spread_draws and gather_draws return results already grouped by variable indices, and methods like median_qi calculate point summaries and intervals for variables and groups simultaneously), and not reinventing too much of the wheel if it is already made easy by functions provided by existing tidyverse packages (unless it makes for much clearer code for a common idiom). For compatibility with other package column names (such as broom::tidy), tidybayes provides transformation functions like to_broom_names that can be dropped directly into data transformation pipelines.

3. Focus on composable operations and plotting primitives, not monolithic plots and operations. Several other packages (notably bayesplot and ggmcmc) already provide an excellent variety of pre-made methods for plotting Bayesian results. tidybayes shies away from duplicating this functionality. Instead, it focuses on providing composable operations for generating and manipulating Bayesian samples in a tidy data format, and graphical primitives for ggplot that allow you to build custom plots easily. Most simply, where bayesplot and ggmcmc tend to have functions with many options that return a full ggplot object, tidybayes tends towards providing primitives (like geoms) that you can compose and combine into your own custom plots. I believe both approaches have their place: pre-made functions are especially useful for common, quick operations that don’t need customization (like many diagnostic plots), while composable operations tend to be useful for more complex custom plots (in my opinion).

4. Sensible defaults make life easy. But options (and the data being tidy in the first place) make it easy to go your own way when you need to.

5. Variable names in models should be descriptive, not cryptic. This principle implies avoiding cryptic (and short) subscripts in favor of longer (but descriptive) ones. This is a matter of readability and accessibility of models to others. For example, a common pattern among Stan users (and in the Stan manual) is to use variables like J to refer to the number of elements in a group (e.g., number of participants) and a corresponding index like j to refer to specific elements in that group. I believe this sacrifices too much readability for the sake of concision; I prefer a pattern like n_participant for the size of the group and participant (or a mnemonic short form like p) for specific elements. In functions where names are auto-generated (like compose_data), tidybayes will (by default) assume you want these sorts of more descriptive names; however, you can always override the default naming scheme.

Supported model types

tidybayes aims to support a variety of models with a uniform interface. Currently supported models include rstan, cmdstanr, brms, rstanarm, runjags, rjags, jagsUI, coda::mcmc and coda::mcmc.list, posterior::draws, MCMCglmm, and anything with its own as.mcmc.list implementation. If you install the tidybayes.rethinking package, models from the rethinking package are also supported.

For an up-to-date list of supported models, see ?"tidybayes-models".

Extracting model variable indices into a separate column in a tidy format data frame

Now that we have our results, the fun begins: getting the draws out in a tidy format! The default methods in Stan for extracting draws from the model do so in a nested format:

str(rstan::extract(m))

## List of 6
##  $ overall_mean     : num [1:4000(1d)] 0.91 0.925 1.485 1.055 0.762 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ condition_zoffset: num [1:4000, 1:5] -0.353 -0.401 -1.11 -0.862 -0.314 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ iterations: NULL
##   .. ..$           : NULL
##  $ response_sd      : num [1:4000(1d)] 0.566 0.555 0.504 0.559 0.555 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ condition_mean_sd: num [1:4000(1d)] 1.48 1.42 1.17 1.15 1.43 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ condition_mean   : num [1:4000, 1:5] 0.3862 0.3568 0.1828 0.0619 0.3131 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ iterations: NULL
##   .. ..$           : NULL
##  $ lp__             : num [1:4000(1d)] 2.39 3.24 1.92 -0.11 4.47 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL

There are also methods for extracting draws as matrices or data frames in Stan (and other model types, such as JAGS and MCMCglmm, have their own formats).

The spread_draws method yields a common format for all model types supported by tidybayes. It lets us instead extract draws into a data frame in tidy format, with a .chain and .iteration column storing the chain and iteration for each row (if available), a .draw column that uniquely indexes each draw, and the remaining columns corresponding to model variables or variable indices. The spread_draws method accepts any number of column specifications, which can include names for variables and names for variable indices. For example, we can extract the condition_mean variable as a tidy data frame, and put the value of its first (and only) index into the condition column, using a syntax that directly echoes how we would specify indices of the condition_mean variable in the model itself:

m %>%
  spread_draws(condition_mean[condition]) %>%
  head(10)

Automatically converting columns and indices back into their original data types

As-is, the resulting variables don’t know anything about where their indices came from. The index of the condition_mean variable was originally derived from the condition factor in the ABC data frame. But Stan doesn’t know this: it is just a numeric index to Stan, so the condition column just contains numbers (1, 2, 3, 4, 5) instead of the factor levels these numbers correspond to ("A", "B", "C", "D", "E").

We can recover this missing type information by passing the model through recover_types before using spread_draws. In itself recover_types just returns a copy of the model, with some additional attributes that store the type information from the data frame (or other objects) that you pass to it. This doesn’t have any useful effect by itself, but functions like spread_draws use this information to convert any column or index back into the data type of the column with the same name in the original data frame. In this example, spread_draws recognizes that the condition column was a factor with five levels ("A", "B", "C", "D", "E") in the original data frame, and automatically converts it back into a factor:

m %>%
  recover_types(ABC) %>%
  spread_draws(condition_mean[condition]) %>%
  head(10)

Because we often want to make multiple separate calls to spread_draws, it is often convenient to decorate the original model using recover_types immediately after it has been fit, so we only have to call it once:

m %<>% recover_types(ABC)

Now we can omit the recover_types call before subsequent calls to spread_draws.

With simple variables, wide format

tidybayes provides a family of functions for generating point summaries and intervals from draws in a tidy format. These functions follow the naming scheme [median|mean|mode]_[qi|hdi], for example, median_qi, mean_qi, mode_hdi, and so on. The first name (before the _) indicates the type of point summary, and the second name indicates the type of interval. qi yields a quantile interval (a.k.a. equi-tailed interval, central interval, or percentile interval) and hdi yields a highest density interval. Custom point or interval functions can also be applied using the point_interval function.

For example, we might extract the draws corresponding to the overall mean and standard deviation of observations:

m %>%
  spread_draws(overall_mean, response_sd) %>%
  head(10)

Like with condition_mean[condition], this gives us a tidy data frame. If we want the median and 95% quantile interval of the variables, we can apply median_qi:

m %>%
  spread_draws(overall_mean, response_sd) %>%
  median_qi(overall_mean, response_sd)

median_qi summarizes each input column using its median. If there are multiple columns to summarize, each gets its own x.lower and x.upper column (for each column x) corresponding to the bounds of the .width% interval. If there is only one column, the names .lower and .upper are used for the interval bounds.

We can specify the columns we want to get medians and intervals from, as above, or if we omit the list of columns, median_qi will use every column that is not a grouping column or a special column (like .chain, .iteration, or .draw). Thus in the above example, overall_mean and response_sd are redundant arguments to median_qi because they are also the only columns we gathered from the model. So we can simplify the previous code to the following:

m %>%
  spread_draws(overall_mean, response_sd) %>%
  median_qi()

With indexed variables

When we have a variable with one or more indices, such as condition_mean, we can apply median_qi (or other functions in the point_interval family) as we did before:

m %>%
  spread_draws(condition_mean[condition]) %>%
  median_qi()

How did median_qi know what to aggregate? Data frames returned by spread_draws are automatically grouped by all index variables you pass to it; in this case, that means it groups by condition. median_qi respects groups, and calculates the point summaries and intervals within all groups. Then, because no columns were passed to median_qi, it acts on the only non-special (.-prefixed) and non-group column, condition_mean. So the above shortened syntax is equivalent to this more verbose call:

m %>%
  spread_draws(condition_mean[condition]) %>%
  group_by(condition) %>%    # this line not necessary (done automatically by spread_draws)
  median_qi(condition_mean)

When given only a single column, median_qi will use the names .lower and .upper for the lower and upper ends of the intervals.

tidybayes also provides an implementation of posterior::summarise_draws() for grouped data frames (tidybayes::summaries_draws.grouped_df()), which you can use to quickly get convergence diagnostics:

m %>%
  spread_draws(condition_mean[condition]) %>%
  summarise_draws()

Using geom_pointinterval

Plotting medians and intervals is straightforward using ggdist::geom_pointinterval() or ggdist::stat_pointinterval(), which are similar to ggplot2::geom_pointrange() but with sensible defaults for multiple intervals. For example:

m %>%
  spread_draws(condition_mean[condition]) %>%
  ggplot(aes(y = fct_rev(condition), x = condition_mean)) +
  stat_pointinterval()

These functions have .width = c(.66, .95) by default (showing 66% and 95% intervals), but this can be changed by passing a .width argument to ggdist::stat_pointinterval().

Intervals with posterior violins (“eye plots”): stat_eye

The ggdist::stat_halfeye() geom provides a shortcut to generating “half-eye plots” (combinations of intervals and densities). This example also demonstrates how to change the interval probability (here, to 90% and 50% intervals):

m %>%
  spread_draws(condition_mean[condition]) %>%
  ggplot(aes(y = fct_rev(condition), x = condition_mean)) +
  stat_halfeye(.width = c(.90, .5))

Or say you want to annotate portions of the densities in color; the fill aesthetic can vary within a slab in all geoms and stats in the ggdist::geom_slabinterval() family, including ggdist::stat_halfeye(). For example, if you want to annotate a domain-specific region of practical equivalence (ROPE), you could do something like this:

m %>%
  spread_draws(condition_mean[condition]) %>%
  ggplot(aes(y = fct_rev(condition), x = condition_mean, fill = after_stat(abs(x) < .8))) +
  stat_halfeye() +
  geom_vline(xintercept = c(-.8, .8), linetype = "dashed") +
  scale_fill_manual(values = c("gray80", "skyblue"))

Other visualizations of distributions: stat_slabinterval

There are a variety of additional stats for visualizing distributions in the ggdist::geom_slabinterval() family of stats and geoms:

See vignette("slabinterval", package = "ggdist") for an overview.

Intervals with multiple probability levels: the .width = argument

If you wish to summarise the data before plotting (sometimes useful for large samples), median_qi() and its sister functions can also produce an arbitrary number of probability intervals by setting the .width = argument:

m %>%
  spread_draws(condition_mean[condition]) %>%
  median_qi(.width = c(.95, .8, .5))

The results are in a tidy format: one row per index (condition) and probability level (.width). This facilitates plotting. For example, assigning -.width to the linewidth aesthetic will show all intervals, making thicker lines correspond to smaller intervals:

m %>%
  spread_draws(condition_mean[condition]) %>%
  median_qi(.width = c(.95, .66)) %>%
  ggplot(aes(
    y = fct_rev(condition), x = condition_mean, xmin = .lower, xmax = .upper,
    # size = -.width means smaller probability interval => thicker line
    # this can be omitted, geom_pointinterval includes it automatically
    # if a .width column is in the input data.
    linewidth = -.width
  )) +  
  geom_pointinterval()

ggdist::geom_pointinterval() includes size = -.width as a default aesthetic mapping to facilitate exactly this usage.

Plotting posteriors as quantile dotplots

Intervals are nice if the alpha level happens to line up with whatever decision you are trying to make, but getting a shape of the posterior is better (hence eye plots, above). On the other hand, making inferences from density plots is imprecise (estimating the area of one shape as a proportion of another is a hard perceptual task). Reasoning about probability in frequency formats is easier, motivating quantile dotplots (Kay et al. 2016, Fernandes et al. 2018), which also allow precise estimation of arbitrary intervals (down to the dot resolution of the plot, 100 in the example below).

Within the slabinterval family of geoms in tidybayes is the dots and dotsinterval family, which automatically determine appropriate bin sizes for dotplots and can calculate quantiles from samples to construct quantile dotplots. ggdist::stat_dots() is the variant designed for use on samples:

m %>%
  spread_draws(condition_mean[condition]) %>%
  ggplot(aes(x = condition_mean, y = fct_rev(condition))) +
  stat_dotsinterval(quantiles = 100)

The idea is to get away from thinking about the posterior as indicating one canonical point or interval, but instead to represent it as (say) 100 approximately equally likely points.

Alternative point summaries and intervals: median, mean, mode; qi, hdi, hdci

The point_interval() family of functions follow the naming scheme [median|mean|mode]_[qi|hdi|hdci], and all work in the same way as median_qi(): they take a series of names (or expressions calculated on columns) and summarize those columns with the corresponding point summary function (median, mean, or mode) and interval (qi, hdi, or hdci). qi yields a quantile interval (a.k.a. equi-tailed interval, central interval, or percentile interval), hdi yields one or more highest (posterior) density interval(s), and hdci yields a single (possibly) highest-density continuous interval. These can be used in any combination desired.

The *_hdi functions have an additional difference: In the case of multimodal distributions, they may return multiple intervals for each probability level. Here are some draws from a multimodal normal mixture:

set.seed(123)
multimodal_draws = tibble(
  x = c(rnorm(5000, 0, 1), rnorm(2500, 4, 1))
)

Passed through mode_hdi(), we get multiple intervals at the 80% probability level:

multimodal_draws %>%
  mode_hdi(x, .width = .80)

This is easier to see when plotted:

multimodal_draws %>%
  ggplot(aes(x = x)) +
  stat_slab(aes(y = 0)) +
  stat_pointinterval(aes(y = -0.5), point_interval = median_qi, .width = c(.95, .80)) +
  annotate("text", label = "median, 80% and 95% quantile intervals", x = 6, y = -0.5, hjust = 0, vjust = 0.3) +
  stat_pointinterval(aes(y = -0.25), point_interval = mode_hdi, .width = c(.95, .80)) +
  annotate("text", label = "mode, 80% and 95% highest-density intervals", x = 6, y = -0.25, hjust = 0, vjust = 0.3) +
  xlim(-3.25, 18) +
  scale_y_continuous(breaks = NULL)

Posterior predictions with posterior distributions of means

Altogether, data, posterior predictions, and posterior distributions of the means:

draws = m %>%
  spread_draws(condition_mean[condition], response_sd)

reps = draws %>%
  mutate(y_rep = rnorm(n(), condition_mean, response_sd))

ABC %>%
  ggplot(aes(y = condition)) +
  stat_interval(aes(x = y_rep), .width = c(.95, .8, .5), data = reps) +
  stat_pointinterval(aes(x = condition_mean), .width = c(.95, .66), position = position_nudge(y = -0.3), data = draws) +
  geom_point(aes(x = response)) +
  scale_color_brewer()

Compatibility of point_interval with broom::tidy: A model comparison example

Combining to_broom_names() with median_qi() (or more generally, the point_interval() family of functions) makes it easy to compare results against models supported by broom::tidy(). For example, let’s compare our model’s fits for conditional means against an ordinary least squares (OLS) regression:

m_linear = lm(response ~ condition, data = ABC)

Combining emmeans::emmeans with broom::tidy, we can generate tidy-format summaries of conditional means from the above model:

linear_results = m_linear %>% 
  emmeans(~ condition) %>% 
  tidy(conf.int = TRUE) %>%
  mutate(model = "OLS")

linear_results

We can derive corresponding fits from our model:

bayes_results = m %>%
  spread_draws(condition_mean[condition]) %>%
  median_qi(estimate = condition_mean) %>%
  to_broom_names() %>%
  mutate(model = "Bayes")

bayes_results

Here, to_broom_names() will convert .lower and .upper into conf.low and conf.high so the names of the columns we need to make the comparison (condition, estimate, conf.low, and conf.high) all line up easily. This makes it simple to combine the two tidy data frames together using bind_rows, and plot them:

bind_rows(linear_results, bayes_results) %>%
  mutate(condition = fct_rev(condition)) %>%
  ggplot(aes(y = condition, x = estimate, xmin = conf.low, xmax = conf.high, color = model)) +
  geom_pointinterval(position = position_dodge(width = .3))

Observe the shrinkage towards the overall mean in the Bayesian model compared to the OLS model.

Compatibility with bayesplot using unspread_draws and ungather_draws

Functions from other packages might expect draws in the form of a data frame or matrix with variables as columns and draws as rows. That is the format returned by tidy_draws(), but not by gather_draws() or spread_draws(), which split indices from variables out into columns.

It may be desirable to use the spread_draws() or gather_draws() functions to transform your draws in some way, and then convert them back into the draw \(\times\) variable format to pass them into functions from other packages, like bayesplot. The unspread_draws() and ungather_draws() functions invert spread_draws() and gather_draws() to return a data frame with variable column names that include indices in them and draws as rows.

As an example, let’s re-do the previous example of compare_levels(), but use bayesplot::mcmc_areas() to plot the results instead of ggdist::stat_eye(). First, the result of compare_levels() looks like this:

m %>%
  spread_draws(condition_mean[condition]) %>%
  compare_levels(condition_mean, by = condition) %>%
  head(10)

To get a version we can pass to bayesplot::mcmc_areas(), all we need to do is invert the spread_draws() call we started with:

m %>%
  spread_draws(condition_mean[condition]) %>%
  compare_levels(condition_mean, by = condition) %>%
  unspread_draws(condition_mean[condition]) %>%
  head(10)

We can pass that into bayesplot::mcmc_areas() directly. The drop_indices = TRUE argument to unspread_draws() indicates that .chain, .iteration, and .draw should not be included in the output:

m %>%
  spread_draws(condition_mean[condition]) %>%
  compare_levels(condition_mean, by = condition) %>%
  unspread_draws(condition_mean[condition], drop_indices = TRUE) %>%
  bayesplot::mcmc_areas()

If you are instead working with tidy draws generated by gather_draws() or gather_variables(), the ungather_draws() function will transform those draws into the draw \(\times\) variable format. It has the same syntax as unspread_draws().

Compatibility with emmeans (formerly lsmeans)

The emmeans::emmeans() function provides a convenient syntax for generating marginal estimates from a model, including numerous types of contrasts. It also supports some Bayesian modeling packages, like MCMCglmm, rstanarm, and brms. However, it does not provide draws in a tidy format. The gather_emmeans_draws() function converts output from emmeans into a tidy format, keeping the emmeans reference grid and adding a .value column with long-format draws.

(Another approach to using emmeans contrast methods is to use emmeans_comparison() to convert emmeans contrast methods into comparison functions that can be used with tidybayes::compare_levels(). See the documentation of emmeans_comparison() for more information).

m_rst = stan_glm(response ~ condition, data = ABC)

m_rst %>%
  emmeans( ~ condition) %>%
  gather_emmeans_draws() %>%
  median_qi()

## Warning: Model has 0 prior weights, but we recovered 50 rows of data.
## So prior weights were ignored.

m_rst %>%
  emmeans( ~ condition) %>%
  contrast(method = "pairwise") %>%
  gather_emmeans_draws() %>%
  median_qi()

## Warning: Model has 0 prior weights, but we recovered 50 rows of data.
## So prior weights were ignored.

m_rst %>%
  emmeans( ~ condition) %>%
  contrast(method = "pairwise") %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = .value, y = contrast)) +
  stat_halfeye()

## Warning: Model has 0 prior weights, but we recovered 50 rows of data.
## So prior weights were ignored.

m_rst %>%
  emmeans(pairwise ~ condition) %>%
  gather_emmeans_draws() %>%
  median_qi()

## Warning: Model has 0 prior weights, but we recovered 50 rows of data.
## So prior weights were ignored.

# MCMCglmm does not support tibbles directly,
# so we convert ABC to a data.frame on the way in
m_glmm = MCMCglmm::MCMCglmm(response ~ condition, data = as.data.frame(ABC))

m_glmm %>%
  emmeans( ~ condition, data = ABC) %>%
  contrast(method = "pairwise") %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = .value, y = contrast)) +
  stat_halfeye()