The posterior R package

Available formats

Because different formats are preferable in different situations, posterior supports multiple formats and easy conversion between them. The currently supported formats are:

• draws_array: An iterations by chains by variables array.
• draws_matrix: A draws (iterations x chains) by variables array.
• draws_df: A draws by variables data frame with addition meta columns .chain, .iteration, .draw.
• draws_list: A list with one sublist per chain. Each sublist is a named list with one vector of iterations per variable.
• draws_rvars: A list of random variable rvar objects, one per variable. See vignette("rvar") for an introduction to this new data type.

These formats are essentially base R object classes and can be used as such. For example, a draws_matrix object is just a matrix with a little more consistency (e.g., no dropping of dimensions with one level when indexing) and additional methods. The exception to this is the draws_rvars format, which contains rvar objects that behave somewhat like arrays but are really a unique data type. See the separate vignette on the rvar and draws_rvars data types for details.

The draws for our example come as a draws_array object with 100 iterations, 4 chains, and 10 variables:

str(eight_schools_array)

##  'draws_array' num [1:100, 1:4, 1:10] 2.01 1.46 5.81 6.85 1.81 ...
##  - attr(*, "dimnames")=List of 3
##   ..$ iteration: chr [1:100] "1" "2" "3" "4" ...
##   ..$ chain    : chr [1:4] "1" "2" "3" "4"
##   ..$ variable : chr [1:10] "mu" "tau" "theta[1]" "theta[2]" ...

Converting between formats

Each of the formats has a method as_draws_<format> (e.g., as_draws_list()) for creating an object of the class from any of the other formats. As a demonstration we can convert the example draws_array to a draws_df, a data frame with additional meta information. To convert to a draws_df we use as_draws_df().

eight_schools_df <- as_draws_df(eight_schools_array)
str(eight_schools_df)

## draws_df [400 × 13] (S3: draws_df/draws/tbl_df/tbl/data.frame)
##  $ mu        : num [1:400] 2.01 1.46 5.81 6.85 1.81 ...
##  $ tau       : num [1:400] 2.77 6.98 9.68 4.79 2.85 ...
##  $ theta[1]  : num [1:400] 3.962 0.124 21.251 14.7 5.96 ...
##  $ theta[2]  : num [1:400] 0.271 -0.069 14.931 8.586 1.156 ...
##  $ theta[3]  : num [1:400] -0.743 0.952 1.829 2.675 3.109 ...
##  $ theta[4]  : num [1:400] 2.1 7.28 1.38 4.39 1.99 ...
##  $ theta[5]  : num [1:400] 0.923 -0.062 0.531 4.758 0.769 ...
##  $ theta[6]  : num [1:400] 1.65 11.26 7.16 8.1 4.66 ...
##  $ theta[7]  : num [1:400] 3.32 9.62 14.8 9.49 1.21 ...
##  $ theta[8]  : num [1:400] 4.85 -8.64 -1.74 5.28 -4.54 ...
##  $ .chain    : int [1:400] 1 1 1 1 1 1 1 1 1 1 ...
##  $ .iteration: int [1:400] 1 2 3 4 5 6 7 8 9 10 ...
##  $ .draw     : int [1:400] 1 2 3 4 5 6 7 8 9 10 ...

print(eight_schools_df)

## # A draws_df: 100 iterations, 4 chains, and 10 variables
##      mu tau theta[1] theta[2] theta[3] theta[4] theta[5] theta[6]
## 1  2.01 2.8     3.96    0.271    -0.74      2.1    0.923      1.7
## 2  1.46 7.0     0.12   -0.069     0.95      7.3   -0.062     11.3
## 3  5.81 9.7    21.25   14.931     1.83      1.4    0.531      7.2
## 4  6.85 4.8    14.70    8.586     2.67      4.4    4.758      8.1
## 5  1.81 2.8     5.96    1.156     3.11      2.0    0.769      4.7
## 6  3.84 4.1     5.76    9.909    -1.00      5.3    5.889     -1.7
## 7  5.47 4.0     4.03    4.151    10.15      6.6    3.741     -2.2
## 8  1.20 1.5    -0.28    1.846     0.47      4.3    1.467      3.3
## 9  0.15 3.9     1.81    0.661     0.86      4.5   -1.025      1.1
## 10 7.17 1.8     6.08    8.102     7.68      5.6    7.106      8.5
## # ... with 390 more draws, and 2 more variables
## # ... hidden reserved variables {'.chain', '.iteration', '.draw'}

Converting regular R objects to draws formats

The example draws already come in a format natively supported by posterior, but we can of course also import the draws from other sources like common base R objects.

x <- matrix(rnorm(50), nrow = 10, ncol = 5)
colnames(x) <- paste0("V", 1:5)
x <- as_draws_matrix(x)
print(x)

