Bayesian Model Comparison with LOO and WAIC
Lesson 6 ended with a confession: the posterior predictive check judged each model against the same sixty days that trained it, and a judge who helped write the answers is lenient. It also left a gap. A check can expose a model (the Normal, caught believing in negative kit sales), but it cannot rank rivals. Asha still has to pick the model that prices next month's moss order.
This lesson removes the conflict of interest. Hide one day, make each model state a probability for what that day actually sold, score the log of that probability, and repeat until all sixty days have taken a turn as the judge. That is leave-one-out cross-validation, LOO for short. Then you will meet WAIC, the shortcut that reads nearly the same answer straight off one matrix, with no refitting at all.
By the end of this lesson you will be able to:
- Score a model's forecast of a held-out day with the log predictive density, and say why the log scale is the honest ruler
- Compute exact leave-one-out scores for two conjugate models in base R, one line of update per left-out day
- Build the log-likelihood matrix, read the flattery in lppd, and charge for it with the WAIC penalty
- Compare two models with an elpd difference and its standard error, and state precisely what the winner did and did not win
Prerequisites: Lessons 1 to 6 of this course (the conjugate updates of Lesson 2, posterior draws, and Lesson 6's kit data, both fitted models, and the double-use-of-data warning), plus base R functions, sapply() and apply().
Below is the engine of the whole lesson. The strip is a dataset cut into folds; step through and watch each fold take one turn as the held-out judge while the rest do the training. Leave-one-out simply turns that dial to its limit: as many folds as days, so every one of Asha's sixty days gets scored by a model that never met it.
Two rivals, one moss order
Asha's question has not moved: sixty logged days of terrarium-kit sales, $39 a kit, live moss that needs re-misting on every day nothing sells, and a moss order to price for next month. What moved is the cast. Lesson 6 caught the rounded Normal red-handed: only 3 of its invented runs in 100 could match her fifteen blank days. The Gamma-Poisson survived. But a colleague looks at the same table and pushes back: the Normal passed on the mean, the spread, and the biggest day, three features out of four. Is one bad feature disqualifying? The debate needs a different kind of answer: not another feature, but one total score per model, earned on days the model never saw.
First, rebuild both rivals exactly as Lesson 6 left them. Each lesson runs in a fresh interactive R session, so run this block once:
Both posteriors agree about the center: demand near 1.68 kits a day. Whatever separates these two models, it is not the average. It is the shape of the days each one expects around that average, and to grade a shape we need a scoring rule.