Observation models

Choosing how your model meets your data

Estimated time: 55 minutes | Requires: MCMC

Slides

Introduction

Every likelihood in this course so far has assumed that observed case counts follow a Poisson distribution:

$$\text{obs}_t \sim \text{Poisson}(\rho \cdot I_t)$$

That is a modelling choice, and it makes a specific claim about how observation works. This session asks whether the claim holds, what happens when it does not, and what a likelihood is really doing when you fit with one.

Objectives

In this session you will:

1. Critically evaluate observation model choices (Poisson, negative binomial)
2. Check whether the Poisson variance assumption holds for real data
3. Understand negative log-likelihood as an implicit distance metric
4. Make informed choices about observation models for your own research

Setup

Source file

The source file of this session is located at sessions/observation_models.qmd.

Load libraries

using Distributions
using DataFrames
using CSV
using Plots
using StatsPlots
using Turing
using FlexiChains ## chain objects and diagnostics
using FlexiChains: @varname ## refer to a model variable by name
using Random
using LinearAlgebra
using DrWatson
using StatsBase  # For Weights in ABC-SMC
using MFIIDD ## the course package, which holds the SIR implementation

Initialisation

We set a random seed for reproducibility.

Random.seed!(1234)

TaskLocalRNG()

Load Tristan da Cunha data

We’ll use the same influenza outbreak data from 1971:

# Load data
flu_tdc = CSV.read(datadir("flu_tdc_1971.csv"), DataFrame)

# Plot
scatter(
    flu_tdc.time, flu_tdc.obs,
    xlabel = "Time (days)", ylabel = "Daily incidence",
    label = "Observed cases", markersize = 4, color = :red,
    title = "Tristan da Cunha 1971 Influenza Outbreak",
    legend = :topright
)

We have 59 daily observations of case counts.

The model we are fitting

The transmission model is the SIR model from the earlier sessions, and $\rho$, the fraction of infections that are reported, is the one the model checking session fitted.

These two functions live in src/sir_model.jl and arrived with using MFIIDD above.

function sir_ode!(du, u, p, t)
    S, I, R = u
    R_0, D_inf = p
    N = S + I + R

    β = R_0 / D_inf
    ν = 1 / D_inf

    du[1] = -β * S * I / N          # dS/dt
    du[2] = β * S * I / N - ν * I   # dI/dt
    du[3] = ν * I                   # dR/dt
end

function simulate_sir(θ, init_state, times)
    u0 = [init_state[:S], init_state[:I], init_state[:R]]
    prob = remake(
        SIR_BASE_PROBLEM;
        u0 = u0,
        p = [θ[:R_0], θ[:D_inf]],
        tspan = (times[1], times[end]),
    )
    sol = solve(prob, Tsit5(), saveat = times)

    df = DataFrame(
        time = sol.t,
        S = [u[1] for u in sol.u],
        I = [u[2] for u in sol.u],
        R = [u[3] for u in sol.u],
    )

    # Daily incidence: new infections, the flow out of S
    df.Inc = [0.0; -diff(df.S)]

    return df
end

## `times` starts at 0.0, the left edge of the first observation day, because the
## incidence over that day is S(0) - S(1) and there is no way to form it without a
## point at 0. `traj.Inc[i + 1]` is then the incidence over the day the i-th
## observation counts
times = [0.0; flu_tdc.time]

## The island as the ship lands, at the start of day 1. Of the five islanders it
## brought back, three had been ill during the eight day voyage and are past their
## infectious period, so they start recovered; the two who fell ill on landing
## start infectious. N = 284, the island population
init_state = Dict(:S => 279.0, :I => 2.0, :R => 3.0)

Dict{Symbol, Float64} with 3 entries:
  :I => 2.0
  :R => 3.0
  :S => 279.0

Each fit below is checked against the data with a posterior predictive interval.

The Poisson likelihood

What we have used so far

This is what we’ve been doing throughout the course:

