Hi Angelo, sorry for the slow reply!

I think it’s easiest to see why the mean and median behave differently if we focus on a single observation. For example, for the model I posted above here’s a plot of the posterior distribution of \lambda_i vs the predictive distribution of y^{\rm rep}_i for i=1, i.e., the first column from `posterior_epred()`

vs the first column from `posterior_predict()`

. Because of the shapes of the distributions (the latter including the impact of the dispersion term), they have the same mean but different medians:

```
fit <- stan_glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = MASS::quine)
lambda_1 <- posterior_epred(fit)[, 1]
yrep_1 <- posterior_predict(fit)[, 1]
make_hist <- function(x) {
hist(x, main = deparse(substitute(x)))
abline(v = mean(x), col = "red")
abline(v = median(x), col = "blue")
legend("right",
legend = c(paste("mean =", round(mean(x), 1)),
paste("median =", round(median(x), 1))),
col = c("red", "blue"),
lty = 1)
}
make_hist(lambda_1)
make_hist(yrep_1)
```

In both cases the mean is about 18 while the medians differ.

Does that help explain the reason for what you’re seeing?