Where did loo difference estimates go?

Interfaces brms
18

Now I’ve had some to play with the new brms/loo workflow, it might be helpful to others to see what it can look like in practice. Here I’ll compare fit1 and fit2 from @avehtari and @jonah’s vignette Bayesian Stacking and Pseudo-BMA weights using the loo package, which are themselves based on models from McElreath’s Statistical Rethinking, Chapter 6.

First, load the packages and format the data.

library(tidyverse)
library(rethinking)
library(brms)
data(milk)

milk <- 
  milk %>%
  drop_na(ends_with("_s")) %>%
  mutate(neocortex = neocortex.perc / 100)

Fit the two models.

fit1 <- 
  brm(data = milk, family = gaussian,
      kcal.per.g ~ 1,
      seed = 1)

fit2 <- 
  brm(data = milk, family = gaussian,
      kcal.per.g ~ 1 + neocortex,
      seed = 1)

Here we’ll compute the LOO and WAIC for both models and them save the results with the brmfit objects.

fit1 <- add_criterion(fit1, c("waic", "loo"))
fit2 <- add_criterion(fit2, c("waic", "loo"))

You can access the results like this.

fit1$loo

Computed from 4000 by 17 log-likelihood matrix

         Estimate  SE
elpd_loo      4.4 1.9
p_loo         1.3 0.3
looic        -8.8 3.7
------
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.

fit2$waic

Computed from 4000 by 17 log-likelihood matrix

          Estimate  SE
elpd_waic      3.6 1.6
p_waic         1.9 0.3
waic          -7.2 3.2

To compare models by their information criteria, use the loo_compare() function. Since it’ll come in handy in a bit, we’ll first save the results as ic.

ic <- loo_compare(fit1, fit2, criterion = "waic")

Here’s the default output.

print(ic)

     elpd_diff se_diff
fit1  0.0       0.0   
fit2 -0.7       0.6

However, there’s more information available.

print(ic, simplify = FALSE)

     elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic se_waic
fit1  0.0       0.0     4.4       1.9          1.3    0.3      -8.7  3.7   
fit2 -0.7       0.6     3.6       1.6          1.9    0.3      -7.3  3.2

Even with all that output, you’ll notice that the difference score and its standard error are in the elpd metric, rather than the conventional information criteria metric. Here’s a quick way to do the conversion.

cbind(waic_diff = ic[, 1] * -2,
      se        = ic[, 2] *  2)

     waic_diff       se
fit1  0.000000 0.000000
fit2  1.416582 1.153017
1 Like