Cannot do save_psis and reloo in the same loo() call?

My perusal of this thread as well as of Gelman, Hill & Vehtari (2020: 174-8) has led me to believe that the LOO objects returned by loo() can be used to extract a PSIS-LOO fitted value for every observation in the dataset. If true, those fitted values would enable one to calculate LOO estimates of any model-performance metric of one’s choosing.

However, the problem I’m running into is that many of my models require at least one refitting in order to obtain a LOO object. And when I run loo(mymodel, reloo = TRUE, save_psis = TRUE) I get the following error:

Error: passing unknown arguments: save_psis.

Is it truly impossible to both deal with the problematic observations through refitting and save the psis objects in the same call? My models are categorical ones with some 30 covariates and 2 group-level terms.

EDIT: The problem is occurring with both brms version 2.14.0 and 2.14.4, the respective rstan versions being 2.21.2 both. StanHeaders versions are 2.21.0-6 and 2.21.0-7 respectively.

I assume you are using brms and changed the category to interfaces with tag brms. Each interface (brms, rstanarm, rstan) have their own loo method, and it’s more likely that brms developer @paul.buerkner can help with this. You could also try moment_match=TRUE if that would be sufficient instead of need to refit.

Thanks for changing the category. I’m happy to hear it’s a brms problem since that implies that this should be possible.

I’ve never had success with moment_match = TRUE. While I can’t go reproduce the error message right now (R is busy refitting a model to get a loo object), I always just assumed that these models are too complex for that kind of shortcuts to work.

If you can share your model, I’d be interested to know what would such a model that moment matching doesn’t help.

It would be useful also to report the brms and loo package versions so it’s easier for @paul.buerkner to answer.

I have updated my first post with the version numbers of brms, rstan and StanHeaders.

@avehtari I will share the model if all that typing becomes necessary (there are a quazillion category-specific priors), but I’m no longer so sure this is a complexity issue. You see, I now went and tried moment_match = TRUE, and the error was as follows:

Error in object@.MISC$stan_fit_instance$unconstrain_pars(pars) : **
** Exception: Variable Intercept_muIndicative missing (in ‘modelcb46d34638c_e5a474e5dd35842a83ebbf6096cc107c’ at line 83)

Error: Moment matching failed. Perhaps you did not set ‘save_all_pars’ to TRUE when fitting your model?

Does this mean that moment matching will likely work as long as one remembers to specify save_pars = save_pars(all = TRUE) when first fitting the model? That could save a lot of time in the future – unless saving all pars significantly slows down the sampling.

It doesn’t slow sampling but makes the loo object much bigger as in addition log_lik values the all posterior draws are saved in the loo object, too.

@avehtari: Sorry. It turns out I was simply mistaken about moment matching. I tried it on two big models today and it worked both times (I got a warning about some pareto-k’s being “slightly” high, which I assume is not the end of the world). I think the reason I had dismissed it earlier was because I wasn’t keen to start re-fitting all these huge models with save_pars(all = TRUE) just to get loo objects. But the fact that it actually works is very good news since it now appears that I can compute custom measures of model quality using the loo fit.

Perhaps it also makes the model object itself much bigger? I do most model fitting with my laptop because it has a fast CPU and works silently overnight, but I have to switch to my desktop machine in order to obtain LOO fits because with only 8GB of memory, the laptop can’t even do posterior_epred() to these models without running into the memory limit.

In summary: moment_match seems to be a viable workaround, but it would surely still be better if reloo = TRUE, save_psis = TRUE was an available option for brms models.

No problem, great that it works!

@paul.buerkner would know, but he seems to be busy.

Be careful to use only one core in laptop, but yes 8GB can be too little even when not running parallely anything.

It’s not just a workaround, it should be also faster.

Would you like to create an issue (assuming you were using the latest versions and it’s not already in more recent code)

Can you please file an issue on GitHub - paul-buerkner/brms: brms R package for Bayesian generalized multivariate non-linear multilevel models using Stan in which you provide a reproducible example? I will try to include the fix in the new brms version to be released soon.

Yes, saving all parameters may make the model much bigger indeed, depending on what kind of model you are fitting exactly.

Alright, an issue has been filed along with what should be a reproducible example using real data.

Do let me know if there’ s a problem with the example.

A word of warning: the model takes over 90 minutes to fit on my laptop (which was high-end in 2018). I suppose I could reproduce the problem in a smaller, mock dataset after enough trial and error, but I figured that a real data example would be in a certain sense more valuable.

As we await the fix to brms, I was wondering what the best strategy is for avoiding “problematic” observations?

It turns out that save_pars(all = TRUE) is no insurance against high pareto-k values. And when those high pareto-k values happen, the only available option is reloo = TRUE which in turn cannot save the PSIS object. The implication is that I cannot compute LOO fits unless I avoid high pareto-k values somehow.

Does increasing the number of simulations reduce the risk of high pareto-k’s? What about increasing adapt_delta? Or some other parameter tweak?

Wait a minute. Isn’t moment_match = TRUE intended as a remedy for high pareto-k’s, i.e. a faster alternative to reloo = TRUE? What I am experiencing with certain models is that when I use save_pars(all = TRUE) and try to extract a loo object with moment_match = TRUE, I still get a warning about high pareto_k’s, which suggests that moment matching has not successfully dealt with them. And if it cannot successfully deal with high pareto_k’s, why use save_pars(all = TRUE) in the first place? It makes the model objects huge, with no apparent benefit from the standpoint of loo(). What am I missing?

Increasing the number of posterior draws allows estimation of Pareto-k more accurately and there can then be less false positives but on the other hand there can then be also more true positives depending on the actual tail shape of the density ratios.



Moment matching works in some cases, but not in all cases. Read more and see the examples in Implicitly adaptive importance sampling.

1 Like

It appears @paul.buerkner has fixed the problem. After updating brms with devtools::install_github("paul-buerkner/brms"), I have now managed to successfully extract my first loofit for a model with high pareto-k’s using loo(mymodel, reloo = TRUE, save_psis = TRUE).

Thanks also to @avehtari and colleagues for continuing to produce research aimed at making these formidable computational tasks easier. Hopefully in the near future, reloo will no longer be needed. But for now, it’s a life-saving feature.

1 Like