- Introduction
- Setup
- Example dataset
- Using
`compose_data`

to prepare a data frame for the model - Extracting
draws from a fit in tidy-format using
`spread_draws`

- Point
summaries and intervals with the
`point_interval`

functions:`[median|mean|mode]_[qi|hdi]`

- Plotting points and intervals
- Using
`geom_pointinterval`

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

- Other
visualizations of distributions:
`stat_slabinterval`

- Intervals
with multiple probability levels: the
`.width =`

argument - Plotting posteriors as quantile dotplots
- Alternative point summaries and intervals: median, mean, mode; qi, hdi, hdci

- Using
- Combining variables with different indices in a single tidy format data frame
- Posterior predictions
- Comparing levels of a factor
- Gathering
all model variable names into a single column:
`gather_draws`

and`gather_variables`

- Selecting variables using regular expressions
- Drawing fit curves with uncertainty
- Compatibility with other packages

This vignette introduces the `tidybayes`

package, which
facilitates the use of tidy data (one observation per row) with Bayesian
models in R. This vignette is geared towards working with tidy data in
general-purpose modeling functions like JAGS or Stan. For a similar
introduction to the use of `tidybayes`

with high-level
modeling functions such as those in `brms`

or
`rstanarm`

, see `vignette("tidy-brms")`

or
`vignette("tidy-rstanarm")`

. This vignette also describes how
to use `ggdist`

(the sister package to
`tidybayes`

) for visualizing model output.

The default output (and sometimes input) data formats of popular
modeling functions like JAGS and Stan often don’t quite conform to the
ideal of tidy
data. For example, input formats might expect a list instead of a
data frame, and for all variables to be encoded as numeric values
(requiring translation of factors to numeric values and the creation of
index variables to store the number of levels per factor or the number
of observations in a data frame). Output formats will often be in matrix
form (requiring conversion for use with libraries like ggplot), and will
use numeric indices (requiring conversion back into factor level names
if the you wish to make meaningfully-labeled plots or tables).
`tidybayes`

automates all of these sorts of tasks.

There are a few core ideas that run through the
`tidybayes`

API that should (hopefully) make it easy to
use:

**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.**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.**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`geom`

s) 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).**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.**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.

`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"`

.

The following libraries are required to run this vignette:

```
library(magrittr)
library(dplyr)
library(forcats)
library(modelr)
library(ggdist)
library(tidybayes)
library(ggplot2)
library(cowplot)
library(emmeans)
library(broom)
library(rstan)
library(rstanarm)
library(brms)
library(bayesplot)
library(RColorBrewer)
theme_set(theme_tidybayes() + panel_border())
```

These options help Stan run faster:

To demonstrate `tidybayes`

, we will use a simple dataset
with 10 observations from 5 conditions each:

```
set.seed(5)
n = 10
n_condition = 5
ABC =
tibble(
condition = factor(rep(c("A","B","C","D","E"), n)),
response = rnorm(n * 5, c(0,1,2,1,-1), 0.5)
)
```

A snapshot of the data looks like this:

condition | response |
---|---|

A | -0.4204277 |

B | 1.6921797 |

C | 1.3722541 |

D | 1.0350714 |

E | -0.1442796 |

A | -0.3014540 |

B | 0.7639168 |

C | 1.6823143 |

D | 0.8571132 |

E | -0.9309459 |

This is a typical tidy format data frame: one observation per row. Graphically:

`compose_data`

to prepare a data frame for the
modelShunting data from a data frame into a format usable in samplers like
JAGS or Stan can involve a tedious set of operations, like generating
index variables storing the number of operations or the number of levels
in a factor. `compose_data`

automates these operations.

A hierarchical model of our example data might fit an overall mean
across the conditions (`overall_mean`

), the standard
deviation of the condition means (`condition_mean_sd`

), the
mean within each condition (`condition_mean[condition]`

) and
the standard deviation of the responses given a condition mean
(`response_sd`

):

```
data {
int<lower=1> n;
int<lower=1> n_condition;
int<lower=1, upper=n_condition> condition[n];
real response[n];
}
parameters {
real overall_mean;
vector[n_condition] condition_zoffset;
real<lower=0> response_sd;
real<lower=0> condition_mean_sd;
}
transformed parameters {
vector[n_condition] condition_mean;
condition_mean = overall_mean + condition_zoffset * condition_mean_sd;
}
model {
response_sd ~ cauchy(0, 1); // => half-cauchy(0, 1)
condition_mean_sd ~ cauchy(0, 1); // => half-cauchy(0, 1)
overall_mean ~ normal(0, 5);
condition_zoffset ~ normal(0, 1); // => condition_mean ~ normal(overall_mean, condition_mean_sd)
for (i in 1:n) {
response[i] ~ normal(condition_mean[condition[i]], response_sd);
}
}
```

We have compiled and loaded this model into the variable
`ABC_stan`

.

This model expects these variables as input:

`n`

: number of observations`n_condition`

: number of conditions`condition`

: a vector of integers indicating the condition of each observation`response`

: a vector of observations

Our data frame (`ABC`

) only has `response`

and
`condition`

, and `condition`

is in the wrong
format (it is a factor instead of numeric). However,
`compose_data`

