Applying loo to multi-variate measurement error model in BRMS

Modeling
,
1

I’m having a strange problem using loo with BRMS: Consider estimating a multi-variate measurement error model (two response variables, both measured with error) in BRMS, something like:

m1 ← bf(y1 | se(se1) ~ pterms + (1 | id | group))
m2 ← bf(y2 | se(se2) ~ pterms + (1 | id | group))

When estimating this model I also allow for residual correlations, that is the formula is m1 + m2 + set_rescor(TRUE)). The model estimates fine, with good convergence diagnostics etc.

However, when I try and calculate PSIS-LOO, I am told 100% of my N = 1180 observations are problematic and the model needs to be re-fit N times.

If I instead estimate the non-measurement error version of the same model, as below, everything works like normal (NB: I have around 7 problematic observations).

m1 ← bf(y1 ~ pterms + (1 | id | group))
m2 ← bf(y2 ~ pterms + (1 | id | group))

I wondered if this had something to do with saving latent variables, although save_pars(latent = TRUE, all = TRUE) doesn’t seem to help.

Thanks in advance for any ideas / tips.

2

@paul.buerkner probably can tell what is the difference in the measurement error model that would make the PSIS more likely to fail (Paul may be on vacation, and thus it may take longer time than usual to get his comment)

1 Like
3

Good question. I don’t really know. Can you provide a minimal reproducible example of this behavior?

1 Like
4

0 Model data.rdata (106.2 KB)

Hi @paul.buerkner (and @avehtari),

Thanks both for your responses. Let me try to provide my first REPREX (apologies in advance if I don’t get the balance quite right or if this should instead be shifted to the OP).

The data I’m using is uploaded with this reply and contains 1,180 observations. Using this data, I estimate the following model (NB: The model below is simpler and runs faster than that in the OP, as it doesn’t (1) model correlations between the group effects and (2) I set “set_rescor(FALSE)”. Nonetheless, it still takes about an hour so you may prefer to work with a sub-sample of the data).

form.1a <- bf(ctYfe_wp | se(se_ctYfe_wp, sigma = TRUE) ~ log(effden) + (1 | year) + (1 | area), family = student)
form.1b <- bf(ctRfe | se(se_ctRfe, sigma = TRUE) ~ log(effden) + (1 | year) + (1 | area), family = student)
form.1 <- form.1a + form.1b + set_rescor(FALSE)

model <- brm(form.1, data = data.brms, seed = 2317,
           iter = 3000, warmup = 1000, cores = 3, chains = 3,
           control = list(max_treedepth = 15, adapt_delta = 0.95),
           save_pars = save_pars(all = TRUE))

loo(model)

The final command returns the following warning:

Found 1180 observations with a pareto_k > 0.7 in model 'model'. It is recommended to set 'moment_match = TRUE' in order to perform moment matching for problematic observations.

When I estimate the same model without measurement error then it works fine with loo, as below, which suggests no observations with pareto_k > 0.7.

loo(model.2)
Computed from 6000 by 1180 log-likelihood matrix

------
Monte Carlo SE of elpd_loo is 0.3.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     1178  99.8%   608       
 (0.5, 0.7]   (ok)          2   0.2%   358       
   (0.7, 1]   (bad)         0   0.0%   <NA>      
   (1, Inf)   (very bad)    0   0.0%   <NA>      

All Pareto k estimates are ok (k < 0.7).
See help('pareto-k-diagnostic') for details
5

Thank you! Can you perhaps even send me the fitted model already? Since the problem happens in the post-processing and the model fits for some time, this might be faster for me, in this specific case (usually, the full reprex as you provided is preferred of course).

1 Like
6

Hi @paul.buerkner – of course, the model object seems to be too large (40MB) to upload to the Stan forums, although you should be able to access it via this link:

Let me know if you have any issues accessing that file.

P.s. I also want to thank you and let you know that (1) I love brms with the passion of a thousand burning Bayesian suns and (2) really appreciate the effort you (and others, like Aki) go to help random people (like myself) on the Stan forums.

7

Thanks . This was a bug in brms that should now be fixed in the github brms version.

2 Likes
8

Thanks, I’m very grateful.

1 Like