How to obtain random draws from grouped residuals?

#1

Suppose I fit the data with a grouped residuals in brms:

fm <- brm(bf(…, sigma ~ 0+factor), data=…))

I can get the grouped \sigma^2 for each level of the factor from

summary(fm)

However, is there a way to obtain random draws for the grouped residuals similar to what ranef(fm) does?

#2

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?

#3

Indeed they follow the same sequence as the input data.

1 Like
#4

Thanks for the confirmation, Paul!

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?

#5

I don’t understand your question I am afraid.

#6

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?

#7

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.

1 Like