@model function sir_poisson(times, init_state, n_obs)
    # Priors
    R_0 ~ Uniform(1.0, 20.0)
    D_inf ~ Uniform(1.0, 10.0)
    ρ ~ Uniform(0.1, 1.0)  # Reporting rate

    # Simulate
    θ = Dict(:R_0 => R_0, :D_inf => D_inf)
    traj = simulate_sir(θ, init_state, times)

    # Likelihood: Poisson observation model
    lambdas = [max(ρ * traj.Inc[i + 1], 1e-10) for i in 1:n_obs]
    obs ~ arraydist(Poisson.(lambdas))

    return traj
end

sir_poisson (generic function with 2 methods)

NoteReminder: arraydist and | conditioning

obs ~ arraydist(Poisson.(lambdas)) declares a vector of independent Poisson observations, one per time point. The | operator, applied when creating the model instance below, conditions the model on the actual data. See the MCMC session for a full explanation of this pattern.

We fit this model:

# Fit with NUTS
model_poisson = sir_poisson(times, init_state, length(flu_tdc.obs)) | (; obs = flu_tdc.obs)
chain_poisson = sample(model_poisson, NUTS(0.65), 500, progress=false)
# The 0.65 is the target acceptance rate for NUTS — higher values (closer to 1)
# make the sampler more conservative, lower values allow larger steps.
# The default 0.65 works well for most models.

┌ Info: Found initial step size
└   ϵ = 0.2

╭─FlexiChain (500 iterations, 1 chain) ────────────────────────────────────────╮
│ ↓ iter  = 251:750                                                            │
│ → chain = 1:1                                                                │
│                                                                              │
│ Parameters (3) ── AbstractPPL.VarName                                        │
│  Float64  R_0, D_inf, ρ                                                      │
│                                                                              │
│ Extras (14)                                                                  │
│  Int64    n_steps, tree_depth                                                │
│  Bool     is_accept, numerical_error                                         │
│  Float64  acceptance_rate, log_density, hamiltonian_energy,                  │
│           hamiltonian_energy_error, max_hamiltonian_energy_error, step_size, │
│           nom_step_size, logprior, loglikelihood, logjoint                   │
╰──────────────────────────────────────────────────────────────────────────────╯

# Plot posteriors
plot(chain_poisson[[:R_0, :D_inf, :ρ]], size=(900, 600))

# Summary statistics
summarystats(chain_poisson)

╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (3) ── AbstractPPL.VarName                                        │
│  Float64  R_0, D_inf, ρ                                                      │
│                                                                              │
│ Extras (14)                                                                  │
│  Float64  n_steps, is_accept, acceptance_rate, log_density,                  │
│           hamiltonian_energy, hamiltonian_energy_error,                      │
│           max_hamiltonian_energy_error, tree_depth, numerical_error,         │
│           step_size, nom_step_size, logprior, loglikelihood, logjoint        │
│                                                                              │
│ Summary                                                                      │
│   param    mean     std    mcse  ess_bulk  ess_tail    rhat      q5  …       │
│     R_0  1.9110  0.2163  0.0254   93.1638   69.4264  1.0035  1.6392  …       │
│   D_inf  3.3647  0.4999  0.0565   92.1495   66.6760  1.0051  2.6957  …       │
│       ρ  0.9900  0.0111  0.0010  186.3279  198.8131  1.0030  0.9689  …       │
╰──────────────────────────────────────────────────────────────────────────────╯

Does it reproduce the data?

Each fit is checked the same way.

function plot_predictive(model, chain, title)
    ## Draw replicate observations from the fitted model, one set per posterior draw
    predictions = predict(decondition(model), chain)
    draws = predictions[@varname(obs)]
    per_day = [[draws[i, 1][t] for i in axes(draws, 1)] for t in 1:length(flu_tdc.obs)]

    lower = [quantile(day, 0.025) for day in per_day]
    middle = [quantile(day, 0.5) for day in per_day]
    upper = [quantile(day, 0.975) for day in per_day]

    p = plot(flu_tdc.time, middle,
             ribbon = (middle .- lower, upper .- middle),
             linewidth = 2, color = :blue, fillalpha = 0.3,
             label = "95% predictive interval",
             xlabel = "Time (days)", ylabel = "Daily cases", title = title)
    scatter!(p, flu_tdc.time, flu_tdc.obs,
             label = "Observed", markersize = 3, color = :red)

    covered = count(lower .<= flu_tdc.obs .<= upper)
    return p, covered
