Model comparison (LOO?) from cmdstanr objects


I am attempting to add in some form of model comparison, weights etc toa project that makes use of cmdstanr.

I have used some of these tools before on brms objects, but never in a project in which I’ve written my own stan code.

Are there any good tutorials that cover current best practices? [There seems to be a range of different packages and approaches]

I am under the impression that the first thing i need to do is to make sure I compute log_lik[] in the generated quantities {} block.

This appears to be working. If I fit my model m I can call m$loo() and I obtain this:

Computed from 4000 by 500 log-likelihood matrix
         Estimate   SE
elpd_loo   -375.2 16.6
p_loo         3.1  0.3
looic       750.4 33.2
Monte Carlo SE of elpd_loo is 0.0.

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

I have now fitted a second version of my model, m2 (which is identical, … I’m just trying to test my code!). Reassuringly, this gives very similar output.

How would I go about computing model weights for these two models?

[Once I have a better handle on what I’m doing, I want to eventually compare my non-linear multi-level model with a linear multi-level model. But, baby steps…]


The loo package has a loo_model_weights function that you can use with a list of loo objects, e.g.

loo1 <- m1$loo()
loo2 <- m2$loo()
loo_model_weights(list("Model1" = loo1, "Model2" = loo2))


That was one of the functions that I had tried, but my syntax was all wrong.

Super-helpful. you are a star.

Looks like your Pareto Ks are good now.
In the case if they become problematic, check out diagnostics for PSIS-LOO


At the moment, I’m just testing trivial simulation examples, and haven’t yet scaled up to real, multi-level data.

I expect the Pareto Ks will go downhill!