diff --git a/.gitignore b/.gitignore index 25fc1a1..94bf63c 100644 --- a/.gitignore +++ b/.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/ diff --git a/CLAUDE.md b/CLAUDE.md index b4c0ae1..fd372b1 100644 --- a/CLAUDE.md +++ b/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//. 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. diff --git a/NEWS.md b/NEWS.md index cc9915a..5185a4c 100644 --- a/NEWS.md +++ b/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//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` diff --git a/R/analysis-model.R b/R/analysis-model.R index 1c6eb78..d1c2bf3 100644 --- a/R/analysis-model.R +++ b/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) + } } diff --git a/R/plot-model-flow.R b/R/plot-model-flow.R index 4cb6864..825b2e8 100644 --- a/R/plot-model-flow.R +++ b/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() ) diff --git a/R/sensitivity/check-autocorrelation.R b/R/sensitivity/check-autocorrelation.R index c23b45f..2750349 100644 --- a/R/sensitivity/check-autocorrelation.R +++ b/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), diff --git a/R/sensitivity/check-link-robustness.R b/R/sensitivity/check-link-robustness.R index d3ab81f..e612d3b 100644 --- a/R/sensitivity/check-link-robustness.R +++ b/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 ~ diff --git a/README.md b/README.md index 25d28fa..94ecdaa 100644 --- a/README.md +++ b/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) diff --git a/output/diagnostics/baseline-included-gaussian-log_log_check.png b/output/diagnostics/baseline-included-gaussian-log_log_check.png new file mode 100644 index 0000000..5e6b9c5 Binary files /dev/null and b/output/diagnostics/baseline-included-gaussian-log_log_check.png differ diff --git a/output/diagnostics/baseline-included-gaussian-log_natural_check.png b/output/diagnostics/baseline-included-gaussian-log_natural_check.png new file mode 100644 index 0000000..6559c9b Binary files /dev/null and b/output/diagnostics/baseline-included-gaussian-log_natural_check.png differ diff --git a/output/diagnostics/fit-summary.csv b/output/diagnostics/fit-summary.csv new file mode 100644 index 0000000..d480b5c --- /dev/null +++ b/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 diff --git a/output/flowchart.png b/output/flowchart.png index abd8749..1e2c0b1 100644 Binary files a/output/flowchart.png and b/output/flowchart.png differ diff --git a/output/log/fit_obs.rds b/output/log/fit_obs.rds index b4d4a80..04aa557 100644 Binary files a/output/log/fit_obs.rds and b/output/log/fit_obs.rds differ diff --git a/output/log/plots/check_joint.png b/output/log/plots/check_joint.png index 4fee9db..5e6b9c5 100644 Binary files a/output/log/plots/check_joint.png and b/output/log/plots/check_joint.png differ diff --git a/output/log/results.rds b/output/log/results.rds index 2fefb37..abff22c 100644 Binary files a/output/log/results.rds and b/output/log/results.rds differ diff --git a/output/natural/fit_obs.rds b/output/natural/fit_obs.rds index 1d63fe9..75f6bbf 100644 Binary files a/output/natural/fit_obs.rds and b/output/natural/fit_obs.rds differ diff --git a/output/natural/plots/check_joint.png b/output/natural/plots/check_joint.png index 422ba91..6559c9b 100644 Binary files a/output/natural/plots/check_joint.png and b/output/natural/plots/check_joint.png differ diff --git a/output/natural/results.rds b/output/natural/results.rds index 20d6b78..5d76bf2 100644 Binary files a/output/natural/results.rds and b/output/natural/results.rds differ diff --git a/report/quarto/_abstract.qmd b/report/quarto/_abstract.qmd index d537cba..69ffe8e 100644 --- a/report/quarto/_abstract.qmd +++ b/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. diff --git a/report/quarto/_discussion.qmd b/report/quarto/_discussion.qmd index 0d4f043..f8f0dad 100644 --- a/report/quarto/_discussion.qmd +++ b/report/quarto/_discussion.qmd @@ -5,11 +5,13 @@ Adjusted estimates for each model structure were no different from the overall a We suggest that when evaluating across multiple forecast targets, it is necessary to account for the factors contributing to variation between targets, before considering factors within the forecasting process. In this work we had limited ability to demonstrate the model-based approach used here, with a small and biased sample reducing any power to detect true differences in performance between model structures. -While we used a large dataset of predictions over time, we had a much smaller effective sample size of 47 independent models. +While we used a large dataset of predictions over time, we had a much smaller effective sample size of 48 independent models. We relied on an opportunistic cohort of models from voluntary contributions to the European Hub. For example, we only observed 3 agent-based models, where all three represented here targeted a single country. This meant our results may be biased by unobserved characteristics that differentially affected participation, such as the resource intensity of model development, or the capacity to incorporate domain expertise. These factors may drive both the choice of model structure and be more useful as targets for identifying potential for modifying forecaster performance. +We included the Hub baseline model, which is a deliberately naive statistical model submitting for every country in almost every week. +It therefore contributes a large share of the statistical group's forecasts, and anchors that group's estimate towards the performance achievable without any epidemiological structure. In addition, we may have misclassified the exposure, with very little metadata to classify model structures. Independent raters disagreed on the classification of most semi-mechanistic models, and we did not track changes to each model's individual methods over time. A more systematic sample of models would support a deeper investigation of factors driving model performance. diff --git a/report/quarto/_methods.qmd b/report/quarto/_methods.qmd index 5bf41dd..7e46cc5 100644 --- a/report/quarto/_methods.qmd +++ b/report/quarto/_methods.qmd @@ -8,7 +8,9 @@ We used forecast data collected via the public European COVID-19 Forecast Hub, d The platform solicited real-time forecasts for between one and four weeks ahead for 32 European countries. Any forecaster was eligible to participate in the Hub, and there were no selection criteria. -We excluded forecasts that did not report the full set of 23 quantiles, to ensure fair comparison among probabilistic results, and excluded baseline and ensemble models created by the Hub team. +We excluded forecasts that did not report the full set of 23 quantiles, to ensure fair comparison among probabilistic results. +We excluded the ensemble model created by the Hub team, because it is a function of the contributed forecasts and would therefore double-count them. +We retained the Hub baseline model, which is an independently specified statistical model and so contributes to the sample on the same terms as any other participant. Full study eligibility and a STROBE flow diagram are detailed in the Supplement. We evaluated forecasts against observed data collated by Johns Hopkins University (JHU). diff --git a/report/quarto/_results.qmd b/report/quarto/_results.qmd index 5c532c9..68478ee 100644 --- a/report/quarto/_results.qmd +++ b/report/quarto/_results.qmd @@ -27,11 +27,11 @@ scores <- process_data(scoring_scale = "log") ensemble <- scores |> filter(grepl("EuroCOVIDhub-ensemble", Model)) scores <- scores |> - filter(!grepl("EuroCOVIDhub-", Model)) + filter(!grepl("EuroCOVIDhub-ensemble", Model)) # natural-scale scores for descriptive error-vs-observation figure scores_natural <- process_data(scoring_scale = "natural") |> - filter(!grepl("EuroCOVIDhub-", Model)) + filter(!grepl("EuroCOVIDhub-ensemble", Model)) ``` @@ -72,6 +72,20 @@ structure_count <- structures |> summarise(agreement = sum(agreement), n = n(), p = agreement / n) +# Named lookups so the text below tracks the data rather than hardcoding counts +n_by_method <- setNames(structure_count$n, as.character(structure_count$Method)) +n_disagree <- sum(!structures$agreement) +p_disagree <- round(n_disagree / nrow(structures) * 100) +p_by_method <- scores |> + count(Method) |> + mutate(p = round(n / sum(n) * 100)) |> + (\(x) setNames(x$p, as.character(x$Method)))() +# Semi-mechanistic models that drew at least one statistical vote +semi_mech_votes <- classify_models(return_majority = FALSE) |> + filter(model %in% structures$Model[structures$Method == "Semi-mechanistic"]) +n_semi_mech_split <- n_distinct( + semi_mech_votes$model[semi_mech_votes$classification == "Statistical"] +) # Country targets targets <- scores |> @@ -113,12 +127,13 @@ Over the study period, forecasters could participate in any of 26,624 available Median forecaster participation was `r model_forecasts[["Median"]]`% (IQR: `r model_forecasts[["1st Qu."]]`-`r model_forecasts[["3rd Qu."]]`%, ranging from `r model_forecasts[["Min."]]` to `r model_forecasts[["Max."]]`) across available targets, reflecting substantial variation in participation across countries and time. -We classified 17 models as using mechanistic structure, 12 as statistical, 12 as semi-mechanistic, 3 as agent-based and 3 that used human judgement forecasting (Table 1; Supplementary Table S2). -Raters disagreed on 17 (36%) model classifications. +We classified `r n_by_method[["Mechanistic"]]` models as using mechanistic structure, `r n_by_method[["Statistical"]]` as statistical, `r n_by_method[["Semi-mechanistic"]]` as semi-mechanistic, `r n_by_method[["Agent-based"]]` as agent-based and `r n_by_method[["Judgement"]]` that used human judgement forecasting (Table 1; Supplementary Table S2). +Raters disagreed on `r n_disagree` (`r p_disagree`%) model classifications. The majority of 2/3 was used as the final classification, with additional manual review which in all cases retained the majority decision. -Raters most often disagreed when classifying semi-mechanistic models, with 8 out of 12 receiving one or more votes for a statistical model structure. -Under the final classification, approximately one-third of all forecast predictions were from mechanistic, semi-mechanistic, and statistical models each. -Agent-based and judgement models provided fewer forecasts, representing <2% of forecasts. +Raters most often disagreed when classifying semi-mechanistic models, with `r n_semi_mech_split` out of `r n_by_method[["Semi-mechanistic"]]` receiving one or more votes for a statistical model structure. +Under the final classification, statistical models contributed the largest share of forecast predictions (`r p_by_method[["Statistical"]]`%), partly because the Hub baseline model submitted for every country in almost every week. +Mechanistic and semi-mechanistic models contributed `r p_by_method[["Mechanistic"]]`% and `r p_by_method[["Semi-mechanistic"]]`% respectively. +Agent-based and judgement models provided fewer forecasts, each representing under 2% of forecasts. We classified `r targets[["Single-country"]]` models as forecasting for a single country and `r targets[["Multi-country"]]` as forecasting for multiple countries; the countries targeted are detailed in the Supplement. Among the `r targets[["Multi-country"]]` multi-country models, only `r nrow(multi_targets |> filter(consistent))` consistently forecast for the same number of countries throughout the study period. diff --git a/report/supplement.qmd b/report/supplement.qmd index abc40b0..3cf0fc2 100644 --- a/report/supplement.qmd +++ b/report/supplement.qmd @@ -71,7 +71,8 @@ Forecasts were validated against minimal formatting requirements for quantile in For this study, we collected all forecasts from between 8 March 2021 to 10 March 2023. We excluded forecasts of hospitalisations, which experienced multiple changes in source data during the study period. We excluded forecasts that did not report the full set of 23 quantiles, in order to ensure fair comparison among probabilistic model results. -We also excluded baseline and ensemble models created by the Hub team. +We also excluded the ensemble model created by the Hub team, because it is constructed from the contributed forecasts and so would double-count them. +We retained the Hub baseline model: it is an independently specified statistical model, and including it gives a naive reference point against which participant forecasts can be read. ```{r} #| label: supplement-model-flow @@ -425,7 +426,7 @@ scores <- process_data(scoring_scale = "natural") ensemble <- scores |> filter(grepl("EuroCOVIDhub-ensemble", Model)) scores <- scores |> - filter(!grepl("EuroCOVIDhub-", Model)) + filter(!grepl("EuroCOVIDhub-ensemble", Model)) print_table1(scores) ```