end

p_poisson, covered_poisson = plot_predictive(model_poisson, chain_poisson, "Poisson")
p_poisson

The interval covers 35 of the 59 days. A 95% interval should cover about 56.

What does Poisson assume?

The Poisson distribution has a key property: variance = mean

# Illustrate Poisson assumption
λ_values = [5, 20, 50]
p = plot(title="Poisson Distribution: Variance = Mean",
         xlabel="Count", ylabel="Probability")

for λ in λ_values
    dist = Poisson(λ)
    x_range = max(0, λ-20):min(λ+20, λ*2)
    plot!(p, x_range, [pdf(dist, x) for x in x_range],
          label="λ=$λ (σ²=$λ)", linewidth=2)
end
p

TipExercise: Does the assumption hold for our data?

The Poisson model claims the variance of the observations equals their mean. Check whether that is plausible for the Tristan da Cunha data.

1. Fit the Poisson model above and extract the posterior mean trajectory.
2. Compute the residuals (observed minus predicted) and their variance.
3. Compare that variance against the mean prediction. If the data were Poisson, the two should be roughly equal.

What do you find, and what would it imply for the parameter estimates?

# Residual variance against mean prediction, under the Poisson fit
θ_fit = Dict(:R_0 => mean(chain_poisson[:R_0]), :D_inf => mean(chain_poisson[:D_inf]))
pred_fit = mean(chain_poisson[:ρ]) .* simulate_sir(θ_fit, init_state, times).Inc[2:end]
resid_fit = flu_tdc.obs .- pred_fit

println("mean prediction:    ", round(mean(pred_fit), digits = 2))
println("residual variance:  ", round(var(resid_fit), digits = 1))
println("ratio:              ", round(var(resid_fit) / mean(pred_fit), digits = 1),
        "  (Poisson demands 1)")

mean prediction:    3.58
residual variance:  87.0
ratio:              24.3  (Poisson demands 1)

TipSolution

The ratio printed above is around 25, where Poisson demands 1. The data are wildly overdispersed relative to Poisson.

Two things produce that, and it is worth keeping them apart. Some of it is genuine observation noise: two islanders infected in the same household are not independent draws, and reporting was not uniform across the outbreak. Most of it here is something else, namely that a single-wave SIR fitted to a two-wave outbreak is wrong about the mean, and everything it gets wrong about the mean lands in the residuals.

Either way the consequence for the Poisson fit is the same. It treats every departure from the trajectory as strong evidence against it, so the posterior comes out far too narrow. Nor can it buy its way out by moving the parameters: the four days around the observed peak are underpredicted by factors of about four to nine and stay that way, because the many low-count days of the second wave outweigh them. The negative binomial in the next section loosens that grip, but notice that it does so by explaining the misfit as noise rather than by fixing the model.

Reasons the assumption might fail for epidemic data:

• Under-reporting might be systematic rather than random
• Behavioural changes might create clustering
• Surveillance quality might vary over time

TipExercise: Predict what the negative binomial will change

Before you fit it, commit to an answer for each of these. Write them down — the point is to be caught out if you are wrong.

1. Will the posterior for $R_0$ get wider, narrower, or stay about the same?
2. Will the posterior mean of $R_0$ move much?
3. The model adds one parameter, $\phi$. If the data were genuinely Poisson, what value would you expect the posterior for $\phi$ to favour?

Then fit it below and check.

TipSolution

1. Wider, and by a lot. The negative binomial can attribute a departure from the trajectory to observation noise, so the data pin the parameters down less tightly. Here the posterior standard deviation for $R_0$ roughly doubles. Loosening the observation model transfers uncertainty back to the parameters, which is the honest direction: the Poisson intervals were narrow because the model was forbidden from admitting the data were noisy.

