Checking and comparing models
A fit that converged tells us the sampler found the posterior of the model we wrote. It does not tell us the model is a good description of what happened. This page takes the converged SIR fit from Fitting the boarding-school outbreak and asks that second question three ways: by comparing its predictions with the data it was fitted to, by fitting a competing model and scoring the two, and by checking an inferred quantity against a number the fit never saw. The first two point to the same weakness and suggest a repair; the third finds a problem the repair does not touch.
The commands run in a copy of docs/guide/first-fit/, as on the first page. Each page builds its own fits, so this one starts by repeating the first page’s fit and its predictive check:
camdl fit run fit-r0.toml --label r0 --seed 4
camdl fit predict @r0...
✓ converged
max R̂ 1.002 (R0), threshold 1.05
...
R̂: below 1.05 counts as converged; above 1.1 draws a finding, a warning up to 1.5
...
fit predict: free_forward horizon — subsampling 200 of 800 posterior draws (raise with --n-draws)
fit predict: one_step horizon — subsampling 200 of 800 posterior draws (raise with --n-draws)
...Where the Poisson model misses
uv run plot_predictive.py ppc-r0.svg Poisson=results/fits/fit_r0-*/predict/*/The two panels answer different questions. The free-forward band (left) simulates each epidemic from day 0 using only posterior draws of the parameters, so it asks whether the model can produce an epidemic like this one. It is wide, because an epidemic that takes off a day early or late stays early or late, and 13 of the 14 counts fall inside its 90% band. The one-day-ahead band (right) predicts each day’s count from the particle filter’s estimate of the epidemic given the counts before it, so it asks whether the model knows how much a single day’s count can vary. It is narrow, and only 10 of the 14 counts fall inside it. The misses are all in the decline, where the observed counts first stay above the model’s median and then fall below it: the school’s epidemic ended more abruptly than the model’s.
A free-forward check can pass while a one-day-ahead check fails, as here, because uncertainty about when the epidemic took off dominates the free-forward band and hides how much any one day may scatter. Both checks test particular consequences of the model; passing one does not show that the model is right.
Two things could produce the one-day-ahead misses. The observation model may be too tight, or the SIR dynamics may be wrong in the decline. The first is easier to test.
A model that allows noisier counts
A Poisson count has variance equal to its mean. With 300 boys in bed, the model expects the recorded count to lie within about \pm 2\sqrt{300}\approx\pm35 of the number infected. Counts read off a graph of a school infirmary are unlikely to be that precise. The negative binomial distribution adds a dispersion parameter k, and the observation model becomes
y_t \mid x, \theta \sim \operatorname{NegBin}\bigl(\text{mean} = I_t(x),\ \text{dispersion} = k\bigr), \qquad \operatorname{Var}(y_t \mid x, \theta) = I_t(x) + \frac{I_t(x)^2}{k}.
A small k allows much more scatter; as k\to\infty the distribution becomes Poisson. The process model, and therefore the distribution of epidemics, is the same as before; only the link from an epidemic to the counts changes.
first-fit/sir-negbin.camdl (excerpt)
parameters {
#' basic reproduction number: boys one case infects in a fully susceptible school
R0 : positive ~ log_normal(mu = 1.0, sigma = 0.5)
#' recovery rate; a boy stays infected 1 / gamma days on average
gamma : rate ~ log_normal(mu = -0.69, sigma = 0.42)
#' dispersion of the daily count: variance = mean + mean^2 / k
k : positive ~ log_normal(mu = 3.0, sigma = 1.0)
#' boys at risk
N0 : count
#' boys ill on day 0
I0 : count
}first-fit/sir-negbin.camdl (excerpt)
observations {
in_bed {
columns { time : time, in_bed : count }
projected = I
emit_schedule = every 1 'days
in_bed ~ neg_binomial(mean = projected, r = k)
}
}The neg_binomial distribution’s r argument is k. Its prior is log-normal with median 20 and 90% of its weight between about 4 and 100. At a mean of 300 boys, that range lets the standard deviation of a count be between about twice and about nine times the Poisson one, since the ratio is \sqrt{1 + 300/k}; at small counts every k in the range is close to Poisson. fit-negbin.toml is fit-r0.toml with the new model and a start value for k.
camdl fit run fit-negbin.toml --label negbin --seed 4...
✓ in_bed I NegBinomial
note: NegBinomial on a prevalence / snapshot projection is valid
but uncommon. Binomial or Poisson is the typical choice for
point-in-time counts. See camdl-run-spec.md §14.4.
...
✓ converged
max R̂ 1.005 (gamma), threshold 1.05
...
R̂: below 1.05 counts as converged; above 1.1 draws a finding, a warning up to 1.5
...
0 32.9% ✓ 7.493e-1 0 / 1000
1 35.8% ✓ 7.493e-1 0 / 1000
2 31.8% ✓ 8.586e-1 0 / 1000
3 31.1% ✓ 6.926e-1 0 / 1000
...
R0 Rhat=1.003 ✓ ESS bulk=575 tail=628 (71.8% of 800 draws)
gamma Rhat=1.005 ✓ ESS bulk=522 tail=654 (65.2% of 800 draws)
k Rhat=1.002 ✓ ESS bulk=562 tail=634 (70.3% of 800 draws)
...The chains converge from the same seed. camdl notes that a negative binomial on a count of people currently infected, rather than on new cases, is uncommon. The usual reason to use one is that transmission itself is overdispersed, which belongs in the process model; here the extra scatter stands for error in how the counts were recorded and read off the graph, which is an observation process, so the choice is deliberate.
camdl fit summary @negbin...
max R̂ 1.005 (gamma), threshold 1.05
0 of 3 parameters above threshold · chains' log-likelihoods agree
0/800 forkable
── posterior ───────────────────────────────────────────────────────────────
param mean ESS bulk ESS tail R̂
gamma 0.5116 522 654 1.005
R0 3.730 575 628 1.003
k 46.32 562 634 1.002
...The two models give nearly the same epidemic. The posterior mean of R_0 is 3.73 here against 3.81 for the Poisson model on the first page, and of \gamma 0.51 against 0.49. The dispersion k has a posterior mean of about 46, above the prior’s median of 20: at 300 boys in bed that makes the standard deviation of a count about \sqrt{1 + 300/46}\approx2.7 times the Poisson one, a modest amount of extra scatter. What differs is how much a day’s count may scatter, and the predictive check shows the consequence:
camdl fit predict @negbin
uv run plot_predictive.py ppc-both.svg Poisson=results/fits/fit_r0-*/predict/*/ \
"Negative binomial"=results/fits/fit_negbin-*/predict/*/fit predict: free_forward horizon — subsampling 200 of 800 posterior draws (raise with --n-draws)
fit predict: one_step horizon — subsampling 200 of 800 posterior draws (raise with --n-draws)
...The negative binomial’s one-day-ahead band is wider and holds 13 of the 14 counts, where the Poisson model’s held 10; its free-forward band holds all 14.
Comparing the two models
camdl compare scores each fit by how well it predicts each day’s count from the counts before it. Its headline score is the expected log predictive density,
\mathrm{elpd} = \sum_{t=0}^{13} \log\Bigl(\frac{1}{M}\sum_{m=1}^{M} \widehat p\bigl(y_t \mid y_{0:t-1}, \theta_m\bigr)\Bigr),
the total log probability the model gave each observed count before seeing it. The inner average is over M = 64 posterior draws \theta_m, each with its own particle filter, so the score is that of the model’s predictive distribution with its parameter uncertainty included, not of one parameter value. Higher is better.
camdl compare @r0 @negbincompare: fit_r0_b7a62cb1: 320 filter passes (64 draw(s) × 5 replicate(s), 1000 particles each; --draws 1 / --replicates 1 for the cheap modes)
compare: fit_negbin_2fc863ba: 320 filter passes (64 draw(s) × 5 replicate(s), 1000 particles each; --draws 1 / --replicates 1 for the cheap modes)
Model T_score elpd Δelpd LR se(Δ) evidence crps Δcrps PIT_cov90
────────────────────────────────────────────────────────────────────────────────────────────────────────
@negbin 14 -58.72 -0.32 0.724 2.17 -1.4 dB, within noise 11.442 +1.280 0.93
@r0 14 -58.39 — — — — 10.162 — 0.71
Scored steps: 14 (t0=0). Baseline: @r0.
Sorted by Δelpd ascending — best-supported model at the bottom.
Conditioning: in-sample — θ fit to all scored data; Δelpd favors the more flexible model.
ⓘ @r0: EssCollapse { step_count: 1, threshold: 100.0 }
ⓘ @r0: StartsAtPrior
ⓘ @negbin: StartsAtPrior
note: filter-noise MC SE of elpd: @r0 ±0.06, @negbin ±0.01
note: lag-1 autocorrelation of the per-step Δelpd vs baseline: @negbin +0.25
note: scores are in-sample (θ fit to all the data → optimistic in absolute level, by roughly the effective number of parameters fit, which does not cancel in a Δelpd and so biases differences toward the more flexible model) — not a leave-future-out forecast score. See LFO (#295).
Definitions and caveats: `camdl compare --explain` (or `camdl docs model-comparison`).Each row is a model. Δelpd is its elpd minus the baseline’s (the Poisson model, @r0), and se(Δ) the standard error of that difference across days. crps, the continuous ranked probability score, is an error in boys in bed, lower being better. PIT_cov90 is the fraction of days whose count fell inside the model’s one-day-ahead 90% interval.
The difference in elpd is 0.3, with a standard error of 2.2, so the counts do not prefer either model, and camdl says so (within noise). The coverage is where they differ. The Poisson model’s one-day-ahead 90% intervals contain 71% of the counts, 10 of the 14 days; the negative binomial’s contain 93%. The Poisson model’s smaller crps does not contradict this. crps rewards narrow predictions and penalizes a miss roughly in proportion to how far off it is, while the log score penalizes a confident miss much more heavily, so a model that is usually close and occasionally overconfident can win on one and not the other.
Fourteen points cannot settle the choice by score alone, and the footer says why the scores flatter: they are in-sample, because the parameters were fitted to the same counts being scored, and that favors the more flexible model. Of the two, the negative binomial is the one to report. It gives the same answer about the epidemic, and its one-day-ahead intervals contain the counts about as often as they claim to. Keep in mind how it gets there. The Poisson model’s misses fall in the decline, where the free-forward check showed the model’s epidemic ending too slowly, and the negative binomial covers those days by allowing more scatter, not by explaining why the decline was faster.
The ⓘ lines are informational. StartsAtPrior says the first day is predicted from the model’s initial state alone, before the filter has seen any count. EssCollapse says that on one day the Poisson model’s particle filter fell below 100 effective particles, a tenth of its 1000, which happens when a count is improbable under almost every simulated epidemic: one more sign of the Poisson model’s overconfidence.
What the fit cannot see
The quantities block reports, for each posterior draw, attack_rate, the fraction of the school ever infected. fit predict wrote its posterior quantiles for both fits:
head results/fits/fit_*/predict/*/quantities/attack_rate.tsv==> results/fits/fit_negbin-2fc863ba/predict/fcca52d2b336/quantities/attack_rate.tsv <==
scenario n_draws rhat ess q05 q25 q50 q75 q95
fitted 200 1.0027 222.6 0.931782 0.963303 0.972477 0.981651 0.990826
==> results/fits/fit_r0-b7a62cb1/predict/1df66e309ce8/quantities/attack_rate.tsv <==
scenario n_draws rhat ess q05 q25 q50 q75 q95
fitted 200 0.9975 236.8 0.948886 0.967235 0.976409 0.982962 0.98827Both models infer that nearly the whole school was infected: the posterior median of attack_rate (q50) is 97% for the negative binomial and 98% for the Poisson model, and even the 5% quantiles (q05) are 93% and 95%. The BMJ note counts 512 of 763 boys ill, 67%.
The fit had no way to notice. The daily counts record how many boys were in bed, never how many escaped, and the model reconciles the counts with an epidemic in which almost everyone was infected. Either about a third of the boys were infected without falling ill enough to be counted, or the model is wrong in a way the counts alone do not expose. One candidate is the simplification flagged on the first page: a boy in the infirmary was isolated, and probably infected fewer others than an I in the model does, which would also explain an epidemic that ended faster than the model’s. Telling those explanations apart needs data this fit did not use, such as a serological survey or a model with a separate isolated compartment fitted to the same counts.
This is the strongest check on the page, and the only one that uses information independent of the fourteen counts. A posterior predictive check compares the model with the data it was fitted to; an external number can contradict a model that fits those data well. The reasoning behind the fit workflow places these checks in order of strength, and the fit workflow lists what a larger fit adds.