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
```