can generate a list containing the above
variables in the correct format automatically. It recognizes that
`condition`

is a factor and converts it to a numeric, adds
the `n_condition`

variable automatically containing the
number of levels in `condition`

, and adds the `n`

column containing the number of observations (number of rows in the data
frame):

```
## $condition
## [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
##
## $n_condition
## [1] 5
##
## $response
## [1] -0.42042774 1.69217967 1.37225407 1.03507138 -0.14427956 -0.30145399 0.76391681 1.68231434 0.85711318
## [10] -0.93094589 0.61381517 0.59911027 1.45980370 0.92123282 -1.53588002 -0.06949307 0.70134345 0.90801662
## [19] 1.12040863 -1.12967770 0.45025597 1.47093470 2.73398095 1.35338054 -0.59049553 -0.14674092 1.70929454
## [28] 2.74938691 0.67145895 -1.42639772 0.15795752 1.55484708 3.10773029 1.60855182 -0.26038911 0.47578692
## [37] 0.49523368 0.99976363 0.11890706 -1.07130406 0.77503018 0.59878841 1.96271054 1.94783398 -1.22828447
## [46] 0.28111168 0.55649574 1.76987771 0.63783576 -1.03460558
##
## $n
## [1] 50
```

This makes it easy to skip right to running the model without munging the data yourself:

The results look like this:

```
## Inference for Stan model: 150269bbfe815db848542eb814742386.
## 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
## overall_mean 0.60 0.03 0.72 -0.75 0.27 0.61 0.96 1.83 491 1.01
## condition_zoffset[1] -0.39 0.02 0.51 -1.40 -0.72 -0.38 -0.05 0.57 795 1.00
## condition_zoffset[2] 0.35 0.02 0.49 -0.59 0.02 0.36 0.69 1.33 889 1.00
## condition_zoffset[3] 1.13 0.02 0.60 0.01 0.72 1.13 1.53 2.33 865 1.00
## condition_zoffset[4] 0.37 0.02 0.49 -0.62 0.05 0.38 0.71 1.32 894 1.00
## condition_zoffset[5] -1.39 0.02 0.69 -2.80 -1.82 -1.36 -0.93 -0.11 782 1.01
## response_sd 0.56 0.00 0.06 0.46 0.52 0.56 0.60 0.70 1603 1.00
## condition_mean_sd 1.25 0.02 0.60 0.61 0.88 1.10 1.43 2.81 592 1.01
## condition_mean[1] 0.20 0.00 0.17 -0.14 0.08 0.20 0.31 0.54 4265 1.00
## condition_mean[2] 1.00 0.00 0.18 0.65 0.89 1.00 1.12 1.34 4434 1.00
## condition_mean[3] 1.84 0.00 0.18 1.48 1.72 1.84 1.96 2.19 4501 1.00
## condition_mean[4] 1.02 0.00 0.18 0.66 0.90 1.02 1.14 1.35 4631 1.00
## condition_mean[5] -0.88 0.00 0.18 -1.23 -1.00 -0.89 -0.77 -0.53 4064 1.00
## lp__ 0.22 0.08 2.49 -5.80 -1.18 0.59 2.03 3.94 873 1.00
##
## Samples were drawn using NUTS(diag_e) at Sat Aug 12 16:10:34 2023.
## 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).
```

`spread_draws`

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:

```
## 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:

condition | condition_mean | .chain | .iteration | .draw |
---|---|---|---|---|

1 | 0.2342760 | 1 | 1 | 1 |

1 | 0.1558843 | 1 | 2 | 2 |

1 | 0.0598095 | 1 | 3 | 3 |

1 | 0.4264011 | 1 | 4 | 4 |

1 | 0.0163756 | 1 | 5 | 5 |

1 | 0.1309990 | 1 | 6 | 6 |

1 | 0.3973789 | 1 | 7 | 7 |

1 | 0.3245946 | 1 | 8 | 8 |

1 | 0.3626703 | 1 | 9 | 9 |

1 | 0.1524349 | 1 | 10 | 10 |

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:

condition | condition_mean | .chain | .iteration | .draw |
---|---|---|---|---|

A | 0.2342760 | 1 | 1 | 1 |

A | 0.1558843 | 1 | 2 | 2 |

A | 0.0598095 | 1 | 3 | 3 |

A | 0.4264011 | 1 | 4 | 4 |

A | 0.0163756 | 1 | 5 | 5 |

A | 0.1309990 | 1 | 6 | 6 |

A | 0.3973789 | 1 | 7 | 7 |

A | 0.3245946 | 1 | 8 | 8 |

A | 0.3626703 | 1 | 9 | 9 |

A | 0.1524349 | 1 | 10 | 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:

Now we can omit the `recover_types`

call before subsequent
calls to `spread_draws`

.

`point_interval`

functions: `[median|mean|mode]_[qi|hdi]`

`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:

.chain | .iteration | .draw | overall_mean | response_sd |
---|---|---|---|---|

1 | 1 | 1 | -0.2899471 | 0.6395042 |

1 | 2 | 2 | -0.2075163 | 0.4545561 |

1 | 3 | 3 | 0.3475639 | 0.7130237 |

