Consider a linear, random intercept model:
$$ y_{ij} \sim \mathcal{N}(\mu_i, \sigma_e^2) $$
where
$$ \mu_i \sim \mathcal{N}(\gamma_0, \sigma_{mu}^2) $$
Estimated in a Bayesian framework using MCMC sampling, I would expect that extracting fitted values for each participant in sample for each MCMC draw of the posterior would give me:
$$ \mu_{i,draw} $$
If I were to average across the $i$ participants, I would get samples from a posterior distribution like:
$$ \bar{\mu}_{draw} $$
which I would expect to have the same value and same uncertainty as $\gamma_0$.
In most software, the model is separated as:
$$ \mu_i = \gamma_0 + u_i $$
where
$$ \gamma_0 \sim \mathcal{N}(0, \sigma_{\mu}^2) $$
When I actually try this, my observation is that the posterior draws for $\gamma_0$ and the average of the fitted values $\mu_i$ have very similar average values, but $\gamma_0$ consistently has a wider uncertainty interval than the uncertainty interval for the average of $\mu_i$, the fitted values for each participant.
Now why I care about this. In the linear model case, I could just use the parameter estimates for $\gamma_0$. However, when fitting say a logistic mixed effects model, the estimate for $\gamma_0$ is on the log odds scale and back transformed to a probability is not the same as getting the averaged probability taking into account the random intercept.
My best guess is that in taking the fitted values:
$$ \hat{\mu_i} = \gamma_0 + u_i $$
and then averaging across all participants and then summarizing that posterior distribution is somehow missing a source of variance, but I've been over this a dozen times and tried to look for articles or books discussing it and have come up empty, so am hoping for some insight here.
R
code if helpful, but it seems to be software agnostic as I tried it through a couple different packages and ended up in the same place. $\endgroup$