2. The mean moves too, and that is the surprise. $R_0$ goes from about 1.9 to about 2.7, several Poisson standard deviations apart. The two models do not merely disagree about the uncertainty, they disagree about the answer.

   The misfit sits almost entirely in the first wave. Days 6 to 9 carry about three quarters of the squared residual, and day 7 alone about a third: 47 cases reported against about five predicted. The fitted trajectory accounts for about 211 of the 312 reported cases, and two thirds of that shortfall falls in the first eighteen days, which carry about 90% of the squared residual. The Poisson fit cannot lift its peak to meet those days, because the outbreak has a second wave and a single SIR curve has only one peak to spend. Days 19 to 59 carry 109 of the 312 cases and rise again to 12 on day 32. Raising $R_0$ makes the epidemic burn out sooner: at $R_0 = 7$, holding $D_\mathrm{inf}$ at its fitted value, the trajectory peaks on day 3 and predicts almost no cases at all after day 19, where the Poisson likelihood then gives those 109 observed cases essentially no probability. Fitting the same model to the first eighteen days alone removes the second wave, and $R_0$ comes out around 7.

   The negative binomial adds a dispersion parameter, so departures from its mean cost much less: over days 19 to 59 its log-likelihood is about 54 units better than the Poisson’s. Less bound to track every day exactly, it settles on a longer infectious period and a larger $R_0$, and so on a larger predicted total, about 245 cases against 211.

   Look at $\rho$ in the comparison table below. The Poisson fit pushes it to about 0.99, against the upper bound of its Uniform(0.1, 1.0) prior. A parameter pinned to the edge of its prior is a model straining against a constraint, and here it is straining for the same reason. Even reporting everybody would not be enough: the model starts with 279 susceptibles, so it can never produce more than 279 cases against the 312 observed. The reporting rate is pinned at a bound that still cannot close the gap, which is about as clear a sign as this session offers that the shortfall is structural.

3. Large $\phi$ is the right answer to the question as asked, because variance is $\mu + \mu^2/\phi$ and the second term vanishes as $\phi \to \infty$. The data say the opposite: the posterior for $\phi$ sits near 0.5, so the variance is roughly $\mu + 2\mu^2$ and the extra flexibility was very much needed.

   Note that a small $\phi$ is not proof of messy reporting. A single-wave SIR cannot reproduce a two-wave outbreak, and the negative binomial is absorbing that structural failure as though it were noise. An observation model will always volunteer to explain a misspecified mean, which is why the model checking session comes first.

The negative binomial likelihood

What does the negative binomial assume?

The negative binomial keeps the same mean and adds a second parameter, so the variance is free to exceed it:

$$\text{obs}_t \sim \text{NegativeBinomial}(\mu_t, \phi), \qquad \text{Var}(\text{obs}_t) = \mu_t + \frac{\mu_t^2}{\phi}$$

where $\mu_t = \rho \cdot I_t$ and $\phi$ controls the overdispersion. Large $\phi$ recovers the Poisson, and small $\phi$ gives a much wider spread at the same mean.

# The same mean, three observation models
counts = 0:80
p = plot(title="Same mean of 20, different spreads",
         xlabel="Count", ylabel="Probability")

plot!(p, counts, pdf.(Poisson(20), counts),
      label="Poisson (σ²=20)", linewidth=2)

for φ in [5.0, 1.0]
    ## Julia parameterises the negative binomial by (r, p), so a mean μ with
    ## overdispersion φ is r=φ and p=φ/(φ+μ)
    dist = NegativeBinomial(φ, φ/(φ+20.0))
    plot!(p, counts, pdf.(dist, counts),
          label="NegBin φ=$(Int(φ)) (σ²=$(round(Int, var(dist))))", linewidth=2)
end

vline!(p, [20], color=:grey, linestyle=:dot, linewidth=2, label="Mean")
p

The extra variance grows with the square of the mean, so a day with a high predicted count is allowed to fall further from it than a day with a low one.

Everything else the Poisson assumed still holds. Observations are conditionally independent given the trajectory, and the same fraction of infections is reported on every day. Overdispersion is the only assumption the negative binomial relaxes.

TipExercise: what would you expect φ to be?

Look back at the ratio of residual variance to mean prediction you computed for the Poisson fit. If that ratio came entirely from overdispersion, and the mean prediction were about 5 cases a day, roughly what $\phi$ would reproduce it?

Use $\text{Var} = \mu + \mu^2/\phi$ and solve for $\phi$.

TipSolution

With $\mu = 5$ and a variance 25 times the mean, $125 = 5 + 25/\phi$, so $\phi \approx 0.2$. The fit below returns a $\phi$ of about 0.5, in the same region.

The calculation assumes the whole ratio is overdispersion, which it is not: the mean itself is wrong, and the fitted $\phi$ is absorbing both.

Fitting it

@model function sir_negbin(times, init_state, n_obs)
    # Priors
    R_0 ~ Uniform(1.0, 20.0)
    D_inf ~ Uniform(1.0, 10.0)
    ρ ~ Uniform(0.1, 1.0)
    φ ~ Exponential(10.0)  # Overdispersion parameter

    # Simulate
    θ = Dict(:R_0 => R_0, :D_inf => D_inf)
    traj = simulate_sir(θ, init_state, times)

    # Likelihood: Negative Binomial observation model
    mus = [max(ρ * traj.Inc[i + 1], 1e-10) for i in 1:n_obs]
    # NegativeBinomial(r, p) where r=φ, p=φ/(φ+μ)
    obs ~ arraydist([NegativeBinomial(φ, φ/(φ+μ)) for μ in mus])

    return traj
end

sir_negbin (generic function with 2 methods)

TipCode explainer

• NegativeBinomial(r, p) - Julia’s parameterisation uses r (number of successes) and p (success probability). To get mean μ with overdispersion φ, use r=φ and p=φ/(φ+μ). This gives variance = μ + μ²/φ, so larger φ approaches Poisson, smaller φ gives more overdispersion.

Now fit it, so that the predictions you just committed to can be checked.

model_negbin = sir_negbin(times, init_state, length(flu_tdc.obs)) | (; obs = flu_tdc.obs)
chain_negbin = sample(model_negbin, NUTS(0.65), 500, progress=false)
summarystats(chain_negbin)

┌ Info: Found initial step size
└   ϵ = 0.2

╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (4) ── AbstractPPL.VarName                                        │
│  Float64  R_0, D_inf, ρ, φ                                                   │
│                                                                              │
│ Extras (14)                                                                  │
│  Float64  n_steps, is_accept, acceptance_rate, log_density,                  │
│           hamiltonian_energy, hamiltonian_energy_error,                      │
│           max_hamiltonian_energy_error, tree_depth, numerical_error,         │
│           step_size, nom_step_size, logprior, loglikelihood, logjoint        │
│                                                                              │
│ Summary                                                                      │
│   param    mean     std    mcse  ess_bulk  ess_tail    rhat      q5  …       │
│     R_0  2.6644  0.5189  0.0385  177.8099  177.8830  0.9988  1.7886  …       │
│   D_inf  6.6800  1.7593  0.1319  171.7470  199.0966  0.9999  3.8495  …       │
│       ρ  0.9631  0.0377  0.0029   99.5079  235.2553  1.0015  0.8812  …       │
│       φ  0.5005  0.1133  0.0064  324.4816  227.8952  0.9992  0.3337  …       │
╰──────────────────────────────────────────────────────────────────────────────╯

p_negbin, covered_negbin = plot_predictive(model_negbin, chain_negbin,
                                           "Negative binomial")
p_negbin

The interval now covers 54 days against the Poisson’s 35. It does so by being wide rather than by following the data: the trajectory underneath it still peaks in the wrong place.

