Visualizing the population distributions for the random intercepts and/or slopes?

General
1

Hello. Not sure where this fits, and this may be an easy question. Suppose I run a hierarchical model with random intercepts. I want to generate a plot of the population distribution, accounting for parameter uncertainty on mu and sigma, the population parameters. How do I do that?
I can extract all the posterior draws using the posterior package in a data frame. Every pair of mu and sigma posterior draws would give me a population distribution. Is there a way to plot all the population distributions in a graph? Does the bayesplot package have a command for this? Or do I need a separate command in R? Any feedback is appreciated.

2

Probably a density plot similar to the style of plots in pp_check() would be effective. You could have a highly opaque density line for the posterior median, and maybe overlay the density from 100 random samples or something. Perhaps Draw a function as a continuous curve — geom_function • ggplot2 would be a simple way to do it?

3

I posted the same on Reddit and got a response from a bayerplot dev. I am sharing it here and I will test this out later. But I believe this can achieve what I asked:

Bayesplot dev here 👋, but bayesplot will be no help for you. ggdist is probably your best bet.

geom_function() does not let you treat function arguments as plotting aesthetics, so you can only plot one function’s curve at a time.

That said, we can fake what geom_function() does by computing the densities ourselves beforehand and using geom_line().

library(dplyr)
library(ggplot2)

mus <- rnorm(10, 4, 1)
sigmas <- abs(rnorm(10, .5, .1))

x <- data.frame(
  draw = seq_along(mus),
  mu = mus,
  sigma = sigmas
)

# Find the overall x limits of the plot
x_lim <- range(c(
  qnorm(p = .0001, x$mu, x$sigma), 
  qnorm(p = .9999, x$mu, x$sigma)
))

# Draw the curves' densities at these x values
x_grid <- seq(x_lim[1], x_lim[2], length.out = 1000)


# Compute the densities on the x-grid for each curve
x <- x |> 
  group_by(draw, mu, sigma) |> 
  summarise(
    x = x_grid,
    y = dnorm(x, mu, sigma)
  ) |> 
  ungroup()

ggplot(x) +
  geom_line(
    aes(x = x, y = y, group = draw)
  )