Model Comparison: Ranking, Not Testing

Fit is not enough — every parameter must pay rent.

When tests run out

The likelihood-ratio test compares exactly two models, one nested inside the other, with a null-hypothesis ritual attached. But the everyday question is usually broader: of several candidate models, which should I use?

Consider an expanded version of the heart-rate study: $n = 150$ participants, and this time three things were measured alongside resting heart rate — weekly exercise hours, age, and (because the research assistant got carried away) shoe size. Four nested candidates present themselves, from intercept-only to everything-in.

Log-likelihood alone cannot referee this contest, because $\ell$ never goes down when a predictor is added — the bigger model can always at worst set the new coefficient to zero. Even shoe size buys a sliver of fit ($+0.0033$ log-likelihood units). Judged on fit alone, the biggest model always wins, junk and all.

The information criteria charge admission:

$$\text{AIC} = -2\ell(\hat\beta) + 2k \qquad\qquad \text{BIC} = -2\ell(\hat\beta) + k\,\ln(n)$$

Rent per parameter. On the $-2\ell$ scale, a pure-noise predictor improves fit by about $1$ on average (half a $\chi^2_1$, doubled) — but AIC charges $2$ for it, and BIC charges $\ln(n)$ ($\approx 5$ at $n = 150$). Useless parameters lose money. Here $k$ counts every estimated parameter including $\sigma$, matching R's AIC(); only differences between models matter, so any consistent convention gives the same ranking.

Interactive: the ladder

Each bar is one model's score (lower is better). The dark segment is the misfit, $-2\ell$; the pale segment on top is the penalty. Adding exercise slashes the misfit — worth every penny of its rent. Adding age helps more modestly. Adding shoe size barely dents the misfit while paying full rent, so its bar goes back up: the signature of a junk predictor. Toggle the criterion and shrink the sample to watch the verdicts shift.

150
misfit (−2ℓ) penalty (rent) winner (lowest total)
Model k logLik AIC ΔAIC BIC ΔBIC
loading…

Rules of thumb (Burnham & Anderson): models within $\Delta \le 2$ of the best are essentially tied; $\Delta \ge 10$ means essentially no support.

Two philosophies, occasionally two answers

AIC and BIC look like siblings but ask different questions. AIC asks: which model will predict new data best? (It approximates out-of-sample prediction error — Kullback–Leibler divergence.) BIC asks: which model most plausibly generated these data? (It approximates a Bayesian posterior model probability.) BIC's rent, $\ln n$, exceeds AIC's flat 2 once $n \ge 8$, so BIC pushes harder towards small models — increasingly so as data accumulate.

The slider above shows a live disagreement: at $n = 50$, AIC scores the exercise-only and exercise-plus-age models as essentially tied ($314.70$ vs $314.77$) — too little data to resolve whether age matters for prediction — while BIC already votes clearly for the simpler model ($320.44$ vs $322.42$). By $n = 150$ both agree that age has earned its place and shoe size hasn't.

Not a test. Nothing here is rejected, no p-value is produced, and the candidates need not be nested. AIC/BIC simply rank — which is exactly what you want when the question is “which model do I take forward?” rather than “is this coefficient zero?”. You met this move in Tutorial 4, where AIC adjudicated Poisson against Negative Binomial for the overdispersed counts — non-nested rivals no LR test could referee directly.

The sampling branch ranks models too, with its own tools — Bayes factors and leave-one-out cross-validation (LOO, WAIC) play the roles that AIC and BIC play here; the Bayesian page shows the machinery those tools are built on.

In code: one line, four verdicts

Outputs below are from real runs on the expanded dataset (docs/data/inference-model-comparison.json, seed 4242) — validation scripts scripts/R/validate-inference.R, scripts/py/generate_inference_data.py.

m0 <- lm(y ~ 1,                      data = hr2)
m1 <- lm(y ~ exercise,               data = hr2)
m2 <- lm(y ~ exercise + age,         data = hr2)
m3 <- lm(y ~ exercise + age + shoe,  data = hr2)

AIC(m0, m1, m2, m3)
BIC(m0, m1, m2, m3)

#    df       AIC            df       BIC
# m0  2 1011.2890        m0  2 1017.3103
# m1  3  933.9958        m1  3  943.0277
# m2  4  921.3732   <--  m2  4  933.4157   <-- winner both ways at n = 150
# m3  5  923.3666        m3  5  938.4198
import numpy as np
import statsmodels.api as sm

for name, cols in [("M0", []), ("M1", ["exercise"]),
                   ("M2", ["exercise", "age"]),
                   ("M3", ["exercise", "age", "shoe"])]:
    X = sm.add_constant(hr2[cols]) if cols else np.ones((len(hr2), 1))
    fit = sm.OLS(hr2["y"], X).fit()
    k = len(cols) + 2          # coefficients + sigma, matching R's AIC()
    print(name, -2 * fit.llf + 2 * k, -2 * fit.llf + np.log(len(hr2)) * k)

# M0 1011.2890 1017.3103
# M1  933.9958  943.0277
# M2  921.3732  933.4157   <-- winner
# M3  923.3666  938.4198
# (statsmodels' own .aic omits sigma from k — a constant offset of 2
#  that never changes a ranking)