I’ve been using the point estimates from `fixef(mod)`

.

So in order to manually obtain the same fitted values as `fitted.brmsfit()`

, I would take the 4000 posterior samples of logit coefficients, use them to construct 4000 matrices of fitted values (transformed back into probabilities), and finally calculate the average over the 4000 matrices for every observation?