How to calculate Bayesian marginal credible interval? I have a Bayesian model that I've fit using Stan, and I'm trying to figure out the best way to calculate the correct credible interval that I am interested in.
The model is a hierarchical GLM with a random intercept for each of $n$ subjects, so I have $k$ samples from the posterior distribution of the linear predictor.
So, I can get the mean and credible interval of the posterior distribution for each individual fairly easily. To construct a credible interval for each individual, that is easy. I just get either the HDPI or use the quantiles and get $n$ different means with credible intervals. But if I want to get the marginal mean and credible interval for the entire population, I'm not sure what to do.
I think I can use the "pooled" data to calculate a CI (e.g., calculate the mean and quantiles of all the samples for each individual put together), but I am not sure if this is the correct approach. I'm mostly using the book "Statistical Rethinking" to develop my model, but the book does not include any guidance on how to get this CI.
Can anyone provide any guidance on how to construct the correct CI in this case?

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

