A single fitted model tells you what it thinks; a set of fitted models makes you choose. The trouble is that the obvious score, the fit to the data you already have, always rewards the more complicated model. Add a parameter and the likelihood can only go up, so the in-sample deviance points you at the most flexible candidate every time, whether or not that flexibility describes anything real.
The frequentist answer is to add a penalty to the maximised fit; that is what AIC does, using a single point estimate of the parameters. The Bayesian versions use the whole posterior instead of a point: the deviance information criterion (DIC), the widely applicable information criterion (WAIC), and leave-one-out cross-validation approximated by Pareto-smoothed importance sampling (PSIS-LOO). This post codes all three by hand from posterior draws, on a small ecological count problem, and shows why the in-sample number is not a model-selection tool at all.
Abundance often peaks at an intermediate value of some gradient and falls away at both ends, which a quadratic term on the log scale captures. We simulate counts from exactly that structure and then pretend we do not know it.
set.seed(0)
n <- 80
b_true <- c(1.0, 0.8, -0.7) # intercept, linear, quadratic
xsim <- sort(runif(n, -2.2, 2.2))
Xfull <- cbind(1, xsim, xsim^2)
y <- rpois(n, exp(as.vector(Xfull %*% b_true)))
c(n = n, mean_count = round(mean(y), 2), zeros = sum(y == 0), max = max(y)) n mean_count zeros max
80.00 1.74 23.00 8.00 The candidates are nested: intercept only (M0), linear (M1), the true quadratic (M2), and a cubic model (M3) that has one parameter too many.
Every criterion below needs posterior draws, so we sample each model with a random-walk Metropolis sampler. The log-posterior is a Poisson log-likelihood plus a weakly-informative N(0, 5^2) prior on each coefficient. The only shortcut is the proposal: we shape it with the inverse Hessian at the maximum, scaled by the standard 2.4^2 / d factor, so the chain mixes without hand-tuning. The Metropolis step itself is coded out in full.
loglik_vec <- function(beta, X) dpois(y, exp(as.vector(X %*% beta)), log = TRUE)
logpost <- function(beta, X, ps = 5) {
lam <- exp(as.vector(X %*% beta))
if (any(!is.finite(lam)) || any(lam > 1e8)) return(-Inf)
sum(dpois(y, lam, log = TRUE)) + sum(dnorm(beta, 0, ps, log = TRUE))
}
run_mcmc <- function(X, iter = 30000, warm = 5000, seed = 1) {
set.seed(seed)
d <- ncol(X)
fit <- optim(rep(0, d), function(b) -logpost(b, X), method = "BFGS", hessian = TRUE)
Cp <- chol(solve(fit$hessian) * (2.4^2 / d)) # proposal Cholesky
beta <- fit$par; lp <- logpost(beta, X)
out <- matrix(NA_real_, iter, d); acc <- 0
for (t in seq_len(iter)) {
prop <- beta + as.vector(crossprod(Cp, rnorm(d)))
lpp <- logpost(prop, X)
if (log(runif(1)) < lpp - lp) { beta <- prop; lp <- lpp; acc <- acc + 1 }
out[t, ] <- beta
}
draws <- out[(warm + 1):iter, , drop = FALSE]
ll <- matrix(NA_real_, nrow(draws), length(y)) # S x n log-likelihood
for (s in seq_len(nrow(draws))) ll[s, ] <- loglik_vec(draws[s, ], X)
list(draws = draws, ll = ll, acc = acc / iter)
}
Xlist <- list(M0 = Xfull[, 1, drop = FALSE], M1 = Xfull[, 1:2],
M2 = Xfull[, 1:3], M3 = cbind(Xfull, xsim^3))
fits <- Map(function(X, s) run_mcmc(X, seed = s), Xlist, c(101, 202, 303, 404))The pointwise log-likelihood matrix ll, with one row per draw and one column per observation, is the single object all three criteria are built from. The acceptance rates sit near the random-walk optimum (0.44, 0.35, 0.32, 0.30 for M0 to M3), and the quadratic model recovers the truth: its posterior mean is 1.065, 0.824, -0.806 against the generating 1, 0.8, -0.7.
The deviance at the maximum, -2 times the maximised log-likelihood, is the fit component that AIC starts from. Because the models are nested, it can only fall as parameters are added.
naive_dev <- sapply(Xlist, function(X) {
Xm <- X; -2 * as.numeric(logLik(glm(y ~ Xm - 1, family = poisson)))
})
round(naive_dev, 1) M0 M1 M2 M3
295.4 281.4 222.6 222.6 The linear-to-quadratic drop is real: leaving out the curvature costs a great deal. But the quadratic-to-cubic step does not move at all (222.6 for both M2 and M3): the extra parameter buys nothing here, yet the in-sample rule is at best indifferent between them and would pick whichever is lower by chance. It never charges you for the flexibility. That is the whole problem, and every criterion below is a way of estimating the price.
DIC uses the posterior distribution of the deviance. Let Dbar be the posterior mean of -2 log-likelihood and let D(thetabar) be the deviance at the posterior mean of the parameters. The effective number of parameters is their difference, pD = Dbar - D(thetabar), and DIC = Dbar + pD (Spiegelhalter 2002).
dic_one <- function(m, X) {
Dbar <- mean(-2 * rowSums(m$ll))
Dhat <- -2 * sum(loglik_vec(colMeans(m$draws), X))
pD <- Dbar - Dhat
c(DIC = Dbar + pD, pD = pD)
}WAIC works observation by observation and never touches a point estimate. The log pointwise predictive density averages the likelihood over the posterior at each point, lppd_i = log mean_s p(y_i | theta_s); the effective parameter count is the posterior variance of the log-likelihood at each point, pWAIC_i = var_s log p(y_i | theta_s). Then WAIC = -2 (sum_i lppd_i - sum_i pWAIC_i) (Watanabe’s criterion in the Bayesian framing of Gelman 2014).
waic_one <- function(m) {
cm <- apply(m$ll, 2, max) # log-sum-exp, stably
lppd <- cm + log(colMeans(exp(sweep(m$ll, 2, cm))))
pw <- apply(m$ll, 2, var)
list(WAIC = -2 * sum(lppd - pw), pWAIC = sum(pw),
pt = lppd - pw) # pointwise elpd
}Leave-one-out cross-validation asks how well the model predicts each point when that point is held out. Refitting n times is expensive, but the held-out predictive density can be estimated from the single full-data fit by importance sampling: the weight for draw s at point i is 1 / p(y_i | theta_s). Those raw weights have heavy tails and occasionally one draw dominates, which is where PSIS-LOO improves on plain importance sampling (Vehtari 2017): it fits a generalised Pareto distribution to the largest weights and replaces them with the fitted quantiles. The estimated Pareto shape, k-hat, is also a diagnostic: k > 0.7 means the importance-sampling estimate for that point is not to be trusted.
The Pareto fit uses the empirical-Bayes estimator of Zhang and Stephens (2009), coded here in full.
gpdfit <- function(x) { # x: exceedances, ascending, > 0
N <- length(x); m <- 30 + floor(sqrt(N))
bs <- 1 - sqrt(m / (seq_len(m) - 0.5))
bs <- bs / (3 * x[floor(N / 4 + 0.5)]) + 1 / x[N]
ks <- vapply(bs, function(b) mean(log1p(-b * x)), numeric(1))
L <- N * (log(-bs / ks) - ks - 1)
w <- 1 / vapply(seq_len(m), function(j) sum(exp(L - L[j])), numeric(1))
b <- sum(bs * w); k <- mean(log1p(-b * x))
list(k = (k * N + 5) / (N + 10), sigma = -k / b) # prior-shrunk k, scale
}
qgpd <- function(p, k, sig)
if (abs(k) < 1e-9) -sig * log1p(-p) else sig / k * ((1 - p)^(-k) - 1)
psis_i <- function(logr) { # log importance ratios, one point
S <- length(logr); r <- exp(logr - max(logr))
M <- ceiling(min(0.2 * S, 3 * sqrt(S)))
ord <- order(r); u <- r[ord][S - M]; tail_i <- ord[(S - M + 1):S]
g <- gpdfit(r[tail_i] - u)
r[tail_i] <- pmin(u + qgpd((seq_len(M) - 0.5) / M, g$k, g$sigma), max(r))
list(w = r / sum(r), k = g$k)
}
loo_one <- function(m) {
nn <- ncol(m$ll); ei <- numeric(nn); kh <- numeric(nn)
for (i in seq_len(nn)) {
ps <- psis_i(-m$ll[, i]) # ratios 1 / p_i
a <- m$ll[, i]; mx <- max(a)
ei[i] <- mx + log(sum(ps$w * exp(a - mx))) # weighted predictive
kh[i] <- ps$k
}
list(LOO = -2 * sum(ei), elpd = ei, khat = kh)
}ics <- lapply(names(Xlist), function(nm) {
m <- fits[[nm]]; X <- Xlist[[nm]]
d <- dic_one(m, X); w <- waic_one(m); l <- loo_one(m)
list(in_sample = naive_dev[[nm]], DIC = d["DIC"], pD = d["pD"],
WAIC = w$WAIC, pWAIC = w$pWAIC, pt = w$pt,
LOO = l$LOO, khat = l$khat, p_loo = w$pWAIC)
})
names(ics) <- names(Xlist)
aic <- sapply(Xlist, function(X) { Xm <- X; AIC(glm(y ~ Xm - 1, family = poisson)) })
tab <- data.frame(
model = names(Xlist), parameters = sapply(Xlist, ncol),
in_sample = round(naive_dev, 1),
DIC = round(sapply(ics, `[[`, "DIC"), 1), pD = round(sapply(ics, `[[`, "pD"), 2),
WAIC = round(sapply(ics, `[[`, "WAIC"), 1),
pWAIC = round(sapply(ics, `[[`, "pWAIC"), 2),
LOO = round(sapply(ics, `[[`, "LOO"), 1), AIC = round(aic, 1),
max_k = round(sapply(ics, function(z) max(z$khat)), 2), row.names = NULL)
tab model parameters in_sample DIC pD WAIC pWAIC LOO AIC max_k
1 M0 1 295.4 297.4 0.99 298.1 1.71 298.1 297.4 0.10
2 M1 2 281.4 285.4 1.99 286.4 2.89 286.4 285.4 0.13
3 M2 3 222.6 228.5 2.95 228.2 2.53 228.3 228.6 0.23
4 M3 4 222.6 230.6 3.94 230.1 3.28 230.1 230.6 0.38Every criterion that estimates out-of-sample loss agrees, and disagrees with the in-sample number. In-sample deviance is flat across the last step; DIC, WAIC, PSIS-LOO and AIC all put the minimum at the true quadratic model and charge the cubic model roughly two deviance units for its spare parameter. The effective parameter counts (pWAIC and pD) track the actual dimensions, and no Pareto k-hat exceeds 0.38, so the importance sampling is stable throughout.
np <- sapply(Xlist, ncol)
seW <- sapply(names(Xlist), function(nm) 2 * sqrt(n * var(ics[[nm]]$pt)))
mk <- function(v, lab) data.frame(npar = np, score = v, kind = lab)
df1 <- rbind(mk(naive_dev, "in-sample (-2 log-lik)"), mk(sapply(ics, `[[`, "DIC"), "DIC"),
mk(sapply(ics, `[[`, "WAIC"), "WAIC"), mk(sapply(ics, `[[`, "LOO"), "PSIS-LOO"))
df1$kind <- factor(df1$kind, levels = c("in-sample (-2 log-lik)", "DIC", "WAIC", "PSIS-LOO"))
best <- which.min(sapply(ics, `[[`, "WAIC"))
dfW <- data.frame(npar = np, WAIC = sapply(ics, `[[`, "WAIC"), se = seW)
ggplot(df1, aes(npar, score, colour = kind)) +
geom_errorbar(data = dfW, inherit.aes = FALSE,
aes(npar, ymin = WAIC - se, ymax = WAIC + se), width = 0.08,
colour = te$forest, alpha = 0.6) +
geom_line(linewidth = 0.9) + geom_point(size = 2.4) +
annotate("point", x = np[best], y = sapply(ics, `[[`, "WAIC")[best],
size = 5, shape = 21, stroke = 1.1, colour = te$ink, fill = NA) +
scale_colour_manual(values = c("in-sample (-2 log-lik)" = te$terra, "DIC" = te$gold,
"WAIC" = te$forest, "PSIS-LOO" = te$sage), name = NULL) +
scale_x_continuous(breaks = np, labels = paste0(names(Xlist), " (", np, "p)")) +
labs(x = "candidate model", y = "deviance scale (lower is better)",
title = "In-sample fit keeps improving; WAIC, LOO and DIC do not",
subtitle = "Poisson counts, true structure is the quadratic model M2") +
theme_te() + theme(legend.position = "top")
The size of a difference matters as much as its sign, and PSIS-LOO carries a standard error for it. The gap between the quadratic and the cubic model is small but real (an elpd difference of 0.93 with a standard error of 0.24), so we prefer the simpler true model without pretending the cubic is badly wrong. Dropping to the linear model is a different matter: that difference is 29.1 with a standard error of 6.9, decisive by any reading.
The k-hat values are worth plotting in their own right, because they tell you when the LOO number can be believed. Here they climb gently with the size of the count, which is expected: a large observation pulls harder on the fit and is therefore closer to being influential. None of them crosses even the cautionary 0.5 line.
kk <- ics[[best]]$khat
ggplot(data.frame(y = y, k = kk), aes(y, k)) +
geom_hline(yintercept = 0.7, linetype = "dashed", colour = te$terra) +
geom_hline(yintercept = 0.5, linetype = "dotted", colour = te$gold) +
geom_jitter(width = 0.12, height = 0, size = 2.3, colour = te$forest, alpha = 0.85) +
annotate("text", x = max(y), y = 0.72, hjust = 1, vjust = 0, size = 3.4,
colour = te$terra, label = "k = 0.7: importance sampling unreliable") +
annotate("text", x = max(y), y = 0.52, hjust = 1, vjust = 0, size = 3.4,
colour = te$gold, label = "k = 0.5") +
coord_cartesian(ylim = c(min(kk) - 0.05, 0.85)) +
labs(x = "observed count y", y = "Pareto k-hat",
title = "PSIS diagnostic for the selected model M2",
subtitle = "every k below 0.5: the LOO approximation is trustworthy here") +
theme_te()
For most models fitted by MCMC, PSIS-LOO is the one to report: it estimates the same predictive quantity as WAIC but degrades more gracefully with weak priors or influential points, and its k-hat values tell you exactly where it is failing. WAIC is a cheaper approximation to the same target and agrees with LOO here, as it should in a well-behaved problem. DIC is the oldest of the three and the least safe, because it collapses the posterior to a single point through D(thetabar); in mixture or weakly-identified models that plug-in can misbehave, and pD can even go negative. All three rest on splitting the data into the units you want to predict, which for these independent counts is one observation each; for grouped or time-structured data the unit is a decision in itself, not a default.
The deeper point is the one the in-sample column makes for free. Fit to the data in hand is not evidence about a model, because it cannot see the cost of flexibility. Only an estimate of out-of-sample loss can, and the Bayesian versions build that estimate out of the same posterior draws you already have.