Good question, with a few approaches depending on who you ask.
Let's set the stage for clarity. Suppose you are fitting a model that looks like
$$ y_{i,j} = \alpha_{j} + \beta_{j}x_i$$
Here, $\alpha_j$ is the intercept and $\beta_j$ is the slope for subject $j$ in our data. If it helps, you could think about this as the linear predictor on the appropriate scale for whatever GLM you're using. If I understand your question, you can compute credible intervals for the $\beta_j$ easily, but want a credible interval for the $\beta$ of a yet-to-be-seen subject. In essence, a prediction interval for that $\beta$.
If that is correct, here are a few ways of doing this.
Marginalizing
If the $\beta_j$ are modelled hierarchy using a centred parameterization, then each of the slopes can be written as
$$ \beta_j = \mu + z_j\sigma $$
where $\mu, \sigma$ are the population level mean and standard deviation. For each subject $j$, we can do inference on the $z_j$ to estimate the $\beta_j$. For a new subject, there is no $z_j$ to estimate. Hence, we can integrate the $z_j$ out to get a population level distribution of $\beta_j$.
To do this, you would:
- Draw a $\mu$ and $\sigma$ from their posteriors
- Draw a $z$ from the prior on $z$ (let's assume $p(z)$ is standard normal)
- Compute $\beta = \mu + z\sigma$
- Repeat $N$ times
The uncertainty in $z$ is integrated out in this appraoch.
Let's demonstrate how this is done using brms (which takes care of writing Stan code so I don't have to). I'll use the sleepstudy data to showcase the approach. To follow along, you'll also need the tidybayes library.
library(tidyverse)
library(lme4)
library(brms)
library(tidybayes)
data("sleepstudy")
fit <- brm(Reaction ~ Days + (Days + 1 | Subject),
data = sleepstudy,
backend = 'cmdstan')
slope = fit %>%
# Samples from the psoterior for population level slope and standard deviation
spread_draws(b_Days, sd_Subject__Days) %>%
mutate(
# Draw from the prior of z
z = rnorm(4000, 0, 1),
# COmpute the slopes
slope = b_Days + sd_Subject__Days)
Whatever brms does
brms has a few approaches to doing this, one (I think) of which is what I describe above. See the documentation here.
The default for brms is interesting, but doesn't provide a rationale for why this should be the default
If "uncertainty" (default), each posterior sample for a new level is drawn from the posterior draws of a randomly chosen existing level. Each posterior sample for a new level may be drawn from a different existing level such that the resulting set of new posterior draws represents the variation across existing levels
