I’ve got a categorical model for 4 response classes, with varying intercepts for the 3 logistic sub-regressions. I’d like to plot the posteriors of the three group-level SDs, either as an mcmc_areas() or as mcmc_hist(), with the exponential(2) prior curve on top of each posterior for reference. Is this possible? The only thing I’ve been able to manage is drawing one single exponential(2) curve on top of the entire plot (see screenshot) rather than one scaled exponential curve for each posterior:
I still haven’t found a simple solution, but I managed to jury-rig it by constructing a data frame with the exponential(2) prior densities for relevant parameter values, then adding three geom_line() calls to add these density curves on top of the mcmc_areas(). Given that there are three PDFs in one plot, these constructed exponential densities had to all be divided by 3 and assigned y intercepts equal to 1, 2, and 3.
Unfortunately, I now realize that my ad-hoc approach is wrong. In actual fact, the second posterior has almost twice as high of a peak as the other two, and this is largely masked by the fact that bayespot::mcmc_areas scales the posteriors to have the same total area. The prior curves are therefore mis-scaled. There’s no option to tell mcmc_areas() to stop scaling things (although you can change how it’s done). I guess I could scale down those prior curves accordingly, but the result would look silly. The true height of the tallest posterior should be accurately reflected in the plot.
I’ll update the topic once I have a better solution.
To expand on what @mike-lawrence said, you can use two layers with ggdist combined with scale_thickness_shared() to ensure that the thickness aesthetic (which is what is used to display the densities) has the same scaling in the layers for the prior and the posterior. Something like this:
library(ggplot2)
library(ggdist) # for stat_halfeye, stat_slab, ...
library(distributional) # for dist_exponential
prior = data.frame(
var = c("a", "b", "c"),
value = dist_exponential(rep(1, 3))
)
posterior = data.frame(
var = c("a","b","c"),
value = rgamma(12000, 1:3, 1:3)
)
posterior |>
ggplot() +
stat_halfeye(aes(y = var, x = value)) +
stat_slab(aes(y = var, xdist = value), fill = NA, color = "black", linetype = "22", data = prior) +
scale_thickness_shared()
(sidebar, ggdist recently got a new default density estimator, which should do a decent job on these lower-bounded densities, which you can see on the example for “a” above…)
Not wanting to wait idly for help, I spent most of yesterday doing more jury-rigging (this was before seeing @mjskay 's response.) It involves the use of geom_density(), geom_line()and facet_wrap(ncol = 3, scales = "fixed"), preceded by manual addition of the desired range of prior values and corresponding prior densities to the data. The vertical bars for both the prior and posterior means had to be added manually as well, using geom_segment(). Here’s how it looks:
After seeing mjskay’s post I applied his approach as best I could. It definitely accomplishes something similar with much less code. The only thing I felt the need to add manually were the vertical bars for the prior means (I don’t want another “eyeball” on the x-axes). Here’s how it looks:
I wonder why the sampling density curves are bumpier in the ggdist version. It seems to me that this bumpiness risks drowning out / obscuring some interesting little peaks in the posterior.
Or was the bumpiness there all along, such that it is actually geom_density() that needlessly levels out detail?
I would t put much credence into these unless you are doing a very high number of sampling iterations; otherwise they’re likely just reflections of Monte Carlo sampling error
Nice, I like it! Curious if you used stat_spike(at = mean) for those? That’s a fairly recent addition to ggdist designed for this kind of thing.
As @mike-lawrence says, this is just noise. The difference stems from different bandwidth estimators: geom_density() uses bw.nrd0() to pick the bandwidth, which is the default bandwidth estimator of stats::density(). In my experience, that estimator does well for distributions that are close to Gaussian, but can perform poorly on distributions with other shapes, particularly ones where certain regions of the distribution really need a smaller bandwidth to correctly estimate their density (e.g., multimodal distributions where the region around one mode is much narrower than the others, or bounded distributions with peaks near the boundary).
Indeed, the documentation for stats::density() recommends against using bw.nrd0() as a default; rather, it defaults to that for historical reasons. As part of the update of the default density estimator in the most recent version of ggdist, I changed the default bandwidth estimator to bw.SJ(method = "dpi"), which is the Sheather-Jones direct plug-in bandwidth estimator. I have found it handles distributions where some regions require a smaller bandwidth better than bw.nrd0(). The tradeoff is that, because of the smaller bandwidth, you can get more noise in flatter regions of the density curve, as you see in the examples above. I think that the improved estimation of high-density regions is worth the slight increase in noise in flatter regions, which is relatively easy for our eyes to “smooth over”.
I didn’t know about stat_spike() so had to use three separate geom_segment() calls instead. But thanks for the tip! I’m trying to adhere to the general look of bayesplot in posterior graphs throughout my analysis chapter, and combining stat_slab() with stat_spike() for the posterior enables me to do that more accurately than before!