1 | 4 | 4 | 0.5457579 | 0.6811170 |

1 | 5 | 5 | 0.3380617 | 0.6259130 |

1 | 6 | 6 | 0.3352196 | 0.5012065 |

1 | 7 | 7 | 0.9522320 | 0.5374218 |

1 | 8 | 8 | 1.0902141 | 0.5430605 |

1 | 9 | 9 | 1.1868684 | 0.5044625 |

1 | 10 | 10 | 1.8497308 | 0.6200319 |

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`

:

overall_mean | overall_mean.lower | overall_mean.upper | response_sd | response_sd.lower | response_sd.upper | .width | .point | .interval |
---|---|---|---|---|---|---|---|---|

0.6122377 | -0.7461883 | 1.828706 | 0.556411 | 0.4568772 | 0.7024135 | 0.95 | median | qi |

`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:

overall_mean | overall_mean.lower | overall_mean.upper | response_sd | response_sd.lower | response_sd.upper | .width | .point | .interval |
---|---|---|---|---|---|---|---|---|

0.6122377 | -0.7461883 | 1.828706 | 0.556411 | 0.4568772 | 0.7024135 | 0.95 | median | qi |

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:

condition | condition_mean | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|

A | 0.1953174 | -0.1413499 | 0.5398734 | 0.95 | median | qi |

B | 0.9992510 | 0.6468258 | 1.3395784 | 0.95 | median | qi |

C | 1.8436750 | 1.4817705 | 2.1910130 | 0.95 | median | qi |

D | 1.0184200 | 0.6625502 | 1.3541714 | 0.95 | median | qi |

E | -0.8860573 | -1.2293855 | -0.5281050 | 0.95 | 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)
```

condition | condition_mean | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|

A | 0.1953174 | -0.1413499 | 0.5398734 | 0.95 | median | qi |

B | 0.9992510 | 0.6468258 | 1.3395784 | 0.95 | median | qi |

C | 1.8436750 | 1.4817705 | 2.1910130 | 0.95 | median | qi |

D | 1.0184200 | 0.6625502 | 1.3541714 | 0.95 | median | qi |

E | -0.8860573 | -1.2293855 | -0.5281050 | 0.95 | median | qi |

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:

condition | variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
---|---|---|---|---|---|---|---|---|---|---|

A | condition_mean | 0.1964118 | 0.1953174 | 0.1741579 | 0.1733574 | -0.0825453 | 0.4849402 | 1.0005379 | 4405.430 | 3092.188 |

B | condition_mean | 1.0012177 | 0.9992510 | 0.1764290 | 0.1719944 | 0.7150631 | 1.2905922 | 1.0004779 | 4462.506 | 2991.338 |

C | condition_mean | 1.8416926 | 1.8436750 | 0.1801959 | 0.1772218 | 1.5441089 | 2.1331651 | 1.0001609 | 4584.273 | 3310.993 |

D | condition_mean | 1.0182492 | 1.0184200 | 0.1756338 | 0.1735525 | 0.7267896 | 1.3014954 | 1.0004832 | 4653.303 | 3307.174 |

E | condition_mean | -0.8846496 | -0.8860573 | 0.1778983 | 0.1745175 | -1.1772792 | -0.5881216 | 0.9998871 | 4188.030 | 2925.845 |

`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()`

.

`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:

`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.

`.width =`

argumentIf 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:

condition | condition_mean | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|

A | 0.1953174 | -0.1413499 | 0.5398734 | 0.95 | median | qi |

B | 0.9992510 | 0.6468258 | 1.3395784 | 0.95 | median | qi |

C | 1.8436750 | 1.4817705 | 2.1910130 | 0.95 | median | qi |

D | 1.0184200 | 0.6625502 | 1.3541714 | 0.95 | median | qi |

E | -0.8860573 | -1.2293855 | -0.5281050 | 0.95 | median | qi |

A | 0.1953174 | -0.0249085 | 0.4140537 | 0.80 | median | qi |

B | 0.9992510 | 0.7786621 | 1.2245413 | 0.80 | median | qi |

C | 1.8436750 | 1.6120169 | 2.0742256 | 0.80 | median | qi |

D | 1.0184200 | 0.7978164 | 1.2478770 | 0.80 | median | qi |

E | -0.8860573 | -1.1113424 | -0.6609991 | 0.80 | median | qi |

A | 0.1953174 | 0.0793111 | 0.3130674 | 0.50 | median | qi |

B | 0.9992510 | 0.8875942 | 1.1199125 | 0.50 | median | qi |

C | 1.8436750 | 1.7217698 | 1.9605066 | 0.50 | median | qi |

D | 1.0184200 | 0.9013932 | 1.1356234 | 0.50 | median | qi |

E | -0.8860573 | -1.0017830 | -0.7658796 | 0.50 | median | qi |

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.

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.

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:

Passed through `mode_hdi()`

, we get multiple intervals at
the 80% probability level:

x | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|

-0.0605292 | -1.455675 | 1.540140 | 0.8 | mode | hdi |

-0.0605292 | 3.106254 | 5.005495 | 0.8 | mode | hdi |

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)
```

`spread_draws()`