## # A draws_matrix: 10 iterations, 1 chains, and 5 variables
##     variable
## draw     V1     V2     V3    V4    V5
##   1   0.983  0.951 -0.029  1.58 -0.56
##   2   1.157  1.049 -0.547  0.14  1.66
##   3  -2.124  0.582  0.921 -2.15 -0.33
##   4  -1.807 -1.366  0.894 -0.57  1.05
##   5  -0.774 -0.237  1.202 -0.16 -0.80
##   6   0.209 -0.020 -1.171 -0.95 -1.89
##   7  -1.166  0.011 -0.110  0.17 -1.24
##   8   0.603 -0.737 -0.337  0.75  0.11
##   9  -0.600 -1.029 -0.565 -1.53  0.27
##   10 -0.072  0.744  0.659 -1.35  0.86

x <- draws_matrix(alpha = rnorm(50), beta = rnorm(50))
print(x)

## # A draws_matrix: 50 iterations, 1 chains, and 2 variables
##     variable
## draw  alpha   beta
##   1  -0.981  0.728
##   2  -0.213 -1.102
##   3  -1.015 -0.831
##   4   2.204 -0.841
##   5  -1.096 -1.138
##   6  -1.055  0.625
##   7  -1.458  1.002
##   8   0.134 -0.163
##   9   0.047 -0.061
##   10 -0.006 -1.304
## # ... with 40 more draws

Subsetting

Subsetting draws objects can be done according to various aspects of the draws (iterations, chains, or variables). The posterior package provides a convenient interface for this purpose via subset_draws(). For example, here is the code to extract the first five iterations of the first two chains of the variable mu.

sub_df <- subset_draws(eight_schools_df, variable = "mu", chain = 1:2, iteration = 1:5)
str(sub_df)

## draws_df [10 × 4] (S3: draws_df/draws/tbl_df/tbl/data.frame)
##  $ mu        : num [1:10] 2.01 1.46 5.81 6.85 1.81 ...
##  $ .chain    : int [1:10] 1 1 1 1 1 2 2 2 2 2
##  $ .iteration: int [1:10] 1 2 3 4 5 1 2 3 4 5
##  $ .draw     : int [1:10] 1 2 3 4 5 6 7 8 9 10

The same call to subset_draws() can be used regardless of the draws format. For example, here is the same code except replacing the draws_df object with the draws_array object.

sub_arr <- subset_draws(eight_schools_array, variable = "mu", chain = 1:2, iteration = 1:5)
str(sub_arr)

##  'draws_array' num [1:5, 1:2, 1] 2.01 1.46 5.81 6.85 1.81 ...
##  - attr(*, "dimnames")=List of 3
##   ..$ iteration: chr [1:5] "1" "2" "3" "4" ...
##   ..$ chain    : chr [1:2] "1" "2"
##   ..$ variable : chr "mu"

We can check that these two calls to subset_draws() (the first with the data frame, the second with the array) produce the same result.

identical(sub_df, as_draws_df(sub_arr))
identical(as_draws_array(sub_df), sub_arr)

## [1] TRUE
## [1] TRUE

It is also possible to use standard R subsetting syntax with draws objects. The following is equivalent to the use of subset_draws() with the array above.

eight_schools_array[1:5, 1:2, "mu"]

## # A draws_array: 5 iterations, 2 chains, and 1 variables
## , , variable = mu
## 
##          chain
## iteration   1    2
##         1 2.0  3.0
##         2 1.5  8.2
##         3 5.8 -1.2
##         4 6.8 10.9
##         5 1.8  9.8

The major difference between how posterior behaves when indexing and how base R behaves is that posterior will not drop dimensions with only one level. That is, even though there is only one variable left after subsetting, the result of the subsetting above is still a draws_array and not a draws_matrix.

Mutating (transformations of variables)

The magic of having obtained draws from the joint posterior (or prior) distribution of a set of variables is that these draws can also be used to obtain draws from any other variable that is a function of the original variables. That is, if we are interested in the posterior distribution of, say, phi = (mu + tau)^2 all we have to do is to perform the transformation for each of the individual draws to obtain draws from the posterior distribution of the transformed variable. This procedure is handled by mutate_variables().

x <- mutate_variables(eight_schools_df, phi = (mu + tau)^2)
x <- subset_draws(x, c("mu", "tau", "phi"))
print(x)

## # A draws_df: 100 iterations, 4 chains, and 3 variables
##      mu tau   phi
## 1  2.01 2.8  22.8
## 2  1.46 7.0  71.2
## 3  5.81 9.7 240.0
## 4  6.85 4.8 135.4
## 5  1.81 2.8  21.7
## 6  3.84 4.1  62.8
## 7  5.47 4.0  88.8
## 8  1.20 1.5   7.1
## 9  0.15 3.9  16.6
## 10 7.17 1.8  79.9
## # ... with 390 more draws
## # ... hidden reserved variables {'.chain', '.iteration', '.draw'}

