Use a simple model as an example,

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

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?