supports extracting variables that have
different indices. It automatically matches up indices with the same
name, and duplicates values as necessary to produce one row per all
combination of levels of all indices. For example, we might want to
calculate the difference between each condition mean and the overall
mean. To do that, we can extract draws from the overall mean and all
condition means:

.chain | .iteration | .draw | overall_mean | condition | condition_mean |
---|---|---|---|---|---|

1 | 1 | 1 | -0.2899471 | A | 0.2342760 |

1 | 1 | 1 | -0.2899471 | B | 1.0707836 |

1 | 1 | 1 | -0.2899471 | C | 1.8049518 |

1 | 1 | 1 | -0.2899471 | D | 0.9431367 |

1 | 1 | 1 | -0.2899471 | E | -0.9705730 |

1 | 2 | 2 | -0.2075163 | A | 0.1558843 |

1 | 2 | 2 | -0.2075163 | B | 1.0581285 |

1 | 2 | 2 | -0.2075163 | C | 1.8634087 |

1 | 2 | 2 | -0.2075163 | D | 0.6074677 |

1 | 2 | 2 | -0.2075163 | E | -1.0961247 |

Within each draw, `overall_mean`

is repeated as necessary
to correspond to every index of `condition_mean`

. Thus, the
`dplyr::mutate()`

function can be used to take the
differences over all rows, then we can summarize with
`median_qi()`

:

```
m %>%
spread_draws(overall_mean, condition_mean[condition]) %>%
mutate(condition_offset = condition_mean - overall_mean) %>%
median_qi(condition_offset)
```

condition | condition_offset | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|

A | -0.4089436 | -1.6494621 | 0.9427468 | 0.95 | median | qi |

B | 0.3907151 | -0.8098025 | 1.6987681 | 0.95 | median | qi |

C | 1.2274271 | 0.0143030 | 2.6095415 | 0.95 | median | qi |

D | 0.4101439 | -0.8225706 | 1.7427413 | 0.95 | median | qi |

E | -1.4865013 | -2.6990185 | -0.1808296 | 0.95 | median | qi |

We can use combinations of variables with difference indices to
generate predictions from the model. In this case, we can combine the
condition means with the residual standard deviation to generate
predictive distributions from the model, then show the distributions
using `ggdist::stat_interval()`

and compare them to the
data:

```
m %>%
spread_draws(condition_mean[condition], response_sd) %>%
mutate(y_rep = rnorm(n(), condition_mean, response_sd)) %>%
median_qi(y_rep, .width = c(.95, .8, .5)) %>%
ggplot(aes(y = fct_rev(condition), x = y_rep)) +
geom_interval(aes(xmin = .lower, xmax = .upper)) + #auto-sets aes(color = fct_rev(ordered(.width)))
geom_point(aes(x = response), data = ABC) +
scale_color_brewer()
```

