Model comparison with WAIC

The Widely Applicable Information Criterion (WAIC) (Watanabe, 2010; Gelman, Hwang & Vehtari, 2014) is a Bayesian model-comparison criterion that estimates out-of-sample predictive accuracy from the posterior log-likelihood matrix. Like PSIS-LOO, it returns an estimate of the expected log predictive density (ELPD), but it uses the variance of the log-likelihood across posterior draws instead of importance-sampling corrections. This makes WAIC cheaper to compute and free of the Pareto k diagnostic concerns that can arise with LOO.

When to use WAIC

Use model.waic() when:

  • you want a quick Bayesian model-comparison metric without the importance-sampling overhead of LOO
  • you are comparing several Bayesian fits (method="VI" or method="MCMC") on the same data
  • none of your observations are highly influential (if some are, LOO’s per-observation diagnostics are more informative)

WAIC and LOO are asymptotically equivalent and typically give very similar results. LOO is generally preferred when Pareto k diagnostics are clean, because it is more robust to influential observations. WAIC is a good first-pass alternative, especially when fitting many candidate models.

For frequentist fits, use model.aic and model.bic instead.

Basic usage

Fit two competing models and call .waic() on each.

import numpy as np
import whittaker as wk
from whittaker.families.poisson import Poisson

rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2 * np.pi, n)
lam = np.exp(1.5 * np.sin(x))
y = rng.poisson(lam).astype(float)

data = {"x": x, "y": y}

# Model 1: smooth effect of x
m1 = wk.GAM("y ~ s(x)", family=Poisson()).fit(data, method="VI")

# Model 2: linear effect of x only
m2 = wk.GAM("y ~ x", family=Poisson()).fit(data, method="VI")

w1 = m1.waic(seed=0)
w2 = m2.waic(seed=0)
print(w1)
WAICResult
  ELPD_WAIC: -405.40  (SE 17.58)
  p_WAIC:    7.77
  WAIC:      810.80

Understanding the WAIC result

The WAICResult object contains:

elpd_waic is the summed expected log pointwise predictive density. Higher values mean better predictive accuracy.

se_elpd_waic is the standard error, computed as \sqrt{n \cdot \text{Var}(\text{pointwise})}.

p_waic is the effective number of parameters, computed as the sum of the per-observation variance of the log-likelihood across posterior draws. It measures how much the posterior predictions vary from observation to observation. When p_waic is much larger than the actual parameter count, it suggests the model may be overfit or misspecified.

waic is the WAIC on the deviance scale: -2 \cdot \text{ELPD}_\text{WAIC}. Lower is better. This scale matches the familiar AIC/BIC convention.

print(f"ELPD_WAIC:     {w1.elpd_waic:.2f} (SE {w1.se_elpd_waic:.2f})")
print(f"p_WAIC:        {w1.p_waic:.2f}")
print(f"WAIC:          {w1.waic:.2f}")
print(f"Observations:  {len(w1.pointwise)}")
ELPD_WAIC:     -405.40 (SE 17.58)
p_WAIC:        7.77
WAIC:          810.80
Observations:  300

Comparing models

waic_compare() computes the paired difference in ELPD between two models, using the pointwise values to account for correlation across observations.

from whittaker import waic_compare

cmp = waic_compare(w1, w2)
print(cmp)
WAICComparison
  ELPD diff: +86.11  (SE 15.57)
  model 1 preferred

A positive elpd_diff means the first model is preferred. As with LOO, a difference larger than two standard errors is strong evidence.

ratio = abs(cmp.elpd_diff) / cmp.se_diff
print(f"|ELPD diff| / SE = {ratio:.1f}")
if ratio > 2:
    preferred = "model 1 (smooth)" if cmp.elpd_diff > 0 else "model 2 (linear)"
    print(f"Strong evidence for {preferred}")
else:
    print("Models are not clearly distinguishable")
|ELPD diff| / SE = 5.5
Strong evidence for model 1 (smooth)

WAIC vs LOO

WAIC and LOO estimate the same quantity (ELPD) and are asymptotically equivalent. For well-behaved models they typically agree closely:

import warnings

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    loo1 = m1.loo(seed=0)

print(f"ELPD_WAIC: {w1.elpd_waic:.2f}  (SE {w1.se_elpd_waic:.2f})")
print(f"ELPD_LOO:  {loo1.elpd_loo:.2f}  (SE {loo1.se_elpd_loo:.2f})")
print(f"p_WAIC:    {w1.p_waic:.2f}")
print(f"p_LOO:     {loo1.p_loo:.2f}")
ELPD_WAIC: -405.40  (SE 17.58)
ELPD_LOO:  -405.48  (SE 17.59)
p_WAIC:    7.77
p_LOO:     7.44

The key differences:

WAIC PSIS-LOO
Speed Fast (just means and variances) Requires PSIS smoothing per observation
Diagnostics No per-observation diagnostics Pareto k flags unreliable observations
Robustness Sensitive to influential observations PSIS stabilizes extreme weights
Recommendation Good first pass Preferred when Pareto k diagnostics are clean

WAIC with MCMC fits

The workflow is the same for MCMC fits. All stored posterior draws are used automatically.

model = wk.GAM("y ~ s(x)", family=Poisson()).fit(data, method="MCMC")
w = model.waic()

Where to go next