predictive_error(fm) seems to generate random draws for the residuals. However, the output matrix rows do not show the grouped structure, nor does the output matrix contain dimension names. Do the random draws of the residuals follow the same sequence as the input data?

I have a conceptual question related to this. The random samples for the population effects, grouped effects and the residuals are separately drawn with fixef(), ranef(), and predictive_error(). It seems to me that each random sample should be constrained by the model together with all the effect components. So is it conceptually appropriate to separately draw those samples without the exact constraint?

Suppose I draw 1000 samples from the posterior distribution for b, \theta_j and \epsilon_{ij}. I assume that fixef(), ranef() and predictive_error() in brms generate random draws for the parameter b, \theta_j and \epsilon_{ij} independently of each other. If this is the case, the combination of each sample (e.g., draw #1) for b, \theta_j and \epsilon_{ij} would not exactly meet the equation above, would it? In other words, would it be better to get the residuals directly from

\epsilon_{ij}=y_{ij} - b - \theta_j

instead of predictive_error()? Am I conceptually missing something too off-base?

brms does not think in response of multiple indices (unless you use a multivariate model, which you don’t). That is, for brms, ij does not exist in the results of predictive_error just a single index over all observations from which to work back your ij indices. Still this is waay easier and much less error prone than doing stuff by hand as you suggest.