If this model is well-calibrated, about 95% of the data should be within the outer intervals, 80% in the next-smallest intervals, and 50% in the smallest intervals.

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()
```

`compare_levels()`

allows us to compare the value of some
variable across levels of some factor. By default it computes all
pairwise differences, though this can be changed using the
`comparison =`

argument:

`gather_draws`

and `gather_variables`

We might also prefer all model variable names to be in a single
column (long-format) instead of as column names. There are two methods
for obtaining long-format data frames with `tidybayes`

, whose
use depends on where and how in the data processing chain you might want
to transform into long-format: `gather_draws()`

and
`gather_variables()`

. There are also two methods for wide (or
semi-wide) format data frame, `spread_draws()`

(described
previously) and `tidy_draws()`

.

`gather_draws()`

is the counterpart to
`spread_draws()`

, except it puts all variable names in a
`.variable`

column and all draws in a `.value`

column:

.variable | condition | .value | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|---|

condition_mean | A | 0.1953174 | -0.1413499 | 0.5398734 | 0.95 | median | qi |

condition_mean | B | 0.9992510 | 0.6468258 | 1.3395784 | 0.95 | median | qi |

condition_mean | C | 1.8436750 | 1.4817705 | 2.1910130 | 0.95 | median | qi |

condition_mean | D | 1.0184200 | 0.6625502 | 1.3541714 | 0.95 | median | qi |

condition_mean | E | -0.8860573 | -1.2293855 | -0.5281050 | 0.95 | median | qi |

overall_mean | NA | 0.6122377 | -0.7461883 | 1.8287064 | 0.95 | median | qi |

Note that `condition = NA`

for the
`overall_mean`

row, because it does not have an index with
that name in the specification passed to
`gather_draws()`

.

While this works well if we do not need to perform computations that
involve multiple columns, the semi-wide format returned by
`spread_draws()`

is very useful for computations that involve
multiple columns names, such as the calculation of the
`condition_offset`

above. If we want to make intermediate
computations on the format returned by `spread_draws`

and
*then* gather variables into one column, we can use
`gather_variables()`

, which will gather all non-grouped
variables that are not special columns (like `.chain`

,
`.iteration`

, or `.draw`

):

```
m %>%
spread_draws(overall_mean, condition_mean[condition]) %>%
mutate(condition_offset = condition_mean - overall_mean) %>%
gather_variables() %>%
median_qi()
```

condition | .variable | .value | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|---|

A | condition_mean | 0.1953174 | -0.1413499 | 0.5398734 | 0.95 | median | qi |

A | condition_offset | -0.4089436 | -1.6494621 | 0.9427468 | 0.95 | median | qi |

A | overall_mean | 0.6122377 | -0.7461883 | 1.8287064 | 0.95 | median | qi |

B | condition_mean | 0.9992510 | 0.6468258 | 1.3395784 | 0.95 | median | qi |

B | condition_offset | 0.3907151 | -0.8098025 | 1.6987681 | 0.95 | median | qi |

B | overall_mean | 0.6122377 | -0.7461883 | 1.8287064 | 0.95 | median | qi |

C | condition_mean | 1.8436750 | 1.4817705 | 2.1910130 | 0.95 | median | qi |

C | condition_offset | 1.2274271 | 0.0143030 | 2.6095415 | 0.95 | median | qi |

C | overall_mean | 0.6122377 | -0.7461883 | 1.8287064 | 0.95 | median | qi |

D | condition_mean | 1.0184200 | 0.6625502 | 1.3541714 | 0.95 | median | qi |

D | condition_offset | 0.4101439 | -0.8225706 | 1.7427413 | 0.95 | median | qi |

D | overall_mean | 0.6122377 | -0.7461883 | 1.8287064 | 0.95 | median | qi |

E | condition_mean | -0.8860573 | -1.2293855 | -0.5281050 | 0.95 | median | qi |

E | condition_offset | -1.4865013 | -2.6990185 | -0.1808296 | 0.95 | median | qi |

E | overall_mean | 0.6122377 | -0.7461883 | 1.8287064 | 0.95 | median | qi |

Note how `overall_mean`

is now repeated here for each
condition, because we have performed the gather after spreading model
variables across columns.

Finally, if we want raw model variable names as columns names instead
of having indices split out as their own column names, we can use
`tidy_draws()`

. Generally speaking
`spread_draws()`

and `gather_draws()`

are
typically more useful than `tidy_draws()`

, but it is provided
as a common method for generating data frames from many types of
Bayesian models, and is used internally by `gather_draws()`

and `spread_draws()`

:

.chain | .iteration | .draw | overall_mean | condition_zoffset[1] | condition_zoffset[2] | condition_zoffset[3] | condition_zoffset[4] | condition_zoffset[5] | response_sd | condition_mean_sd | condition_mean[1] | condition_mean[2] | condition_mean[3] | condition_mean[4] | condition_mean[5] | lp__ | accept_stat__ | stepsize__ | treedepth__ | n_leapfrog__ | divergent__ | energy__ |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

1 | 1 | 1 | -0.2899471 | 0.4831306 | 1.2540664 | 1.9306850 | 1.1364255 | -0.6272734 | 0.6395042 | 1.085055 | 0.2342760 | 1.0707836 | 1.804952 | 0.9431367 | -0.9705730 | 0.5108208 | 0.9979826 | 0.044034 | 5 | 31 | 0 | 9.683216 |

1 | 2 | 2 | -0.2075163 | 0.1399514 | 0.4874201 | 0.7975465 | 0.3138634 | -0.3422174 | 0.4545561 | 2.596620 | 0.1558843 | 1.0581285 | 1.863409 | 0.6074677 | -1.0961247 | -0.8374519 | 0.9889947 | 0.044034 | 7 | 127 | 0 | 2.236431 |

1 | 3 | 3 | 0.3475639 | -0.2265616 | 0.2707671 | 1.2311989 | 0.5052696 | -0.8482141 | 0.7130237 | 1.270093 | 0.0598095 | 0.6914634 | 1.911301 | 0.9893035 | -0.7297473 | -0.7632896 | 0.9987001 | 0.044034 | 6 | 79 | 0 | 4.686032 |

1 | 4 | 4 | 0.5457579 | -0.0857109 | 0.3658783 | 0.9256314 | 0.3378245 | -1.2429921 | 0.6811170 | 1.392551 | 0.4264011 | 1.0552620 | 1.834747 | 1.0161956 | -1.1851718 | 0.4434354 | 0.9988460 | 0.044034 | 5 | 31 | 0 | 2.496561 |

1 | 5 | 5 | 0.3380617 | -0.2204426 | 0.3971486 | 1.1010441 | 0.4589120 | -0.8072970 | 0.6259130 | 1.459274 | 0.0163756 | 0.9176103 | 1.944786 | 1.0077399 | -0.8400057 | 2.8357268 | 0.9972974 | 0.044034 | 6 | 127 | 0 | 1.452722 |

1 | 6 | 6 | 0.3352196 | -0.1467775 | 0.4757227 | 1.2092398 | 0.5614834 | -0.9003408 | 0.5012065 | 1.391362 | 0.1309990 | 0.9971222 | 2.017710 | 1.1164465 | -0.9174806 | 3.9459284 | 0.9952214 | 0.044034 | 6 | 63 | 0 | -1.484267 |

1 | 7 | 7 | 0.9522320 | -0.3245970 | 0.0292202 | 0.5217779 | -0.0166794 | -1.2163993 | 0.5374218 | 1.709360 | 0.3973789 | 1.0021798 | 1.844138 | 0.9237210 | -1.1270323 | 3.3700958 | 0.9920297 | 0.044034 | 7 | 159 | 0 | -1.450467 |

1 | 8 | 8 | 1.0902141 | -0.4231561 | -0.0262401 | 0.7844110 | 0.0389688 | -1.1048486 | 0.5430605 | 1.809308 | 0.3245946 | 1.0427377 | 2.509455 | 1.1607207 | -0.9087971 | -2.6385535 | 0.8829194 | 0.044034 | 3 | 15 | 0 | 5.682317 |

1 | 9 | 9 | 1.1868684 | -0.4837738 | 0.0986432 | 0.4706185 | 0.0653861 | -1.2278957 | 0.5044625 | 1.703685 | 0.3626703 | 1.3549252 | 1.988654 | 1.2982656 | -0.9050788 | 0.2801249 | 0.9991187 | 0.044034 | 4 | 15 | 0 | 5.287994 |

1 | 10 | 10 | 1.8497308 | -1.3883017 | -0.8230590 | -0.2118896 | -0.7794664 | -2.2329638 | 0.6200319 | 1.222570 | 0.1524349 | 0.8434835 | 1.590681 | 0.8967786 | -0.8802236 | -1.0706447 | 1.0000000 | 0.044034 | 7 | 127 | 0 | 3.482374 |

Combining `tidy_draws()`

with
`gather_variables()`

also allows us to derive similar output
to `ggmcmc::ggs()`

, if desired:

.chain | .iteration | .draw | .variable | .value |
---|---|---|---|---|

1 | 1 | 1 | overall_mean | -0.2899471 |

1 | 2 | 2 | overall_mean | -0.2075163 |

1 | 3 | 3 | overall_mean | 0.3475639 |

1 | 4 | 4 | overall_mean | 0.5457579 |

1 | 5 | 5 | overall_mean | 0.3380617 |

1 | 6 | 6 | overall_mean | 0.3352196 |

1 | 7 | 7 | overall_mean | 0.9522320 |

1 | 8 | 8 | overall_mean | 1.0902141 |

1 | 9 | 9 | overall_mean | 1.1868684 |

1 | 10 | 10 | overall_mean | 1.8497308 |

But again, this approach does not handle variable indices for us
automatically, so using `spread_draws()`

and
`gather_draws()`

is generally recommended unless you do not
have variable indices to worry about.

You can use regular expressions in the specifications passed to
`spread_draws()`

and `gather_draws()`

to match
multiple columns by passing `regex = TRUE`

. Our example fit
contains variables named `condition_mean[i]`

and
`condition_zoffset[i]`

. We could extract both using a single
regular expression:

condition | condition_mean | condition_zoffset | .chain | .iteration | .draw |
---|---|---|---|---|---|

A | 0.2342760 | 0.4831306 | 1 | 1 | 1 |

A | 0.1558843 | 0.1399514 | 1 | 2 | 2 |

A | 0.0598095 | -0.2265616 | 1 | 3 | 3 |

A | 0.4264011 | -0.0857109 | 1 | 4 | 4 |

A | 0.0163756 | -0.2204426 | 1 | 5 | 5 |

A | 0.1309990 | -0.1467775 | 1 | 6 | 6 |

A | 0.3973789 | -0.3245970 | 1 | 7 | 7 |

A | 0.3245946 | -0.4231561 | 1 | 8 | 8 |

A | 0.3626703 | -0.4837738 | 1 | 9 | 9 |

A | 0.1524349 | -1.3883017 | 1 | 10 | 10 |

This result is equivalent in this case to
`spread_draws(c(condition_mean, condition_zoffset)[condition])`

,
but does not require us to list each variable explicitly—this can be
useful, for example, in models with naming schemes like
`b_[some name]`

for coefficients.

To demonstrate drawing fit curves with uncertainty, let’s fit a
slightly naive model to part of the `mtcars`

dataset using
`brms::brm()`

:

```
m_mpg = brm(
mpg ~ hp * cyl,
data = mtcars,
file = "models/tidybayes_m_mpg.rds" # cache model (can be removed)
)
```

We can draw fit curves with probability bands using
`add_epred_draws()`

and
`ggdist::stat_lineribbon()`

:

```
mtcars %>%
group_by(cyl) %>%
data_grid(hp = seq_range(hp, n = 51)) %>%
add_epred_draws(m_mpg) %>%
ggplot(aes(x = hp, y = mpg, color = ordered(cyl))) +
stat_lineribbon(aes(y = .epred)) +
geom_point(data = mtcars) +
scale_fill_brewer(palette = "Greys") +
scale_color_brewer(palette = "Set2")
```

Or we can sample a reasonable number of fit lines (say 100) and overplot them:

```
mtcars %>%
group_by(cyl) %>%
data_grid(hp = seq_range(hp, n = 101)) %>%
# NOTE: this shows the use of ndraws to subsample within add_epred_draws()
# ONLY do this IF you are planning to make spaghetti plots, etc.
# NEVER subsample to a small sample to plot intervals, densities, etc.
add_epred_draws(m_mpg, ndraws = 100) %>%
ggplot(aes(x = hp, y = mpg, color = ordered(cyl))) +
geom_line(aes(y = .epred, group = paste(cyl, .draw)), alpha = .1) +
geom_point(data = mtcars) +
scale_color_brewer(palette = "Dark2")
```

For more examples of fit line uncertainty, see the corresponding
sections in `vignette("tidy-brms")`

or
`vignette("tidy-rstanarm")`

.

`point_interval`

with
`broom::tidy`

: A model comparison exampleCombining `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:

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
```

condition | estimate | std.error | df | conf.low | conf.high | statistic | p.value | model |
---|---|---|---|---|---|---|---|---|

A | 0.1815842 | 0.173236 | 45 | -0.1673310 | 0.5304993 | 1.048190 | 0.3001485 | OLS |

B | 1.0142144 | 0.173236 | 45 | 0.6652993 | 1.3631296 | 5.854526 | 0.0000005 | OLS |

C | 1.8745839 | 0.173236 | 45 | 1.5256687 | 2.2234990 | 10.820985 | 0.0000000 | OLS |

D | 1.0271794 | 0.173236 | 45 | 0.6782642 | 1.3760946 | 5.929366 | 0.0000004 | OLS |

E | -0.9352260 | 0.173236 | 45 | -1.2841411 | -0.5863108 | -5.398567 | 0.0000024 | OLS |

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
```