# Compare the two observation models on the parameter they disagree about most
comparison = DataFrame(
    parameter = ["R_0", "D_inf", "ρ"],
    poisson_mean = [mean(chain_poisson[p]) for p in (:R_0, :D_inf, :ρ)],
    poisson_sd = [std(chain_poisson[p]) for p in (:R_0, :D_inf, :ρ)],
    negbin_mean = [mean(chain_negbin[p]) for p in (:R_0, :D_inf, :ρ)],
    negbin_sd = [std(chain_negbin[p]) for p in (:R_0, :D_inf, :ρ)],
)
comparison.sd_ratio = round.(comparison.negbin_sd ./ comparison.poisson_sd, digits = 1)
comparison

3×6 DataFrame

Row	parameter	poisson_mean	poisson_sd	negbin_mean	negbin_sd	sd_ratio
	String	Float64	Float64	Float64	Float64	Float64
1	R_0	1.91101	0.216349	2.66437	0.518879	2.4
2	D_inf	3.36472	0.49987	6.67999	1.75933	3.5
3	ρ	0.99004	0.0110836	0.963125	0.0377182	3.4

WarningCheck the chain before trusting the comparison

Both fits here are 500 NUTS draws, which is enough to see the difference in question but not enough to report. The negative binomial chain in particular comes back with Monte Carlo standard errors of about 0.17 on $D_\mathrm{inf}$ and 0.05 on $R_0$, so the first decimal of one and the second of the other are already noise, however many digits the table prints. Those two are the parameters the comparison is about. The diagnostics session is the reason you should have looked at summarystats before reading the table, rather than after.

# What did the data say about the overdispersion parameter?
φ_post = vec(chain_negbin[:φ])
println("φ posterior: median $(round(median(φ_post), digits = 2)), ",
        "95% interval [$(round(quantile(φ_post, 0.025), digits = 2)), ",
        "$(round(quantile(φ_post, 0.975), digits = 2))]")

φ posterior: median 0.49, 95% interval [0.31, 0.73]

Why the negative binomial?

A negative binomial is a Poisson whose rate varies from day to day with a gamma distribution, so it is what you get when the reporting rate is not constant and you have not modelled why. Clustered transmission produces the same shape: household outbreaks and superspreading events deliver several cases at once, so counts vary more than independent arrivals would.

Two weaker reasons also apply. It fits these data better than the Poisson, and it is cheap to evaluate inside MCMC.

What it does not do is tell you which of those mechanisms was at work, or distinguish any of them from a mean that is simply wrong.

The two fits differ in one thing: what they charge the model for a day that misses the prediction. That charge is what moved $R_0$, so it is worth looking at on its own.

Interlude: likelihood as distance

Likelihood as distance

The negative log-likelihood is actually a distance function between observed and predicted data:

$$\text{Distance} = -\log p(\text{obs} | \text{predicted})$$

We can see what distance each likelihood implies:

# Compare different observation models as distance functions
obs_val = 10.0
predicted_range = 1.0:0.5:20.0

# Poisson: -log P(obs | λ)
poisson_distance = [-logpdf(Poisson(λ), obs_val) for λ in predicted_range]

# Negative Binomial with different φ (convert to r,p parametrization)
# NegativeBinomial(r, p) where r=φ, p=φ/(φ+μ)
negbin_distance_phi5 = [-logpdf(NegativeBinomial(5.0, 5.0/(5.0+λ)), obs_val) for λ in predicted_range]
negbin_distance_phi20 = [-logpdf(NegativeBinomial(20.0, 20.0/(20.0+λ)), obs_val) for λ in predicted_range]

# Gaussian (for comparison)
gaussian_distance = [-logpdf(Normal(λ, 3.0), obs_val) for λ in predicted_range]

# Plot
plot(predicted_range, poisson_distance,
     label="Poisson", linewidth=2,
     xlabel="Predicted value", ylabel="-log(likelihood)",
     title="Different Likelihoods as Distance Functions (obs=10)",
     legend=:topright)
plot!(predicted_range, negbin_distance_phi5,
      label="NegBin(φ=5)", linewidth=2)
plot!(predicted_range, negbin_distance_phi20,
      label="NegBin(φ=20)", linewidth=2)
plot!(predicted_range, gaussian_distance,
      label="Gaussian(σ=3)", linewidth=2, linestyle=:dash)

