# load only the turnout data frame and create rescaled columns
turnout <- ZeligData::turnout |>
mutate(across(age:income, arm::rescale, .names = "rs_{.col}"))
# linear in age
f_linear <- vote ~ rs_age + rs_educate + rs_income + race
# cubic in age
f_cubic <- vote ~ rs_age + I(rs_age^2) + I(rs_age^3) +
rs_educate + rs_income + race23 LOOIC
In Information Criteria, we compared Poisson and NB for Holland’s Santiago data with AIC and BIC. However, when we have posterior simulations, the AIC and BIC work poorly because we do not have a single point estimate. Instead, we have a distribution. When we have a sample from the posterior, we can instead use the LOOIC.
The AIC also works poorly with posterior simulations for two other reasons. First, plugging in a point estimate (e.g., the posterior mean) ignores the uncertainty in the posterior, while the LOOIC averages over all the draws. Second, the AIC’s penalty counts parameters, but a prior that pulls the estimates toward zero makes a model less flexible than its count of parameters suggests. With a strong prior or a hierarchical model, \(k\) overstates the model’s effective number of parameters.
23.1 Example: Linear or cubic age?
When we have posterior simulations via Stan, we can use the LOOIC as an alternative to the AIC or BIC.
First, the data and formulas.
Next fit the models.
fit_linear_brm <- brm(f_linear, data = turnout, family = bernoulli,
backend = "cmdstanr", chains = 4, cores = 4, seed = 1234)
fit_cubic_brm <- brm(f_cubic, data = turnout, family = bernoulli,
backend = "cmdstanr", chains = 4, cores = 4, seed = 1234)Finally, compute the LOOIC.
# looic for linear model
loo_linear_brm <- loo(fit_linear_brm)
loo_linear_brm
Computed from 4000 by 2000 log-likelihood matrix.
Estimate SE
elpd_loo -1017.2 23.1
p_loo 5.3 0.2
looic 2034.4 46.1
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.3]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.# looic for cubic model
loo_cubic_brm <- loo(fit_cubic_brm)
loo_cubic_brm
Computed from 4000 by 2000 log-likelihood matrix.
Estimate SE
elpd_loo -1011.4 23.3
p_loo 7.3 0.3
looic 2022.8 46.5
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.5]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.# compare the two models
loo_compare(loo_linear_brm, loo_cubic_brm) model elpd_diff se_diff p_worse diag_diff diag_elpd
fit_cubic_brm 0.0 0.0 NA
fit_linear_brm -5.8 4.1 0.92 # AIC and BIC from glm() fits, for comparison
fit_linear_glm <- glm(f_linear, family = binomial, data = turnout)
fit_cubic_glm <- glm(f_cubic, family = binomial, data = turnout)
AIC(fit_linear_glm, fit_cubic_glm) df AIC
fit_linear_glm 5 2033.981
fit_cubic_glm 7 2022.419BIC(fit_linear_glm, fit_cubic_glm) df BIC
fit_linear_glm 5 2061.986
fit_cubic_glm 7 2061.626We haven’t described what the LOOIC is yet, but the cubic model is ahead by 5.8 elpd with an se_diff of 4.1. This is modest evidence for the cubic. (The AIC is 2022.4 for the cubic and 2034.0 for the linear, and the BIC is 2061.6 and 2062.0.)
Let’s dig into the LOOIC.
23.2 Recall the AIC
\[ \text{AIC} = -2 \ell(\hat{\theta}) + 2k \]
Here, we have two parts:
- \(-2\ell(\hat{\theta})\) measures in-sample fit.
- We can think of \(2k\) as the penalty for using the same data to fit and to evaluate the model.
23.3 The LOOIC idea
LOOIC is based on the logic of leave-one-out (LOO) cross-validation. We could simply leave out observation \(i\), fit to the other \(N - 1\) observations, and then ask: how much predictive density does the fit put on the held-out \(y_i\)? Then we could repeat for every \(i\).
Let’s see brute force leave-one-out cross-validation with ML. Refit the linear turnout model 2,000 times, dropping one respondent each time. The log-likelihoods of the held-out vote at the refit’s \(\hat{\theta}\) sum to −1017.3.
N <- nrow(turnout)
loo_ml <- numeric(N)
for (i in 1:N) {
fit <- glm(f_linear,
family = binomial,
data = turnout[-i, ])
p <- predict(fit,
newdata = turnout[i, ],
type = "response")
loo_ml[i] <- dbinom(
turnout$vote[i],
size = 1, prob = p,
log = TRUE)
}
sum(loo_ml)[1] -1017.293We formalize this intuition with the \(\text{elpd}_{\text{loo}}\).
\[ \text{elpd}_{\text{loo}} = \sum_{i=1}^{N} \log f(y_i \mid y_{-i}), \quad \text{where} \quad f(y_i \mid y_{-i}) = \int f(y_i \mid \theta)\, f(\theta \mid y_{-i})\, d\theta \]
- \(y_{-i}\) = every observation except \(i\)
- elpd = “expected log predictive density” (higher is better)
- But there’s a problem—we have no \(\hat{\theta}\) to plug in. So instead, we average \(y_i\)’s likelihood over the posterior fit without \(y_i\). This requires computing the posterior without \(y_i\), which means refitting the model \(N\) times.
Then we can define the LOOIC as follows.
\[ \text{LOOIC} = -2 \times \text{elpd}_{\text{loo}} \]
- The \(-2\) only puts LOO on AIC’s scale (i.e., lower is better!)
- No explicit penalty term needed, because of the cross-validation logic.
23.3.1 In-sample fit and \(p_{\text{loo}}\)
- Averaging over draws from the full-data posterior gives the log pointwise predictive density (lppd) \(= \sum_{i=1}^N \log \left( \frac{1}{S} \sum_{s=1}^S f(y_i \mid \theta^{(s)}) \right)\). This is just the posterior version of \(\ell(\hat{\theta})\). Each \(y_i\) helped produce the draws, so, like \(\ell(\hat{\theta})\), it is in-sample and too optimistic.
- We can compute \(p_{\text{loo}} = \text{lppd} - \text{elpd}_{\text{loo}}\), so \(\text{LOOIC} = -2\,\text{lppd} + 2\,p_{\text{loo}}\). Notice that this has the same form as \(\text{AIC} = -2\ell(\hat{\theta}) + 2k\). While AIC assumes the penalty, LOO measures it.
loo()calls \(p_{\text{loo}}\) the “effective number of parameters”.
# S × N matrix of log f(y_i | beta^(s)): 4,000 draws by 2,000 respondents
ll <- log_lik(fit_linear_brm)
dim(ll)[1] 4000 2000# lppd: average the likelihood over the draws, then take the log and sum
lppd <- sum(log(colMeans(exp(ll))))
lppd[1] -1011.864# p_loo: lppd minus elpd_loo
lppd - loo_linear_brm$estimates["elpd_loo", "Estimate"][1] 5.316335This matches the p_loo of 5.3 from loo() above.
But exact LOO is computationally expensive. Sampling from one posterior is usually time-consuming. Can we avoid sampling from N of them? Yes!
23.3.2 PSIS: LOO from one fit
Instead of dropping observations one-by-one, loo() uses Pareto-smoothed importance sampling (PSIS) (Vehtari, Gelman, and Gabry 2017) on the draws from the one full-data fit.
- Here’s the intuition of importance sampling: the full-data draws are “close” to draws from \(f(\theta \mid y_{-i})\), so we can weight each by \(1/f(y_i \mid \theta^{(s)})\). This weighting gives \(f(y_i \mid y_{-i}) \approx \left[ \frac{1}{S} \sum_{s=1}^S \frac{1}{f(y_i \mid \theta^{(s)})} \right]^{-1}\) (average the reciprocal of the likelihood, then take the reciprocal).1
1 This is exact, not a Jensen’s-inequality approximation. When the observations are independent given \(\theta\), \(f(\theta \mid y) \propto f(y_i \mid \theta)\, f(\theta \mid y_{-i})\), so the average of \(1/f(y_i \mid \theta)\) over the full-data posterior is exactly \(1/f(y_i \mid y_{-i})\). The \(\approx\) comes only from using \(S\) draws.
# average 1/f over the draws, take the reciprocal, take the log, and sum
elpd_is <- sum(-log(colMeans(exp(-ll))))
elpd_is[1] -1017.173This matches loo()’s −1017.2… with no refits!
But this raw version can be unstable. If a few draws give \(y_i\) a tiny likelihood, their weights \(1/f(y_i \mid \theta^{(s)})\) are huge, and those few draws dominate the average. PSIS replaces the largest weights with smoothed values from a Pareto distribution fit to them (Vehtari, Gelman, and Gabry 2017).
The estimated shape of that Pareto distribution, \(\hat{k}\), is a diagnostic, one per observation. If \(\hat{k} \le 0.7\), the PSIS estimate for that observation is reliable. If \(\hat{k} > 0.7\), it is not. For both turnout models, every \(\hat{k}\) is below 0.7.
A large \(\hat{k}\) usually means the observation is influential: the posterior changes a lot when we leave it out, so the full-data draws are a poor stand-in for draws from \(f(\theta \mid y_{-i})\) (Vehtari, Gelman, and Gabry 2017).
Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.
23.4 Example: Poisson vs. negative binomial, again
# load data for santiago
sant <- crdata::holland2015 |>
filter(city == "santiago") |>
# create rescaled columns (see below)
mutate(across(c(lower, vendors, budget, population), arm::rescale,
.names = "rs_{.col}"))
# formula corresponds to model 1 for each city in holland (2015) table 2
f <- operations ~ rs_lower + rs_vendors + rs_budget + rs_population
# poisson regression
pois_fit <- glm(f, family = poisson, data = sant)
# nb regression
nb_fit <- MASS::glm.nb(f, data = sant)
# AIC from wk05; df is k
AIC(pois_fit, nb_fit) df AIC
pois_fit 5 226.7905
nb_fit 6 133.6848# fit the same two models with brm(), using the default priors
pois_brm <- brm(f, data = sant, family = poisson,
iter = 3000, # 1,500 draws per chain after warmup
backend = "cmdstanr", chains = 4, cores = 4, seed = 1234)
nb_brm <- brm(f, data = sant, family = negbinomial,
iter = 3000,
backend = "cmdstanr", chains = 4, cores = 4, seed = 1234)We rescaled the predictors above. On their original scales (budget runs up to 1,118), {brms}’s random starting values put \(\log \mu\) in the hundreds or thousands, and some chains never recover (R-hat above 3). Rescaling changes the coefficients but not the fit, so the AIC is the same as before. Check convergence before computing LOO.
# check convergence (R-hat and ESS) before computing LOO
summary(pois_brm) Family: poisson
Links: mu = log
Formula: operations ~ rs_lower + rs_vendors + rs_budget + rs_population
Data: sant (Number of observations: 34)
Draws: 4 chains, each with iter = 3000; warmup = 1500; thin = 1;
total post-warmup draws = 6000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.62 0.14 0.32 0.88 1.00 3610 3422
rs_lower -1.42 0.34 -2.10 -0.76 1.00 2842 3014
rs_vendors -1.18 0.60 -2.45 -0.12 1.00 3051 3174
rs_budget -0.31 0.25 -0.79 0.18 1.00 3385 3093
rs_population 0.56 0.34 -0.10 1.25 1.00 3057 3199
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).summary(nb_brm) Family: negbinomial
Links: mu = log
Formula: operations ~ rs_lower + rs_vendors + rs_budget + rs_population
Data: sant (Number of observations: 34)
Draws: 4 chains, each with iter = 3000; warmup = 1500; thin = 1;
total post-warmup draws = 6000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.77 0.40 0.05 1.64 1.00 3401 3482
rs_lower -2.02 1.20 -4.53 0.15 1.00 3445 2925
rs_vendors -1.11 1.51 -4.09 1.98 1.00 3469 3447
rs_budget 0.09 1.28 -2.25 2.87 1.00 2758 2722
rs_population 1.42 1.70 -1.38 5.31 1.00 3014 2843
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 0.31 0.11 0.14 0.58 1.00 3025 4029
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Notice that R-hat is 1.00 for every parameter in both fits, and that the smallest bulk ESS is 2,758. The Poisson posterior means are also close to the glm() estimates. The NB posterior means differ more from the glm.nb() estimates, but by less than half a posterior SD (see the margin).
| term | ML | Posterior |
|---|---|---|
| Intercept | 0.66 | 0.62 |
| rs_lower | -1.41 | -1.42 |
| rs_vendors | -1.06 | -1.18 |
| rs_budget | -0.30 | -0.31 |
| rs_population | 0.54 | 0.56 |
| term | ML | Posterior |
|---|---|---|
| Intercept | 0.60 | 0.77 |
| rs_lower | -1.74 | -2.02 |
| rs_vendors | -0.90 | -1.11 |
| rs_budget | -0.14 | 0.09 |
| rs_population | 0.72 | 1.42 |
23.4.1 loo()
# LOO for each model
loo_pois <- loo(pois_brm)
loo_pois
Computed from 6000 by 34 log-likelihood matrix.
Estimate SE
elpd_loo -129.2 31.9
p_loo 29.2 10.5
looic 258.4 63.8
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 0.9]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 27 79.4% 438
(0.7, 1] (bad) 4 11.8% <NA>
(1, Inf) (very bad) 3 8.8% <NA>
See help('pareto-k-diagnostic') for details.loo_nb <- loo(nb_brm)
loo_nb
Computed from 6000 by 34 log-likelihood matrix.
Estimate SE
elpd_loo -68.5 12.0
p_loo 6.3 2.8
looic 136.9 24.1
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 0.8]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 32 94.1% 724
(0.7, 1] (bad) 2 5.9% <NA>
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.We understand the following:
elpd_loo(higher is better)p_loo(effective number of parameters)looic= −2 \(\times\)elpd_loo
However, for the Poisson, 7 of 34 observations have Pareto \(\hat{k}\) above 0.7. And the NB has 2 of 34. This means that our PSIS approximation might not be working well.
23.4.2 Fixing high \(\hat{k}\): refit for those observations
reloo_pois <- loo(pois_brm,
1 reloo = TRUE,
2 reloo_args = list(seed = 1234)
)
reloo_nb <- loo(nb_brm,
reloo = TRUE,
reloo_args = list(seed = 1234)
)- 1
- For each observation with \(\hat{k} > 0.7\), refit the model without that observation and compute exact LOO; use PSIS for the rest.
- 2
- Pass a seed to the refits, so the numbers reproduce.
reloo_pois
Computed from 6000 by 34 log-likelihood matrix.
Estimate SE
elpd_loo -131.0 32.5
p_loo 31.0 11.3
looic 262.1 65.0
------
MCSE of elpd_loo is 0.5.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 0.9]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.reloo_nb
Computed from 6000 by 34 log-likelihood matrix.
Estimate SE
elpd_loo -69.4 12.6
p_loo 7.2 3.7
looic 138.8 25.3
------
MCSE of elpd_loo is 0.2.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 0.8]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.Here, reloo = TRUE means that we compute LOO for the flagged observations and PSIS for the rest. It doesn’t take much time here (7 refits for the Poisson, and 2 for the NB).
As we can see, the LOOIC has changed a little. Before, the difference was 121.5, suggesting that the NB is much better. Now it’s 123.3, suggesting the same.
23.4.3 loo_compare()
# compare models with PSIS-LOO
loo_compare(loo_pois, loo_nb) model elpd_diff se_diff p_worse diag_diff diag_elpd
nb_brm 0.0 0.0 NA 2 k_psis > 0.7
pois_brm -60.8 22.6 1.00 N < 100 7 k_psis > 0.7# compare models after refitting the high-k observations
loo_compare(reloo_pois, reloo_nb) model elpd_diff se_diff p_worse diag_diff diag_elpd
nb_brm 0.0 0.0 NA
pois_brm -61.7 22.4 1.00 N < 100
ELPD or LOOIC?
loo_compare() reports differences on the ELPD scale. ELPD differences are half the size of LOOIC differences and have the opposite sign (higher ELPD is better; lower LOOIC is better). It’s fine to use either. But the LOOIC connects nicely to the AIC and BIC, so that’s the one I’ve emphasized here.
loo_compare()prints the best model in the top row, andelpd_diffis relative to that best model.elpd_diffis on the elpd scale, so multiply by (−2) for the LOOIC scale (e.g., −61.65 elpd becomes 123.3 LOOIC).se_diffis the SE of the difference.
Observations:
- After using
reloo = TRUE, NB beats Poisson by 61.7 elpd (se_diff of 22.4, so about 2.8 SEs). With PSIS alone, the difference was 60.8 (se_diff of 22.6, so about 2.7 SEs). This is a huge difference on the elpd scale, but the SE is also really wide. p_worseis a normal-approximation probability that a model predicts worse than the top model (i.e.,pnorm(0, elpd_diff, se_diff)). It’s 1.00 for the Poisson, which is a bad sign for that model.- The
diag_diffcolumn warns “N < 100.” With only 34 observations,se_diffis probably too small, sop_worseis probably too close to 1. The loo documentation suggests doublingse_diffas a rough, conservative fix (see?loo-glossary). Doubled, the difference is about 61.7 / 44.8 ≈ 1.4 SEs.
Raftery’s table in Information Criteria calls a BIC difference larger than 10 “very strong” evidence. The BIC, like the AIC and LOOIC, is on the \(-2 \times \log\) scale, so a BIC difference of 10 is the same size as a difference of 5 on the elpd scale.
23.4.4 Interpreting elpd_diff and se_diff
Some rules of thumb:
- For |
elpd_diff| < 4, the models predict about equally well, regardless of se_diff. - For |
elpd_diff| > 4 AND large compared tose_diff: a meaningful difference (see?loo-glossary).2
2 Compare Raftery’s table for BIC differences in Information Criteria, where we established “about 10” as a big difference between the models. Remember that LOOIC differences are twice elpd differences, so 4 elpd is the same as 8 LOOIC.
23.4.5 LOOIC next to AIC
# AIC and k from the ML fits; LOOIC and p_loo from the brm() fits (after reloo)
tibble(
model = c("Poisson", "Negative binomial"),
k = AIC(pois_fit, nb_fit)$df,
p_loo = c(reloo_pois$estimates["p_loo", "Estimate"],
reloo_nb$estimates["p_loo", "Estimate"]),
AIC = AIC(pois_fit, nb_fit)$AIC,
LOOIC = c(reloo_pois$estimates["looic", "Estimate"],
reloo_nb$estimates["looic", "Estimate"])
) |>
mutate(across(p_loo:LOOIC, \(x) round(x, 1))) |>
as.data.frame() model k p_loo AIC LOOIC
1 Poisson 5 31.0 226.8 262.1
2 Negative binomial 6 7.2 133.7 138.8NB is much better by both criteria. The AIC difference is 93.1. The LOOIC difference is 123.3.
Notice that the NB’s \(p_{\text{loo}}\) (7.2) is close to its \(k\) (6), but the Poisson’s \(p_{\text{loo}}\) (31.0) is about six times its \(k\) (5). The AIC assumes the penalty, and LOO measures it.