Model stacking and LOO (brms models)

[jrrae]

@avehtari

Based on what you posted in your first post, it seems they agree quite well. m5 is best, m4 the second best, m3 the third best, and the uncertainties in pairwise comparisons are such that comparisons agree with the weights.

jrrae:

That is, can I gain any useful insight from knowing that, for example, model 5 is given a majority of the weight by the Pseudo-BMA+ (and stacking) method but only represents a small improvement (about 1 SE) when I compare it to model 3 when using LOO?

When you look at compare and about 1 SE improvement you are considering only whether m3-m5 (btw I hate IC scale multiplication by -2, ping @paul.buerkner) is positive, while Pseudo-BMA is taking into account also the magnitude of the difference. Thus it’s not easy to calculate weights in your head from the pairwise differences and SEs for multiple models. In both cases you should consider these comparisons to have lot of uncertainty and only trust differences of several SEs or weights 0.000.

I understand my error now. In my head I was thinking that a weight of .60 to .70 felt like much stronger support for a model than a 1 SE LOO difference. I now see that this conceptualization is incorrect.

@avehtari

Model 5 includes model 3, and thus an equivalent model for the weighted average of model 5 and 3 could be formed by increasing the coefficient for x2 in model 5. Two possibilities why model 5 is not doing that already: 1) prior for coefficients is narrow, 2) model 5 has also additive term x1, which is not that useful alone based on predictive performance of model 2 (but seems to be helpful in the interaction) and if the prior for all coefficients is vague then coefficient for x1 is fitting to the noise in model 5, but which is not happening in the model 3. So we might say that averaginf model 5 and 3 corresponds to adding more mass to zero coefficient for covariate x1 and the interaction term x1:x2.

This is super helpful. I understand now.

@avehtari

Show your priors, and tell how much data you have, and I can tell more. It would also help a lot to see loo output for model 5 alone (as for rstan and rstanarm models with p_loo and Pareto diagnostics, I’m not certain what loo prints for brms models).

I only set priors on regression intercepts and and coefficients (both N(0,5)). I scaled my regression inputs by 2 sd’s (to facilitate comparisons between the coefficients since x1 is dichotomous and x2 is continuous). For my 5 outcome measures, they are logit transformed probabilities (doing multivariate beta regression isn’t possible or I would have done it) that haven’t otherwise been transformed. I have a sample size of 4,524.

Here’s the loo output for model 5:

Computed from 4000 by 4524 log-likelihood matrix

         Estimate    SE
elpd_loo -13563.7 187.1
p_loo        45.2   2.0
looic     27127.4 374.3
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.