Renaming

To rename variables use rename_variables(). Here we rename the scalar mu to mean and the vector theta to alpha.

# mu is a scalar, theta is a vector
x <- rename_variables(eight_schools_df, mean = mu, alpha = theta)
variables(x)

##  [1] "mean"     "tau"      "alpha[1]" "alpha[2]" "alpha[3]" "alpha[4]"
##  [7] "alpha[5]" "alpha[6]" "alpha[7]" "alpha[8]"

In the call to rename_variables() above, mu and theta can be quoted or unquoted.

It is also possible to rename individual elements of non-scalar parameters, for example we can rename just the first element of alpha:

x <- rename_variables(x, a1 = `alpha[1]`)
variables(x)

##  [1] "mean"     "tau"      "a1"       "alpha[2]" "alpha[3]" "alpha[4]"
##  [7] "alpha[5]" "alpha[6]" "alpha[7]" "alpha[8]"

Binding

The bind_draws() method can be used to combine draws objects along different dimensions. As an example, suppose we have several different draws_matrix objects:

x1 <- draws_matrix(alpha = rnorm(5), beta = rnorm(5))
x2 <- draws_matrix(alpha = rnorm(5), beta = rnorm(5))
x3 <- draws_matrix(theta = rexp(5))

We can bind x1 and x3 together along the 'variable' dimension to get a single draws_matrix with the variables from both x1 and x3:

x4 <- bind_draws(x1, x3, along = "variable")
print(x4)

## # A draws_matrix: 5 iterations, 1 chains, and 3 variables
##     variable
## draw alpha  beta theta
##    1 -1.57  0.40 1.720
##    2  0.88 -1.11 1.330
##    3  1.08  0.62 0.135
##    4  0.62  1.67 0.996
##    5  0.30  0.92 0.049

Because x1 and x2 have the same variables, we can bind them along the 'draw' dimension to create a single draws_matrix with more draws:

x5 <- bind_draws(x1, x2, along = "draw")
print(x5)

## # A draws_matrix: 10 iterations, 1 chains, and 2 variables
##     variable
## draw  alpha  beta
##   1  -1.574  0.40
##   2   0.876 -1.11
##   3   1.076  0.62
##   4   0.620  1.67
##   5   0.296  0.92
##   6  -0.378 -0.46
##   7   0.015 -0.60
##   8  -0.137  0.74
##   9  -1.298  0.19
##   10 -0.359  0.37

As with all posterior methods, bind_draws() can be used with all draws formats and depending on the format different dimensions are available to bind on. For example, we can bind draws_array objects together by iteration, chain, or variable, but a 2-D draws_matrix with the chains combined can only by bound by draw and variable.

summarise_draws() basic usage

Computing summaries of posterior or prior draws and convergence diagnostics for posterior draws are some of the most common tasks when working with Bayesian models fit using Markov Chain Monte Carlo (MCMC) methods. The posterior package provides a flexible interface for this purpose via summarise_draws() (or summarize_draws()), which can be passed any of the formats supported by the package.

# summarise_draws or summarize_draws
summarise_draws(eight_schools_df)

## # A tibble: 10 × 10
##    variable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 mu         4.2    4.2   3.4   3.6  -0.85   9.4   1.0     558.     322.
##  2 tau        4.2    3.1   3.6   2.9   0.31  11.    1.0     246.     202.
##  3 theta[1]   6.7    6.0   6.3   4.9  -1.2   19.    1.0     400.     254.
##  4 theta[2]   5.3    5.1   4.6   4.3  -2.0   13.    1.0     564.     372.
##  5 theta[3]   3.0    4.0   6.8   4.9 -10.    12.    1.0     312.     205.
##  6 theta[4]   4.9    5.0   4.9   4.5  -3.6   12.    1.0     695.     252.
##  7 theta[5]   3.2    3.7   5.1   4.4  -5.9   11.    1.0     523.     306.
##  8 theta[6]   4.0    4.1   5.2   4.8  -4.3   12.    1.0     548.     205.
##  9 theta[7]   6.5    5.9   5.3   4.5  -1.2   15.    1.0     434.     308.
## 10 theta[8]   4.6    4.6   5.3   4.9  -3.8   12.    1.0     355.     146.

The result is a data frame with one row per variable and one column per summary statistic or convergence diagnostic. The summaries rhat, ess_bulk, and ess_tail are described in Vehtari et al. (2020). We can choose which summaries to compute by passing additional arguments, either functions or names of functions. For instance, if we only wanted the mean and its corresponding Monte Carlo Standard Error (MCSE) we could use either of these options:

# the function mcse_mean is provided by the posterior package
s1 <- summarise_draws(eight_schools_df, "mean", "mcse_mean") 
s2 <- summarise_draws(eight_schools_df, mean, mcse_mean) 
identical(s1, s2)

## [1] TRUE

print(s1)

## # A tibble: 10 × 3
##    variable  mean mcse_mean
##    <chr>    <dbl>     <dbl>
##  1 mu         4.2      0.15
##  2 tau        4.2      0.21
##  3 theta[1]   6.7      0.32
##  4 theta[2]   5.3      0.20
##  5 theta[3]   3.0      0.45
##  6 theta[4]   4.9      0.19
##  7 theta[5]   3.2      0.23
##  8 theta[6]   4.0      0.22
##  9 theta[7]   6.5      0.25
## 10 theta[8]   4.6      0.27

Changing column names

The column names in the output can be changed by providing the functions as name-value pairs, where the name is the name to use in the output and the value is a function name or definition. For example, here we change the names mean and sd to posterior_mean and posterior_sd.

summarise_draws(eight_schools_df, posterior_mean = mean, posterior_sd = sd)

## # A tibble: 10 × 3
##    variable posterior_mean posterior_sd
##    <chr>             <dbl>        <dbl>
##  1 mu                  4.2          3.4
##  2 tau                 4.2          3.6
##  3 theta[1]            6.7          6.3
##  4 theta[2]            5.3          4.6
##  5 theta[3]            3.0          6.8
##  6 theta[4]            4.9          4.9
##  7 theta[5]            3.2          5.1
##  8 theta[6]            4.0          5.2
##  9 theta[7]            6.5          5.3
## 10 theta[8]            4.6          5.3

Using custom functions

For a function to work with summarise_draws(), it needs to take a vector or matrix of numeric values and return a single numeric value or a named vector of numeric values. Additional arguments to the function can be specified in a list passed to the .args argument.

weighted_mean <- function(x, wts) {
  sum(x * wts)/sum(wts)
}
summarise_draws(
  eight_schools_df, 
  weighted_mean, 
  .args = list(wts = rexp(ndraws(eight_schools_df)))
)

## # A tibble: 10 × 2
##    variable weighted_mean
##    <chr>            <dbl>
##  1 mu                 4.2
##  2 tau                4.0
##  3 theta[1]           6.6
##  4 theta[2]           5.3
##  5 theta[3]           3.3
##  6 theta[4]           5.0
##  7 theta[5]           3.3
##  8 theta[6]           4.2
##  9 theta[7]           7.0
## 10 theta[8]           4.9

Specifying functions using lambda-like syntax

It is also possible to specify a summary function using a one-sided formula that follows the conventions supported by rlang::as_function(). For example, the function

function(x) quantile(x, probs = c(0.4, 0.6))

can be simplified to

# for multiple arguments `.x` and `.y` can be used, see ?rlang::as_function
~quantile(., probs = c(0.4, 0.6))

Both can be used with summarise_draws() and produce the same output:

summarise_draws(eight_schools_df, function(x) quantile(x, probs = c(0.4, 0.6)))

## # A tibble: 10 × 3
##    variable `40%` `60%`
##    <chr>    <dbl> <dbl>
##  1 mu         3.4   5.4
##  2 tau        2.5   4.0
##  3 theta[1]   4.9   7.0
##  4 theta[2]   4.3   6.1
##  5 theta[3]   2.5   5.3
##  6 theta[4]   3.8   6.1
##  7 theta[5]   2.7   4.7
##  8 theta[6]   2.9   5.5
##  9 theta[7]   4.8   7.3
## 10 theta[8]   3.5   5.9

summarise_draws(eight_schools_df, ~quantile(.x, probs = c(0.4, 0.6)))

## # A tibble: 10 × 3
##    variable `40%` `60%`
##    <chr>    <dbl> <dbl>
##  1 mu         3.4   5.4
##  2 tau        2.5   4.0
##  3 theta[1]   4.9   7.0
##  4 theta[2]   4.3   6.1
##  5 theta[3]   2.5   5.3
##  6 theta[4]   3.8   6.1
##  7 theta[5]   2.7   4.7
##  8 theta[6]   2.9   5.5
##  9 theta[7]   4.8   7.3
## 10 theta[8]   3.5   5.9

See help("as_function", "rlang") for details on specifying these functions.

Other diagnostics

In addition to the default diagnostic functions used by summarise_draws() (rhat(), ess_bulk(), ess_tail()), posterior also provides additional diagnostics like effective sample sizes and Monte Carlo standard errors for quantiles and standard deviations, an experimental new diagnostic called R*, and others. For a list of available diagnostics and links to their individual help pages see help("diagnostics", "posterior").

If you have suggestions for additional diagnostics that should be implemented in posterior, please open an issue at https://github.com/stan-dev/posterior/issues.