condition | estimate | conf.low | conf.high | .width | .point | .interval | model |
---|---|---|---|---|---|---|---|

A | 0.1953174 | -0.1413499 | 0.5398734 | 0.95 | median | qi | Bayes |

B | 0.9992510 | 0.6468258 | 1.3395784 | 0.95 | median | qi | Bayes |

C | 1.8436750 | 1.4817705 | 2.1910130 | 0.95 | median | qi | Bayes |

D | 1.0184200 | 0.6625502 | 1.3541714 | 0.95 | median | qi | Bayes |

E | -0.8860573 | -1.2293855 | -0.5281050 | 0.95 | median | qi | Bayes |

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.

`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)
```

.chain | .iteration | .draw | condition | condition_mean |
---|---|---|---|---|

1 | 1 | 1 | B - A | 0.8365076 |

1 | 2 | 2 | B - A | 0.9022442 |

1 | 3 | 3 | B - A | 0.6316538 |

1 | 4 | 4 | B - A | 0.6288609 |

1 | 5 | 5 | B - A | 0.9012346 |

1 | 6 | 6 | B - A | 0.8661232 |

1 | 7 | 7 | B - A | 0.6048009 |

1 | 8 | 8 | B - A | 0.7181431 |

1 | 9 | 9 | B - A | 0.9922549 |

