Reinstate baseline mode; preserve alternative model fitting diagnostics

.gitignore

@@ -76,5 +76,8 @@ README_files/
 # supplement embeds those, so they stay tracked.
 /output/**/plots/*.pdf
 
+# Stray default device output from non-interactive Rscript runs
+/Rplots.pdf
+
 # Tool caches
 .panache-cache/

CLAUDE.md

@@ -11,7 +11,7 @@ However, when many forecasters each predict many different targets, it becomes d
 Here we use a regression model to separate the effect of the forecasting method, from the difficulty of the target, in forecast performance.
 
 We evaluated forecasts of weekly COVID-19 cases and deaths over two years across 32 European countries, scoring them against observed data with the Weighted Interval Score (WIS).
-We expected a model's structure to shape how well it predicted, so we classified 47 models by structure (agent-based, mechanistic, semi-mechanistic, statistical, or human judgement) and estimated how much structure alone affected performance.
+We expected a model's structure to shape how well it predicted, so we classified 48 models by structure (agent-based, mechanistic, semi-mechanistic, statistical, or human judgement) and estimated how much structure alone affected performance.
 A generalised additive mixed model let us adjust for everything that makes a target easier or harder to predict: the outcome being forecast, its level and trend, the dominant variant, the country, the forecast horizon, and differences between individual models.
 
 Once we accounted for the difficulty of the target, no single type of model performed best.
@@ -39,7 +39,8 @@ As infectious disease forecasting grows, we encourage evaluators to choose from
   - Isolates impact of Method (model structure) and CountryTargets (geographic specificity)
   - Uses `mgcv`, `gammit`, and `gratia` packages
   - Outputs (per scale, under scale-named subdirs of `output/` — `log/`, `natural/`): `results.rds` (includes fitted `data`), `fit_obs.rds`, and diagnostic plots (`plots/check_joint.png`)
-  - Defines `model_wis(scoring_scale, family_link = "log", output_dir)`; sourcing alone fits nothing. Call it once per scale (`log`, `natural`) to write outputs. Must be run separately before rendering — it is **not** sourced by `report/quarto/_results.qmd`, which only reads `output/log/results.rds`
+  - Defines `model_wis(scoring_scale, family_link = "log", output_dir, spec_label)`; sourcing alone fits nothing. Call it once per scale (`log`, `natural`) to write outputs. Must be run separately before rendering — it is **not** sourced by `report/quarto/_results.qmd`, which only reads `output/log/results.rds`
+  - `spec_label` archives that fit's `appraise()` panel and residual/fit statistics under `output/diagnostics/`, upserting a row in `fit-summary.csv` keyed on (`spec_label`, `scale`) so successive model specifications stay comparable. Pass a new label whenever the specification changes.
 
 - **analysis-descriptive.R**: Descriptive statistics and summary tables
   - Bootstrap confidence intervals
@@ -89,10 +90,10 @@ The manuscript prose lives in per-section Quarto files under `report/quarto/`, a
 
 ### Rendered analysis (code and outputs)
 
-- `report/manuscript.qmd` — top-level Quarto document; includes the `report/quarto/_*.qmd` sections and the supplement. This is the render target.
+- `report/manuscript.qmd` — top-level Quarto document; includes the `report/quarto/_*.qmd` sections. The supplement is a separate page, not included here. This is the render target for the manuscript alone.
 - `report/quarto/_results.qmd` — results section; sources R scripts and renders figures/tables.
-- `report/quarto/supplement/_supplement.qmd` — supplementary materials (with `future-work.qmd` and figure assets alongside it).
-- Root render wrappers (`quarto render` uses `_quarto.yml` → `_site/`): `index.qmd` wraps the manuscript, `supplement.qmd` wraps the supplement; two-page site with navbar. Bibliography `report/references.bib`, style `report/plos-computational-biology.csl`.
+- `report/supplement.qmd` — supplementary materials, including its own setup chunk and sensitivity analyses.
+- Site build (`quarto render` uses `_quarto.yml` → `_site/`): renders `index.qmd` (a thin wrapper including `report/manuscript.qmd`) and `report/supplement.qmd` directly; two-page site with navbar. Bibliography `report/references.bib`, style `report/plos-computational-biology.csl`.
 - Pre-print: [medRxiv 10.1101/2025.04.10.25325611](https://doi.org/10.1101/2025.04.10.25325611)
 
 **Note**: manuscript prose and rendered analysis are separate. The section `.qmd` files are not auto-generated — changes to analysis code and changes to manuscript text must be coordinated manually.
@@ -121,14 +122,16 @@ source(here("R", "process-data.R"))
 # 3. Fit GAMM to weighted interval scores (run before rendering).
 #    Sourcing only DEFINES model_wis(); call it per scale to write output/<scale>/.
 source(here("R", "analysis-model.R"))
-model_wis(scoring_scale = "log",     output_dir = here("output", "log"))
-model_wis(scoring_scale = "natural", output_dir = here("output", "natural"))
+model_wis(scoring_scale = "log",     output_dir = here("output", "log"),
+          spec_label = "baseline-included-gaussian-log")
+model_wis(scoring_scale = "natural", output_dir = here("output", "natural"),
+          spec_label = "baseline-included-gaussian-log")
 
 # 4. Render the manuscript alone (results section only; supplement is a separate page)
 # quarto::quarto_render("report/manuscript.qmd")
 
 # Or build the full two-page website (manuscript + supplement, with navbar):
-# quarto render   # uses root _quarto.yml; index.qmd + supplement.qmd wrap the report/ content
+# quarto render   # uses root _quarto.yml; renders index.qmd and report/supplement.qmd
 ```
 
 ## Making Changes
@@ -138,7 +141,7 @@ model_wis(scoring_scale = "natural", output_dir = here("output", "natural"))
 | Change manuscript prose (wording, framing, conclusions) | Relevant `report/quarto/_*.qmd` section file |
 | Change analysis, model, or figures | Relevant `R/` script. `_results.qmd` sources `process-data.R`, `analysis-descriptive.R`, `plot-model-results.R` at render. But `analysis-model.R` and `plot-model-flow.R` are **not** sourced — re-run `model_wis()` per scale (and regenerate the flowchart) to refresh `output/` before rendering |
 | Respond to a reviewer comment | Check `submission/Revision_reviews-response.md`, update `R/` script if needed, then update the relevant `report/quarto/_*.qmd`, mark as completed in `submission/Revision_reviews-response.md`, and close the relevant Github Issue with a note |
-| Add or change a supplementary figure | Relevant `R/` script + `report/quarto/supplement/_supplement.qmd` |
+| Add or change a supplementary figure | Relevant `R/` script + `report/supplement.qmd` |
 | All changes | Update `NEWS.md` (change log; newest first) |
 
 ## Dependencies
@@ -157,6 +160,8 @@ Major R packages:
 
 Outstanding issues.
 Status: [ ] not started, [x] done.
+- [ ] Update manuscript text to clarify: the random-effect variance penalised any effect from model structure to zero, i.e. the joint model is reporting no information about structure - use this phrasing instead of "overlapping uncertainty"
+- [ ] Untrack the Quarto freeze cache on branch `supplementary-descriptive`. Commit 8b7ede2 is the first to track `_freeze/` (10.1 MB: binaries, minified JS, vendored `site_libs` like `jquery-3.5.1/`), and it contains a stale duplicate of the same 12 figures under `_freeze/report/descriptive-scores/` alongside the current `_freeze/report/quarto/_supplementary-descriptive/`. CI (`.github/workflows/render-report.yaml`) installs R via renv and runs `quarto render` from source, so the cache is not load-bearing — `main` tracks none of it. Fix: `git rm -r --cached _freeze`, add `/_freeze/` to `.gitignore` next to the other Quarto render artefacts, note in `NEWS.md`.
 
 ### Verification
 After changing `analysis-model.R` (or upstream scoring/data), regenerate the saved model outputs first — the manuscript reads `output/log/results.rds` and `fit_obs.rds`, it does not re-fit. Stale outputs render silently wrong, or break (e.g. the supplement density chunk needs `results$data`). Then render `report/manuscript.qmd` (or `quarto render` for the full site) and check figures render correctly.

NEWS.md

@@ -3,6 +3,27 @@
 Notable changes to the analysis, manuscript, and repository.
 Newest first.
 
+## Unreleased — Include the Hub baseline model; archive fit diagnostics per specification
+
+`R/analysis-model.R`, `R/plot-model-flow.R`, `R/sensitivity/check-autocorrelation.R`, `R/sensitivity/check-link-robustness.R`, `report/quarto/_abstract.qmd`, `report/quarto/_methods.qmd`, `report/quarto/_results.qmd`, `report/quarto/_discussion.qmd`, `report/supplement.qmd`, `CLAUDE.md`
+
+`EuroCOVIDhub-baseline` is now included in the analysis; only `EuroCOVIDhub-ensemble` remains excluded.
+The two had been dropped together by a single `grepl("EuroCOVIDhub-")` filter, but they are not equivalent: the ensemble is a function of the contributed forecasts and would double-count them, whereas the baseline is an independently specified statistical model and belongs in the sample on the same terms as any other participant.
+Narrowed the filter at every remaining site, and relabelled the flow-diagram exclusion step from "Not created by Hub" to "Not the Hub ensemble".
+
+The sample grows from 47 to 48 models across 38 teams, with the statistical structure group going from 12 to 13 models.
+Because the baseline submitted for every country in almost every week, the statistical group's share of forecasts rises to roughly 40%, so the previous claim that mechanistic, semi-mechanistic and statistical models each contributed about a third no longer holds and has been replaced.
+Model counts, rater-disagreement counts, and per-structure forecast shares in the results text are now computed inline from the data rather than hardcoded, so they track future changes to the sample.
+Added a sentence to the Discussion noting that the baseline anchors the statistical group towards the performance achievable without epidemiological structure.
+
+`model_wis()` gains a `spec_label` argument. When supplied it writes an archived copy of the `appraise()` panel plus a row of fit statistics (family, link, formula, n, AIC, deviance explained, deviance-residual skew and kurtosis) to `output/diagnostics/`, upserted on (`spec_label`, `scale`).
+The supplement still reads the stable `output/<scale>/plots/check_joint.png` path.
+This exists so the model specifications planned next — a skew-tolerant error family, and a method-by-target interaction — can be compared against this fit rather than silently overwriting it.
+
+Corrected `CLAUDE.md`, which pointed at `report/quarto/supplement/_supplement.qmd` and described root render wrappers that do not exist; the supplement is at `report/supplement.qmd` and `_quarto.yml` renders it directly.
+
+Note: `R/sensitivity/check-autocorrelation.R` sources `R/sensitivity/model-logresp.R`, deleted in the change below, so it does not currently run. Its filter was narrowed for consistency but the script needs rebasing onto another residual source before it is usable again.
+
 ## Unreleased — Reorganise supplement; drop double-log and log-response sensitivity arms
 
 `report/quarto/supplement/_supplement.qmd`, `R/analysis-model.R`, `R/sensitivity/model-logresp.R`, `R/plot-model-flow.R`, `CLAUDE.md`

R/analysis-model.R

@@ -36,12 +36,58 @@ m.formula_joint <- wis ~
   s(Horizon, by = Model, k = 3, bs = "sz") +
   s(Model, bs = "re")
 
+# Archive a fit's diagnostics under a stable label so successive model
+# specifications can be compared rather than overwriting each other.
+# Appends to output/diagnostics/fit-summary.csv, upserting on
+# (spec_label, scale) so re-running a spec replaces its own row.
+archive_diagnostics <- function(fit, spec_label, scoring_scale, plot,
+                                dir = here("output", "diagnostics")) {
+  dir.create(dir, recursive = TRUE, showWarnings = FALSE)
+  ggsave(file.path(dir, paste0(spec_label, "_", scoring_scale, "_check.png")),
+         plot, dpi = 300)
+
+  # Deviance residuals put every family on a comparable footing; the raw
+  # response residuals of a log-link Gaussian would confound family choice
+  # with the skew of WIS itself.
+  r <- residuals(fit, type = "deviance")
+  r <- r[is.finite(r)]
+  centred <- r - mean(r)
+  # Population (biased) moment estimates: with n ~ 5e5 the small-sample
+  # correction is negligible and this keeps the values comparable to the
+  # skewness already quoted in the supplement.
+  m2 <- mean(centred^2)
+  row <- tibble::tibble(
+    spec_label = spec_label,
+    scale = scoring_scale,
+    family = fit$family$family,
+    link = fit$family$link,
+    formula = paste(deparse(formula(fit)), collapse = " "),
+    n = length(fit$y),
+    aic = AIC(fit),
+    dev_expl = summary(fit)$dev.expl,
+    resid_skew = mean(centred^3) / m2^(3 / 2),
+    resid_kurtosis = mean(centred^4) / m2^2,
+    fitted_on = as.character(Sys.Date())
+  )
+
+  path <- file.path(dir, "fit-summary.csv")
+  if (file.exists(path)) {
+    row <- read_csv(path, show_col_types = FALSE) |>
+      # coerce so a previously-written column type can't block the bind
+      mutate(across(everything(), as.character)) |>
+      filter(!(spec_label == row$spec_label & scale == row$scale)) |>
+      bind_rows(mutate(row, across(everything(), as.character)))
+  }
+  write_csv(row, path)
+  invisible(row)
+}
+
 model_wis <- function(scoring_scale = "log", family_link = "log",
- output_dir = "output") {
+ output_dir = "output", spec_label = NULL) {
   # --- Data handling ---
   m.data <- process_data(scoring_scale = scoring_scale)
   m.data <- m.data |>
-    filter(!grepl("EuroCOVIDhub-", Model)) |>
+    filter(!grepl("EuroCOVIDhub-ensemble", Model)) |>
     filter(!is.na(wis)) |> # drop unscored forecasts explicitly (bam would drop these silently)
     mutate(Epi_target = as.factor(epi_target))
 
@@ -154,4 +200,11 @@ model_wis <- function(scoring_scale = "log", family_link = "log",
   # to ~20MB per scale. PNG keeps it under 1MB with no loss of legibility.
   p <- appraise(m.fits_joint)
   ggsave(here(output_dir, "plots", "check_joint.png"), p, dpi = 300)
+
+  # Keep a labelled copy plus summary statistics, so this fit stays comparable
+  # against the specifications tried in later work. The path above is the one
+  # the supplement reads, so it deliberately stays stable.
+  if (!is.null(spec_label)) {
+    archive_diagnostics(m.fits_joint, spec_label, scoring_scale, p)
+  }
 }

R/plot-model-flow.R

@@ -53,8 +53,8 @@ create_model_flow <- function() {
   models3 <- distinct(fc_clean, target_variable, model) |>
     mutate(inc_target = TRUE)
 
-  # (4) Exclude Hub-created models
-  fc_clean <- filter(fc_clean, !grepl("EuroCOVIDhub-", model))
+  # (4) Exclude Hub-created ensemble
+  fc_clean <- filter(fc_clean, !grepl("EuroCOVIDhub-ensemble", model))
   models4 <- distinct(fc_clean, target_variable, model) |>
     mutate(inc_xhub = TRUE)
 
@@ -79,7 +79,7 @@ create_model_flow <- function() {
                            label = "Provided 1:4 week predictions",
                            show_exc = TRUE) |>
                  fc_filter(inc_xhub,
-                           label = "Not created by Hub",
+                           label = "Not the Hub ensemble",
                            show_exc = TRUE) |>
                  fc_draw()
   )

R/sensitivity/check-autocorrelation.R

@@ -46,7 +46,7 @@ check_autocorrelation <- function() {
   # Reconstruct m.data with the SAME filters/mutate as model-logresp.R:29-37 so
   # rows align 1:1 with fit$y (verified by stopifnot below).
   m.data <- process_data(scoring_scale = "log") |>
-    filter(!grepl("EuroCOVIDhub-", Model)) |>
+    filter(!grepl("EuroCOVIDhub-ensemble", Model)) |>
     filter(!is.na(wis)) |>
     mutate(
       Epi_target = as.factor(epi_target),

R/sensitivity/check-link-robustness.R

@@ -19,7 +19,7 @@ source(here("R", "process-data.R"))
 check_link_robustness <- function() {
   # --- Data: mirror model_wis(scoring_scale = "log") ---
   m.data <- process_data(scoring_scale = "log") |>
-    filter(!grepl("EuroCOVIDhub-", Model)) |>
+    filter(!grepl("EuroCOVIDhub-ensemble", Model)) |>
     mutate(Incidence = log(Incidence + 1))
 
   m.formula_joint <- wis ~

README.md

@@ -1,4 +1,5 @@
 
+
 [![Zenodo](https://img.shields.io/badge/Code%20DOI-10.5281/zenodo.14903161-blue)](https://doi.org/10.5281/zenodo.14903161)
 [![medRxiv](https://img.shields.io/badge/medRxiv-10.1101/2025.04.10.25325611-blue)](https://doi.org/10.1101/2025.04.10.25325611)
 
@@ -32,7 +33,7 @@ method, from the difficulty of the target, in forecast performance.
 We evaluated forecasts of weekly COVID-19 cases and deaths over two
 years across 32 European countries, scoring them against observed data
 with the Weighted Interval Score (WIS). We expected a model’s structure
-to shape how well it predicted, so we classified 47 models by structure
+to shape how well it predicted, so we classified 48 models by structure
 (agent-based, mechanistic, semi-mechanistic, statistical, or human
 judgement) and estimated how much structure alone affected performance.
 A generalised additive mixed model let us adjust for everything that
@@ -66,8 +67,9 @@ Read the work as it stands:
   section files in [report/quarto/](./report/quarto/).
 - The supplement is in
   [report/quarto/supplement/\_supplement.qmd](./report/quarto/supplement/_supplement.qmd).
-- The pre-print is on
-  [medRxiv](https://doi.org/10.1101/2025.04.10.25325611).
+- The former pre-print is on
+  [medRxiv](https://doi.org/10.1101/2025.04.10.25325611); note, this is
+  out of date.
 
 #### Reproduce
 
@@ -135,10 +137,6 @@ A guide to the codebase:
 
 ### Project docs
 
-- Current [working
-  draft](https://docs.google.com/document/d/1OOVUHR_BGWcviSNxvaHvbXD16Bb3Y_zhw--7gAGBqMk/edit#)
-- Submitted
-  [manuscript](https://docs.google.com/document/d/1B_HviobjSIak4c1FKoSOqgFwRF_DxMrJM6zsdbdth4E/edit?tab=t.0)
 - [Authorship](https://docs.google.com/spreadsheets/d/18mt6c47MCzLdMAKth6Bv2PC7b10KYgWieExHHWLnP3Q/edit?gid=0#gid=0)
 - [Slide
   deck](https://docs.google.com/presentation/d/1BSdTEuZ_zKdU8tBFuRMmP7GwHht1D0oZSkaFWovz9ao/edit?slide=id.p#slide=id.p)

output/diagnostics/fit-summary.csv

@@ -0,0 +1,3 @@
+spec_label,scale,family,link,formula,n,aic,dev_expl,resid_skew,resid_kurtosis,fitted_on
+baseline-included-gaussian-log,log,gaussian,log,"wis ~ Epi_target + s(Method, bs = ""re"") + s(CountryTargets, bs = ""re"") +      s(Incidence) + s(Trend, bs = ""re"") + s(Location, bs = ""re"") +      s(VariantPhase, bs = ""re"") + s(Horizon, by = Model, k = 3,      bs = ""sz"") + s(Model, bs = ""re"")",207713,262977.342419348,0.285523068894277,5.84369437663132,77.4753183536698,2026-07-28
+baseline-included-gaussian-log,natural,Gamma,log,"wis ~ Epi_target + s(Method, bs = ""re"") + s(CountryTargets, bs = ""re"") +      s(Incidence) + s(Trend, bs = ""re"") + s(Location, bs = ""re"") +      s(VariantPhase, bs = ""re"") + s(Horizon, by = Model, k = 3,      bs = ""sz"") + s(Model, bs = ""re"")",207713,2447775.88168385,0.642710457169962,4.58628582506588,122.817170132029,2026-07-28

report/quarto/_abstract.qmd

@@ -5,7 +5,7 @@ However, when many forecasters each predict many different targets, it becomes d
 Here we use a regression model to separate the effect of the forecasting method, from the difficulty of the target, in forecast performance.
 
 We evaluated forecasts of weekly COVID-19 cases and deaths over two years across 32 European countries, scoring them against observed data with the Weighted Interval Score (WIS).
-We expected a model's structure to shape how well it predicted, so we classified 47 models by structure (agent-based, mechanistic, semi-mechanistic, statistical, or human judgement) and estimated how much structure alone affected performance.
+We expected a model's structure to shape how well it predicted, so we classified 48 models by structure (agent-based, mechanistic, semi-mechanistic, statistical, or human judgement) and estimated how much structure alone affected performance.
 A generalised additive mixed model let us adjust for everything that makes a target easier or harder to predict: the outcome being forecast, its level and trend, the dominant variant, the country, the forecast horizon, and differences between individual models.
 
 Once we accounted for the difficulty of the target, no single type of model performed best.