Thanks very much, Andrew for your help. I’m getting exactly the plots I want and expect, but would like to display them in one plot. The two plotting lines are

These plot two types of prior distributions, with one being reasonably correct domain knowledge as I am trying to display prior predictive checking. The two plots are attached. I’ve been trying a number of methods to get them on the same plot, to no avail.

Do you mean you want to lay out Rplot1 and Rplot2 next to each other (two separate plots layed out together on the same page) or that you want to combine Rplot1 and Rplot2 into a single panel with two sets of prior predictive densities?

Thanks for your response. I would prefer a single panel with two sets of prior predictive densities. If that’s not easily done, I can put them side by side with LaTeX.

The plots returned by bayesplot are simply ggplot objects, so you can use any of the existing frameworks for combining/arranging plots (such as cowplot, gridExtra, and egg)

To lay out the plots side by side, one option is the patchwork package, which has a sort of “grammar” of plot layouts. For example, the following code will lay out the two plots side by side:

To plot posterior predictions from both models in a single panel takes more work (though there may be some convenience functions I’m not aware of). For example, if you have two models fit with brm, let’s call them fit1 and fit2 (these could be from two different prior predictive checks), then you can do this:

library(tidyverse)
# Get posterior predictions
fit1.post = posterior_predict(fit1)
fit2.post = posterior_predict(fit2)
# Get 10 draws from each posterior. Convert them to long data frames
# and stack them into a single data frame. The "model" column is a
# grouping variable marking whether the posterior came from fit1 or fit2.
bind_rows(
fit1.post[sample(1:nrow(fit1.post), 10), ] %>%
t %>%
as.data.frame() %>%
pivot_longer(everything()) %>%
mutate(model="fit1"),
fit2.post[sample(1:nrow(fit2.post), 10), ] %>%
t %>%
as.data.frame() %>%
pivot_longer(everything()) %>%
mutate(model="fit2")
) %>%
# Plot densities of both posterior predictions together in the same plot
ggplot(aes(value, group=interaction(model, name), colour=model)) +
geom_density(alpha=0.3, size=0.2) +
scale_y_continuous(expand=expansion(c(0,0.02))) +
scale_x_continuous(expand=expansion(c(0.01,0.0))) +
theme_classic()