1 | 10 | 10 | B - A | 0.6910487 |

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)
```

.chain | .iteration | .draw | condition_mean[B - A] | condition_mean[C - A] | condition_mean[C - B] | condition_mean[D - A] | condition_mean[D - B] | condition_mean[D - C] | condition_mean[E - A] | condition_mean[E - B] | condition_mean[E - C] | condition_mean[E - D] |
---|---|---|---|---|---|---|---|---|---|---|---|---|

1 | 1 | 1 | 0.8365076 | 1.570676 | 0.7341682 | 0.7088607 | -0.1276469 | -0.8618150 | -1.2048490 | -2.041357 | -2.775525 | -1.913710 |

1 | 2 | 2 | 0.9022442 | 1.707524 | 0.8052802 | 0.4515835 | -0.4506608 | -1.2559409 | -1.2520090 | -2.154253 | -2.959533 | -1.703593 |

1 | 3 | 3 | 0.6316538 | 1.851492 | 1.2198381 | 0.9294940 | 0.2978401 | -0.9219980 | -0.7895568 | -1.421211 | -2.641049 | -1.719051 |

1 | 4 | 4 | 0.6288609 | 1.408346 | 0.7794846 | 0.5897945 | -0.0390664 | -0.8185510 | -1.6115729 | -2.240434 | -3.019919 | -2.201367 |

1 | 5 | 5 | 0.9012346 | 1.928411 | 1.0271762 | 0.9913643 | 0.0901297 | -0.9370466 | -0.8563813 | -1.757616 | -2.784792 | -1.847746 |

1 | 6 | 6 | 0.8661232 | 1.886711 | 1.0205881 | 0.9854475 | 0.1193243 | -0.9012638 | -1.0484796 | -1.914603 | -2.935191 | -2.033927 |

1 | 7 | 7 | 0.6048009 | 1.446759 | 0.8419585 | 0.5263421 | -0.0784588 | -0.9204173 | -1.5244112 | -2.129212 | -2.971171 | -2.050753 |

1 | 8 | 8 | 0.7181431 | 2.184860 | 1.4667173 | 0.8361261 | 0.1179830 | -1.3487343 | -1.2333917 | -1.951535 | -3.418252 | -2.069518 |

1 | 9 | 9 | 0.9922549 | 1.625984 | 0.6337288 | 0.9355953 | -0.0566596 | -0.6903884 | -1.2677492 | -2.260004 | -2.893733 | -2.203345 |

1 | 10 | 10 | 0.6910487 | 1.438246 | 0.7471974 | 0.7443437 | 0.0532951 | -0.6939024 | -1.0326585 | -1.723707 | -2.470905 | -1.777002 |

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()`

