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.

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.

@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.

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.

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.