Each curve is at its lowest where the prediction meets the observation and rises as the two part company. What differs is how steeply, and on which side:

• Poisson: penalises deviations relative to the predicted value rather than in absolute cases, so the same gap costs more at a low prediction than at a high one
• Negative binomial: the same shape with a much flatter upper tail, so a prediction well below the observation costs little
• Gaussian: penalises the absolute gap symmetrically, at a cost that does not depend on the size of the counts

That flat tail is what moved the fit. Day 7 of the outbreak reports 47 cases where the Poisson fit predicts about 5:

# What each observation model charges for day 7
obs_day7, pred_day7 = 47, 5.0
φ_fit = median(φ_post)

println("Poisson:  ", round(-logpdf(Poisson(pred_day7), obs_day7), digits = 1))
println("NegBin:   ", round(-logpdf(NegativeBinomial(φ_fit, φ_fit/(φ_fit+pred_day7)),
                                    obs_day7), digits = 1))

Poisson:  66.2
NegBin:   8.1

The Poisson is charged about eight times as much for that single day as the negative binomial. That is the shift in the comparison table above, seen from the observation model rather than from the chain: a likelihood that charges less for a badly missed day leaves the parameters freer to go elsewhere.

Learning points

TipLearning points

• An observation model is a claim about how cases reach the data. The common choices have a mechanism behind them, and a negative binomial is a Poisson whose rate varies, which is what clustered transmission and uneven reporting both produce
• The Poisson assumption that variance equals the mean is testable, and surveillance data routinely breaks it
• A negative log-likelihood is a distance function, so choosing an observation model is choosing how to measure the gap between prediction and data
• An observation model that is too tight makes the posterior overconfident. Loosening it widens the intervals and moves the mean, but a mean that is wrong for structural reasons is not fixed by the observation model
• Some problems look like observation-model problems and are not. A delay between an event and its report, or a reporting fraction that changes part-way through an outbreak, needs a change to the model rather than a wider likelihood

TipExercise: Which observation model for which system?

The choice is not abstract. For each of these surveillance systems, say which observation model you would start with and what feature of the system drives the answer.

1. A sentinel network of GP practices reporting influenza-like illness weekly, with stable participation and no weekend effect.
2. Notifications of a rare disease, where a single case triggers an investigation that often finds several more in the same household.
3. Deaths in a national registry, complete but reported with a delay that varies between one and three weeks.
4. A new pathogen where testing capacity tripled part-way through the outbreak.

TipSolution

1. Poisson is defensible. Stable reporting, no obvious clustering, and cases arriving independently is roughly what Poisson assumes. Check the variance against the mean anyway rather than assuming.

2. Negative binomial. Case-finding around an index case produces exactly the clustering that breaks the independence assumption, and the visible symptom is variance far above the mean.

3. Neither, on its own. The problem here is not the spread but the timing: a case counted this week may belong to three weeks ago. No observation model fixes that. You need a delay distribution between the event and its report, which is a change to the model rather than to the likelihood family.

4. Neither, again. Tripling testing changes the proportion of infections you observe partway through, so a constant reporting fraction $\rho$ is wrong. You need $\rho$ to vary with time. Fitting a fixed $\rho$ here will make the model infer a change in transmission that never happened.

The general point: 1 and 2 are observation-model choices, 3 and 4 look like them and are not. Asking “is this about how much noise, or about what the data mean?” separates the two.

Going further

Other observation models

The choice of observation model affects inference more than many researchers realise. Consider exploring:

• Gaussian likelihood: Appropriate when counts are large and approximately symmetric
• Zero-inflated models: When surveillance systems miss entire days
• Reporting delay models: When cases are reported with variable lag

Computing the variance/mean ratio in your data can help diagnose whether Poisson (ratio ≈ 1) or Negative Binomial (ratio > 1) is more appropriate.

Next session

If a likelihood is really a distance function that you committed to in advance by naming a distribution, an obvious question follows: what if we chose the distance directly, and never wrote a likelihood at all?

That is Approximate Bayesian computation, tomorrow morning. Worth turning over tonight: which features of these data would you choose to match?