.

`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).

`rstanarm`

or `brms`

Both `rstanarm`

and `brms`

behave similarly
when used with `emmeans::emmeans()`

. The example below uses
`rstanarm`

, but should work similarly for
`brms`

.

Given this `rstanarm`

model:

We can use `emmeans::emmeans()`

to get conditional means
with uncertainty:

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

condition | .value | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|

A | 0.1872551 | -0.1916132 | 0.5332921 | 0.95 | median | qi |

B | 1.0163400 | 0.6705062 | 1.3722505 | 0.95 | median | qi |

C | 1.8750251 | 1.5210300 | 2.2300568 | 0.95 | median | qi |

D | 1.0262675 | 0.6682517 | 1.3754806 | 0.95 | median | qi |

E | -0.9405232 | -1.2748228 | -0.5901560 | 0.95 | median | qi |

Or `emmeans::emmeans()`

with
`emmeans::contrast()`

to do all pairwise comparisons:

```
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.
```

contrast | .value | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|

A - B | -0.8329464 | -1.3416973 | -0.3433737 | 0.95 | median | qi |

A - C | -1.6885546 | -2.1906835 | -1.1999622 | 0.95 | median | qi |

A - D | -0.8408330 | -1.3353113 | -0.3206278 | 0.95 | median | qi |

A - E | 1.1225019 | 0.5992254 | 1.6161816 | 0.95 | median | qi |

B - C | -0.8609638 | -1.3673337 | -0.3440654 | 0.95 | median | qi |

B - D | -0.0088381 | -0.4878745 | 0.4841102 | 0.95 | median | qi |

B - E | 1.9508883 | 1.4587024 | 2.4416107 | 0.95 | median | qi |

C - D | 0.8464505 | 0.3731628 | 1.3487615 | 0.95 | median | qi |

C - E | 2.8086519 | 2.3257780 | 3.3000857 | 0.95 | median | qi |

D - E | 1.9641383 | 1.4606453 | 2.4563015 | 0.95 | median | qi |

See the documentation for `emmeans::pairwise.emmc()`

for a
list of the numerous contrast types supported by
`emmeans::emmeans()`

.

As before, we can plot the results instead of using a table:

```
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.
```

`gather_emmeans_draws()`

also supports
`emm_list`

objects, which contain estimates from different
reference grids (see `emmeans::ref_grid()`

for more
information on reference grids). An additional column with the default
name of `.grid`

is added to indicate the reference grid for
each row in the output:

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

condition | contrast | .grid | .value | .lower | .upper | .width | .point | .interval |
---|---|---|---|---|---|---|---|---|

A | NA | emmeans | 0.1872551 | -0.1916132 | 0.5332921 | 0.95 | median | qi |

B | NA | emmeans | 1.0163400 | 0.6705062 | 1.3722505 | 0.95 | median | qi |

C | NA | emmeans | 1.8750251 | 1.5210300 | 2.2300568 | 0.95 | median | qi |

D | NA | emmeans | 1.0262675 | 0.6682517 | 1.3754806 | 0.95 | median | qi |

E | NA | emmeans | -0.9405232 | -1.2748228 | -0.5901560 | 0.95 | median | qi |

NA | A - B | contrasts | -0.8329464 | -1.3416973 | -0.3433737 | 0.95 | median | qi |

NA | A - C | contrasts | -1.6885546 | -2.1906835 | -1.1999622 | 0.95 | median | qi |

NA | A - D | contrasts | -0.8408330 | -1.3353113 | -0.3206278 | 0.95 | median | qi |

NA | A - E | contrasts | 1.1225019 | 0.5992254 | 1.6161816 | 0.95 | median | qi |

NA | B - C | contrasts | -0.8609638 | -1.3673337 | -0.3440654 | 0.95 | median | qi |

NA | B - D | contrasts | -0.0088381 | -0.4878745 | 0.4841102 | 0.95 | median | qi |

NA | B - E | contrasts | 1.9508883 | 1.4587024 | 2.4416107 | 0.95 | median | qi |

NA | C - D | contrasts | 0.8464505 | 0.3731628 | 1.3487615 | 0.95 | median | qi |

NA | C - E | contrasts | 2.8086519 | 2.3257780 | 3.3000857 | 0.95 | median | qi |

NA | D - E | contrasts | 1.9641383 | 1.4606453 | 2.4563015 | 0.95 | median | qi |

`MCMCglmm`

Let’s do the same example as above again, this time using
`MCMCglmm::MCMCglmm()`

instead of `rstanarm`

. The
only different when using `MCMCglmm()`

is that to use
`MCMCglmm()`

with `emmeans()`

you must also pass
the original data used to fit the model to the `emmeans()`

call (see `vignette("models", package = "emmeans"))`

for more
information).

First, we’ll fit the model:

```
# 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))
```

Now we can use `emmeans()`

and
`gather_emmeans_draws()`

exactly as we did with
`rstanarm`

, but we need to include a `data`

argument in the `emmeans()`

call: