Model structure by epidemiological outcome (#158), DAG update (#162)

[kathsherratt]

Draft. Stacked on #166, which is stacked on #164 — review those first, and this diff will shrink once they merge.

Closes #158 and #162.

What changed

The model assumed each structure predicted cases and deaths equally well. Replaced s(Method, bs = "re") with s(Method, Epi_target, bs = "re"), so a structure may predict one outcome relatively better than the other.

There is deliberately no structure main effect alongside it. mgcv's bs = "re" interaction is an unconstrained zero-mean prior over all cells, so its average across outcomes is exactly what a main effect represents. With both terms penalised, the split between them follows the relative variance estimates rather than the data — and fitted together, mgcv gave the main effect 0.001 effective degrees of freedom against 4.9 for the crossed term. Dropping it changed AIC, deviance explained and residuals by nothing on either scale.

The pooled per-structure effect is instead recovered as a contrast averaging a structure's two cells, evaluated against the full covariance matrix.

Why this matters beyond answering the issue

The old main effect was shrunk flat, so Table 2 was reporting a term carrying no information:

Structure	before	after
Agent-based	1.000 (0.994–1.007)	0.982 (0.881–1.094)
Judgement	1.000 (0.994–1.007)	0.977 (0.878–1.088)
Mechanistic	1.000 (0.994–1.007)	1.001 (0.912–1.100)
Semi-mechanistic	1.000 (0.994–1.008)	1.046 (0.946–1.155)
Statistical	1.000 (0.994–1.007)	0.995 (0.904–1.095)

Adding the crossed term improves AIC by 93 on the log scale and by 16,724 on the natural scale, where it also cuts residual skew from 4.60 to 4.38.

Findings

Human judgement models perform relatively better on cases than deaths (0.96 against 0.99), the same direction as Bosse et al. (2022). The widest separation is among agent-based models (cases 1.09, deaths 0.88) — from three models, flagged as such in text. Every interval spans the grand mean.

The direction of these contrasts is stable across error families; the magnitude is not. Under a Gaussian family the judgement separation widens and the agent-based separation disappears, because a Gaussian fit is dominated by the largest scores, which fall disproportionately among case forecasts. Recorded as a supplementary sensitivity with a pointer from the main text. Including or excluding the Hub baseline makes no material difference.

Required detour: replacing the effect extraction

gammit::extract_ranef() cannot handle a factor-by-factor random effect. It reads only the last variable name of the interaction and looks up that factor's levels, so it collects 5 labels for a 10-coefficient term and fails — taking down extraction for every term in the fit, not just the interaction. #158 could not be delivered without replacing this layer.

New R/utils-effects.R rebuilds each smooth's design matrix from the formula mgcv stores on the smooth object ($form), so labels map to coefficients exactly without assuming how mgcv orders interactions. Validated against gammit on a no-interaction fit: identical on every column to the last decimal, 48 rows. A stopifnot guards the mapping in case mgcv's construction changes.

DAG (#162)

The epidemiological outcome is added as a confounder, not merely a covariate: forecasters chose which outcomes to submit for, and that choice is associated with structure, opening a backdoor through modeller strategy.

This earned more than expected. Querying the updated diagram returns our exact covariate set as a minimal sufficient adjustment set for the direct effect, and returns no valid set for the total effect, because latent modeller strategy cannot be blocked. That turns the paper's partial-effect framing from an assertion into a derivation. The crossed term is effect modification, which a DAG does not encode, so it does not alter the adjustment set.

Prose

Corrected the discussion claim that adjusted estimates "were no different from the overall average" — that described the collapsed term; they are imprecise, not identical to the mean.

Added a discussion paragraph identifying this interaction as the one specification choice that materially changed the results, against an approach otherwise robust to covariate selection, link function and error family.

Review notes

• Model outputs under output/ are regenerated; both scales refit with no convergence warnings.
• output/diagnostics/fit-summary.csv keeps the both-terms fit alongside the final one, so the 0.001-edf claim is checkable.
• Site renders clean, both pages.
• Still open, unrelated: R/sensitivity/check-autocorrelation.R sources the deleted model-logresp.R and does not run.

🤖 Generated with Claude Code

[kathsherratt]

Added: the offsetting result, brought out through the text

Follow-up commits 9c08409 and f3da10d.

While reviewing why the shared structure effect was always flat, the underlying reason turned out to be worth stating explicitly, and it had never been written down:

Structural differences pointed in opposite directions for cases and deaths, so any term averaging over outcomes recovers close to zero.

Structure	cases	deaths	mean
Agent-based	+0.086	−0.123	−0.018
Mechanistic	−0.043	+0.046	+0.001
Statistical	−0.034	+0.024	−0.005
Judgement	−0.037	−0.008	−0.023
Semi-mechanistic	+0.028	+0.061	+0.045

Cell effects of ±0.04–0.12 average to ±0.001–0.045. A shared effect is that average, so it was near zero before any penalty applied — and a penalised term near zero is shrunk to zero exactly. This explains both why the main effect was always flat (it was flat in every earlier specification too, before the interaction existed) and why the crossed term carries 4.9 edf.

Worth noting what this is not: the signal did not migrate out of s(Model). Individual-model effects are essentially unchanged by adding the interaction (correlation 0.995, sd 0.199 → 0.194, largest change 0.049), and s(Model) is one number per model, so it cannot represent a within-model case/death difference for the 34 of 48 models that forecast both. The variance was simply unexplained before.

Where the text changed

• Discussion, prior work — the most substantive addition. Offers a mechanism for null findings in earlier structure comparisons: a pooled null is consistent with either no differences or offsetting ones, and they cannot be distinguished without letting the effect vary by target.
• Discussion, specification — why assuming a shared effect is not neutral.
• Results, pooled and per-outcome — the two possible causes of a near-mean pooled estimate; then the offsetting, illustrated by the widest-separated structure.
• Methods — the substantive cause alongside the technical one.
• Abstract — one sentence.

Deliberate limits

Every per-outcome interval still spans the grand mean, so the text says we demonstrate the mechanism rather than establish any particular contrast. The widest separation rests on three models and does not survive a Gaussian family; all three caveats are stated together rather than scattered.

The widest-separation claim is computed rather than hardcoded, so a refit cannot silently falsify it.

[kathsherratt]

Above messages generated by LLM. This sequence of changes was developed iteratively based on LLM and detailed human review.

My TLDR from a long sequence of working out what means what where:

• In this version, we set up structure to interact with the epidemiological outcome, with a fixed effect of epi outcome and the other terms remaining as they were. This specification is better at bringing out any real underlying signal from structure, which was previously getting left in the residuals.

• Because the epi outcome is specified as a fixed effect, it's not penalised, while the interaction cells are. So the fixed term absorbs the common difference between cases and deaths (the heftiest effect of any of the terms we're interested in). This leaves the interaction cells with only the structure-specific variation.

• The structure-specific interpretation is justified because the cells average to zero within each epi outcome. So they carry no component shared across structures. mgcv left all 10 coefficients (5 structures * 2 outcomes) unconstrained, so it never imposed centring, which I take as showing this is a meaningful property of the fit.

• At the same time, we drop the model-structure-only term. This was providing no distinct information (0.001 edf), and only made it more annoying to extract results from the interaction. Removing it changed no results. We now recover the pooled per-structure effect as a contrast across cells, which is now an informative pooled estimate rather than an uninformative one.

[kathsherratt]

Notes for review:

PRs are stacked on top of each other - so should be merged in order - but I've left them open hoping this helps to review the three major changes one at a time.

#164, Include the Hub baseline model

• We'd been excluding EuroCOVIDhub-baseline alongside the ensemble, but I wanted to add it back in; I think it's fine to include it as a statistical model (which the raters all classified it as).

#166, Change the error family

• Compared Gaussian, Gamma and Tweedie. Tweedie fits decently, Gamma fits just as well but doesn't converge. Detail in the PR body.

#167, Let effect of model structure differ by outcome

• Add interaction term between model structure and epi outcome. Also updates the DAG. See PR comments for detail of what changed here, but main impact shown in manuscript Table 2. Where the main estimate used to be a flat 1.00 for all five model types, specifying the interaction (and dropping the isolated term for model structure) shows a more interpretable result. Still overlapping 1 for all, but I think this actually means something now rather than just a term that shrunk to nothing.
• The manuscript is updated to reflect this, just currently rather peppered with LLM speak. Still worth looking at as the single source of truth across all these changes (but I'll also take another pass over the text).