Recently, I joined the development team for
bayesplot, an R package by the Stan team for
plotting Bayesian models. Because visualizing Bayesian models in ggplot2 is a
recurring
topic
here, it was a natural fit. So from time to
time, I’ll post about some programming techniques and new features we develop in
the bayesplot package.
For this post, I describe one of strategies I have been
using to clean up and simplify some of the plotting code in the package: avoid
if-branches by sometimes plotting nothing.
Warm-up example
Let’s start with a non-plotting example. We consider a function
that takes a list of counts of people and returns the total number of people.
That sum, however, is controlled by some of the function’s arguments:
head_count [1] 1 2 3
# data-frame binding
rbind(head(iris), list())
#__ Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#__ 1 5.1 3.5 1.4 0.2 setosa
#__ 2 4.9 3.0 1.4 0.2 setosa
#__ 3 4.7 3.2 1.3 0.2 setosa
#__ 4 4.6 3.1 1.5 0.2 setosa
#__ 5 5.0 3.6 1.4 0.2 setosa
#__ 6 5.4 3.9 1.7 0.4 setosa
# magrittr pipes
f % f()
#__ [1] 1 2
Below, I demonstrate two examples of how I applied a similar strategy in
bayesplot. In ggplot2, we assemble plots by adding layers and options to a plot
by stringing together statements with +. To handle conditional plotting
elements and construct the plot all in one pass, I took advantage of ways to
adding nothing to a ggplot2 plot.1
Posterior interval plots
Here are some posterior samples for a Bayesian model, in this case the “eight schools” meta-analysis example. We want to create a plot that
shows the 90% and 50% uncertainty intervals for each parameter.
We first use a helper function to clean up the data to prepare it for plotting.
(At the time of writing, mcmc_intervals_data() is not yet in the CRAN version
of bayesplot.)
library(dplyr, warn.conflicts = FALSE)
library(bayesplot)
#__ This is bayesplot version 1.4.0.9000
#__ - Plotting theme set to bayesplot::theme_default()
#__ - Online documentation at mc-stan.org/bayesplot
library(ggplot2)
theme_set(theme_grey())
interval_data % mcmc_intervals_data(pars = &quot;mu&quot;, regex_pars = &quot;theta&quot;) %__% mutate_if(is.numeric, round, 2) interval_data
#__ # A tibble: 9 x 9
#__ parameter outer_width inner_width point_est ll l m h
#__ #__ 1 mu 0.9 0.5 median -0.66 4.74 8.11 11.54
#__ 2 theta[1] 0.9 0.5 median -0.40 6.63 11.24 17.05
#__ 3 theta[2] 0.9 0.5 median -3.16 3.63 7.95 12.20
#__ 4 theta[3] 0.9 0.5 median -8.86 1.29 6.33 11.02
#__ 5 theta[4] 0.9 0.5 median -3.69 3.47 7.84 12.01
#__ 6 theta[5] 0.9 0.5 median -6.84 0.47 5.10 8.99
#__ 7 theta[6] 0.9 0.5 median -6.53 1.76 6.27 10.42
#__ 8 theta[7] 0.9 0.5 median 0.33 6.67 10.84 15.71
#__ 9 theta[8] 0.9 0.5 median -4.97 3.79 8.53 13.50
#__ # ... with 1 more variables: hh Below is the basic plot we want to
create. (Normally, in bayesplot, we would use
mcmc_intervals()
to make this plot, but for the sake of illustration, we will write a version of
it from scratch.
ggplot(interval_data) + aes_(y = ~ parameter, yend = ~ parameter) +
geom_segment(aes_(x = ~ ll, xend = ~ hh), size = 1) + geom_segment(aes_(x = ~ l, xend = ~ h), size = 2) + scale_y_discrete(limits = rev(interval_data$parameter)) + labs(x = NULL, y = NULL)
Wait, what’s going on with the tildes? We use aes_(y = ~ parameter)
because in an R package, writing aes(y = parameter) would raise an undefined
global variable warning :warning: during CRAN package checks. The code looks like it uses
a variable called parameter, but a variable with that name has not been defined in
the code yet. The formula form ~ parameter “quotes” the variable name, so
it doesn’t appear as a global variable.
We are going to wrap this code in an R function along with two annotation
options:
whether to draw point estimates (medians) over the intervals
whether to draw a vertical reference line at x = 0
Here is how we might write the function with branching code.
plot_intervals