diff --git a/CLAUDE.md b/CLAUDE.md index fd372b1..135137c 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -123,9 +123,9 @@ source(here("R", "process-data.R")) # 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"), - spec_label = "baseline-included-gaussian-log") + spec_label = "tweedie-log") model_wis(scoring_scale = "natural", output_dir = here("output", "natural"), - spec_label = "baseline-included-gaussian-log") + spec_label = "tweedie-log") # 4. Render the manuscript alone (results section only; supplement is a separate page) # quarto::quarto_render("report/manuscript.qmd") diff --git a/NEWS.md b/NEWS.md index 5185a4c..19f3d9a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,35 @@ Notable changes to the analysis, manuscript, and repository. Newest first. +## Unreleased — Model WIS with a Tweedie family (#159) + +`R/analysis-model.R`, `R/sensitivity/check-family.R`, `report/quarto/_methods.qmd`, `report/quarto/_results.qmd`, `report/supplement.qmd` + +The primary model used `gaussian(link = "log")`, which left deviance residuals with skew 5.8 and kurtosis 77. +Modelling `log(WIS)` directly would fix the residuals but loses propriety of the score, so the fix had to come from the error family instead. + +Compared Gaussian, Gamma and Tweedie families on the joint specification, holding formula and data fixed. +Both scales now use `tw(link = "log")`, replacing `gaussian(log)` on the log scale and `Gamma(log)` on the natural scale. + +On the log scale this is a large improvement: residual skew falls from 5.84 to 0.58, kurtosis from 77.5 to 9.2, and deviance explained rises from 0.286 to 0.380. +Gamma fits the same data almost identically (skew 0.52, deviance explained 0.378) but does not converge on either scale, which is the reason for preferring Tweedie. +The Tweedie power parameter is estimated at 1.99, the upper limit `mgcv` permits, so the fitted family is a Gamma in all but numerical behaviour. + +This also resolves the natural-scale non-convergence recorded in the previous entry: that was a Gamma problem, not a scale problem. +On the natural scale the change fixes convergence but not the fit — residual skew is unchanged at 4.59, because natural-scale WIS is skewed beyond what any Tweedie can absorb. +Nothing in the rendered manuscript or supplement reads `output/natural/`, so this affects no reported result. + +Several adjusted estimates moved materially under the new family, most notably the deaths-versus-cases contrast (ratio 0.17 to 0.38). +All substantive conclusions hold: no model structure differs from the grand mean, stable trends remain the most predictable, increasing trends the least, and Omicron BA.1 the hardest variant phase. +Delta's interval now excludes 1, where previously it did not. + +Investigated whether the `1e-7` constant that `process-data.R` adds to every score was driving the skew, since 553 forecasts (0.27%) score exactly zero and the constant parks them 11 log-units below the next smallest value. +It was not: refitting with the constant removed and the exact zeros retained changes residual skew by 0.01. +`process-data.R` is therefore unchanged, and the result is recorded in the supplement as a negative finding. + +Fixed `archive_diagnostics()`, which reassigned its accumulator and so returned the whole `fit-summary.csv` rather than the row just written. +Fixed a non-standard-evaluation trap: `tw()` deparses its `link` argument, so passing a variable sent the literal string `"family_link"`. + ## 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` diff --git a/R/analysis-model.R b/R/analysis-model.R index d1c2bf3..de34ace 100644 --- a/R/analysis-model.R +++ b/R/analysis-model.R @@ -9,7 +9,8 @@ # Horizon: forecast horizon (smooth, by model) # Model: individual model (random effect) # -# Response: WIS (log-transformed, Gaussian family with log link) +# Response: WIS, modelled with a Tweedie family and log link on both scales. +# See R/sensitivity/check-family.R library(here) library(dplyr) @@ -71,19 +72,24 @@ archive_diagnostics <- function(fit, spec_label, scoring_scale, plot, ) path <- file.path(dir, "fit-summary.csv") + summary_table <- mutate(row, across(everything(), as.character)) if (file.exists(path)) { - row <- read_csv(path, show_col_types = FALSE) |> + summary_table <- 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))) + bind_rows(summary_table) } - write_csv(row, path) + write_csv(summary_table, path) invisible(row) } -model_wis <- function(scoring_scale = "log", family_link = "log", - output_dir = "output", spec_label = NULL) { +model_wis <- function( + scoring_scale = "log", + family_link = "log", + output_dir = "output", + spec_label = NULL +) { # --- Data handling --- m.data <- process_data(scoring_scale = scoring_scale) m.data <- m.data |> @@ -91,17 +97,17 @@ model_wis <- function(scoring_scale = "log", family_link = "log", filter(!is.na(wis)) |> # drop unscored forecasts explicitly (bam would drop these silently) mutate(Epi_target = as.factor(epi_target)) - # Settings for log or natural scale + # Settings for log or natural scale. Both scales use the same family if (scoring_scale == "log") { # log-transform incidence to match scoring on log scale m.data <- m.data |> mutate(Incidence = log(Incidence + 1)) - m.family <- gaussian(link = family_link) - } else if (scoring_scale == "natural") { - m.family <- Gamma(link = family_link) - } else { + } else if (scoring_scale != "natural") { stop("scoring_scale must be either 'log' or 'natural'") } + # tw() deparses its `link` argument, so passing the variable directly would + # send the literal string "family_link". do.call forces the value through. + m.family <- do.call(tw, list(link = family_link)) # --- Model formula --- # Univariate for each @@ -151,7 +157,8 @@ model_wis <- function(scoring_scale = "log", family_link = "log", transmute( group_var = "Epi_target", group = "Deaths", - value, se, + value, + se, lower_2.5 = .data[[ci_cols[grepl("^lower", ci_cols)]]], upper_97.5 = .data[[ci_cols[grepl("^upper", ci_cols)]]], model = model_label @@ -159,7 +166,9 @@ model_wis <- function(scoring_scale = "log", family_link = "log", } # Univariate random effects (exclude smooth-only and the fixed target fit) - random_effects_uni <- m.fits_uni[!grepl("horizon|incidence|epi_target", names(m.fits_uni))] |> + random_effects_uni <- m.fits_uni[ + !grepl("horizon|incidence|epi_target", names(m.fits_uni)) + ] |> map(extract_ranef) |> list_rbind() |> mutate(model = "Unadjusted") |> @@ -196,14 +205,11 @@ model_wis <- function(scoring_scale = "log", family_link = "log", ) saveRDS(fit_obs, here(output_dir, "fit_obs.rds")) - # Raster, not vector: appraise() plots ~150k residuals, and a PDF of that runs - # to ~20MB per scale. PNG keeps it under 1MB with no loss of legibility. + # appraise() plots 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/sensitivity/check-family.R b/R/sensitivity/check-family.R new file mode 100644 index 0000000..f0c5afa --- /dev/null +++ b/R/sensitivity/check-family.R @@ -0,0 +1,113 @@ +# Sensitivity: choice of error family for the log-scale (primary) GAMM. +# +# +# Two considerations: +# +# 1. WIS on the log scale is continuous, positive, and strongly right-skewed. +# Gamma and Tweedie are the natural candidates. Symmetric heavy-tailed +# families (e.g. scat()) are deliberately excluded: they could only +# downweight the tail, not represent the skew, and they put support on +# negative values, which is wrong for a strictly positive score. +# +# 2. 553 forecasts (0.27%) have WIS exactly 0, as perfect predictions of +# zero-incidence targets, almost all deaths in small countries (Iceland, +# Liechtenstein, Malta). process-data.R adds 1e-7 to every score so these +# are representable on a log link, which parks them at log(1e-7) = -16.1, +# roughly 11 log-units below the next smallest score. +# +# Tweedie with 1 < p < 2 has a genuine point mass at zero, so it can model +# those forecasts as what they are instead of displacing them. The +# "tweedie-nooffset" arm therefore removes the 1e-7 and keeps the exact +# zeros. +# Using gaussian(link = "log") leaves strongly skewed +# deviance residuals (skew ~5.8, kurtosis ~77), which is a poor description of +# the outcome. +# +# Run: source(here::here("R", "sensitivity", "check-family.R")); check_family() + +library(here) +library(dplyr) +library(readr) +library(purrr) +library(mgcv) +library(ggplot2) +library(gratia) +source(here("R", "analysis-model.R")) # m.formula_joint, archive_diagnostics() + +# Candidate families. `offset` records whether the 1e-7 added in process-data.R +# is retained; the no-offset arm is only meaningful for a family that admits +# exact zeros. +.family_candidates <- list( + list(label = "gaussian-log", family = quote(gaussian(link = "log")), offset = TRUE), + list(label = "gamma-log", family = quote(Gamma(link = "log")), offset = TRUE), + list(label = "tweedie-log", family = quote(tw(link = "log")), offset = TRUE), + list(label = "tweedie-nooffset", family = quote(tw(link = "log")), offset = FALSE) +) + +check_family <- function(candidates = .family_candidates, + spec_prefix = "family") { + m.data <- process_data(scoring_scale = "log") |> + filter(!grepl("EuroCOVIDhub-ensemble", Model)) |> + filter(!is.na(wis)) |> + mutate( + Epi_target = as.factor(epi_target), + Incidence = log(Incidence + 1) + ) + + results <- map(candidates, \(cand) { + message("-------- fitting family: ", cand$label) + dat <- m.data + if (!cand$offset) { + # Undo the constant added in process-data.R, restoring the exact zeros. + dat <- mutate(dat, wis = pmax(wis - 1e-7, 0)) + } + + # bam() signals non-convergence through a warning + warnings_seen <- character() + fit <- withCallingHandlers( + bam( + formula = m.formula_joint, + data = dat, + family = eval(cand$family), + method = "fREML", + discrete = TRUE + ), + warning = function(w) { + warnings_seen <<- c(warnings_seen, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + + p <- appraise(fit) + row <- archive_diagnostics( + fit, + spec_label = paste(spec_prefix, cand$label, sep = "-"), + scoring_scale = "log", + plot = p + ) + + converged <- !any(grepl("did not converge", warnings_seen)) + message( + " skew ", + signif(as.numeric(row$resid_skew), 3), + " kurtosis ", + signif(as.numeric(row$resid_kurtosis), 3), + " converged: ", + converged + ) + if (length(warnings_seen)) { + message(" warnings: ", paste(unique(warnings_seen), collapse = "; ")) + } + + tibble::tibble( + label = cand$label, + offset = cand$offset, + converged = converged, + warnings = paste(unique(warnings_seen), collapse = "; ") + ) + }) + + # AIC is comparable only within an offset arm: the no-offset fit has a + # different response vector, so its likelihood is on a different scale. + bind_rows(results) +} diff --git a/R/sensitivity/model-building.qmd b/attic/model-building.qmd similarity index 100% rename from R/sensitivity/model-building.qmd rename to attic/model-building.qmd diff --git a/output/diagnostics/baseline-included-gaussian-log_natural_check.png b/output/diagnostics/baseline-included-gaussian-log_natural_check.png deleted file mode 100644 index 6559c9b..0000000 Binary files a/output/diagnostics/baseline-included-gaussian-log_natural_check.png and /dev/null differ diff --git a/output/diagnostics/family-gamma-log_log_check.png b/output/diagnostics/family-gamma-log_log_check.png new file mode 100644 index 0000000..b5f8057 Binary files /dev/null and b/output/diagnostics/family-gamma-log_log_check.png differ diff --git a/output/diagnostics/baseline-included-gaussian-log_log_check.png b/output/diagnostics/family-gaussian-log_log_check.png similarity index 100% rename from output/diagnostics/baseline-included-gaussian-log_log_check.png rename to output/diagnostics/family-gaussian-log_log_check.png diff --git a/output/diagnostics/family-tweedie-log_log_check.png b/output/diagnostics/family-tweedie-log_log_check.png new file mode 100644 index 0000000..6c23a41 Binary files /dev/null and b/output/diagnostics/family-tweedie-log_log_check.png differ diff --git a/output/diagnostics/family-tweedie-nooffset_log_check.png b/output/diagnostics/family-tweedie-nooffset_log_check.png new file mode 100644 index 0000000..d0bbdc8 Binary files /dev/null and b/output/diagnostics/family-tweedie-nooffset_log_check.png differ diff --git a/output/diagnostics/fit-summary.csv b/output/diagnostics/fit-summary.csv index d480b5c..fa3ddb0 100644 --- a/output/diagnostics/fit-summary.csv +++ b/output/diagnostics/fit-summary.csv @@ -1,3 +1,7 @@ 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 +family-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 +family-tweedie-log,log,Tweedie(p=1.99),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,-103072.05282003,0.380078174238279,0.579174455925947,9.22170906116645,2026-07-28 +family-tweedie-nooffset,log,Tweedie(p=1.931),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,-130379.800593153,0.382245908293812,0.56727987358893,9.90754617063616,2026-07-28 +family-gamma-log,log,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,-104260.314306789,0.378042615339602,0.519239757977383,9.39292725143449,2026-07-28 +tweedie-log,log,Tweedie(p=1.99),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,-103072.004819446,0.38007819340237,0.57917415910616,9.22170793637798,2026-07-28 +tweedie-log,natural,Tweedie(p=1.99),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,1544350.87312651,0.645351609661281,4.59696577231419,123.232938667622,2026-07-28 diff --git a/output/diagnostics/tweedie-log_log_check.png b/output/diagnostics/tweedie-log_log_check.png new file mode 100644 index 0000000..8b1931d Binary files /dev/null and b/output/diagnostics/tweedie-log_log_check.png differ diff --git a/output/diagnostics/tweedie-log_natural_check.png b/output/diagnostics/tweedie-log_natural_check.png new file mode 100644 index 0000000..3709db0 Binary files /dev/null and b/output/diagnostics/tweedie-log_natural_check.png differ diff --git a/output/log/fit_obs.rds b/output/log/fit_obs.rds index 04aa557..abb4151 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 5e6b9c5..8b1931d 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 abff22c..89e3a0d 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 75f6bbf..7f2033c 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 6559c9b..3709db0 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 5d76bf2..c8aaac3 100644 Binary files a/output/natural/results.rds and b/output/natural/results.rds differ diff --git a/report/quarto/_methods.qmd b/report/quarto/_methods.qmd index 7e46cc5..790edb4 100644 --- a/report/quarto/_methods.qmd +++ b/report/quarto/_methods.qmd @@ -58,12 +58,16 @@ We did not attempt to adjust for unobserved characteristics of the forecasting t **Analysis** We used a single hierarchical model structure fit across both epidemiological targets, with target included as a fixed factor. -We treated the WIS of each forecast $i$ as Gaussian with a log link, so that the linear predictor $\eta_i$ acts multiplicatively on expected score: +The WIS is continuous, non-negative, and strongly right-skewed, so we treated the WIS of each forecast $i$ as Tweedie-distributed with a log link, letting the linear predictor $\eta_i$ act multiplicatively on the expected score: $$ -\mathrm{WIS}_i \sim \mathcal{N}(\mu_i,\, \sigma^2), \qquad \log \mu_i = \eta_i. +\mathrm{WIS}_i \sim \mathrm{Tw}_p(\mu_i,\, \phi), \qquad \log \mu_i = \eta_i, $$ +where $\phi$ is the dispersion and the power parameter $p \in (1, 2)$ is estimated alongside the smoothing parameters, giving a variance function $\mathrm{Var}(\mathrm{WIS}_i) = \phi\,\mu_i^{p}$. +We selected this family by comparing it against a Gaussian and a Gamma family on the same specification, reported in the Supplement. +The estimated $p$ reached 1.99, the upper limit of the permitted range, so the fitted family is effectively a Gamma; we report the Tweedie fit because the Gamma parameterisation did not converge on these data. + The linear predictor decomposed into a fixed intercept, a fixed effect for the epidemiological target, a sum of random-effect contributions from categorical covariates, and two smooth terms for continuous covariates: $$ diff --git a/report/quarto/_results.qmd b/report/quarto/_results.qmd index 68478ee..ab0edf6 100644 --- a/report/quarto/_results.qmd +++ b/report/quarto/_results.qmd @@ -192,13 +192,13 @@ effects_comp <- results$effects |> We fitted a generalised additive mixed model to give adjusted estimates of the partial effect of model structure, while controlling for varying forecaster target selection and epidemic dynamics between targets. In this structure, partial effects are deviations from the grand mean WIS under a sum-to-zero constraint, so a negative value indicates better-than-average performance. -The WIS was highly right-skewed, and we used a log link to account for this (diagnostics in Supplementary Figures S2-S3). -The log link models the score multiplicatively but leaves some skew in the residuals. -A sensitivity analysis modelling the log-transformed score directly substantially improves the residual distribution while preserving the direction of all effects and leaving the model-structure conclusions unchanged (Supplement). +The WIS was highly right-skewed, and we modelled it with a Tweedie family and a log link to account for this (diagnostics in Supplementary Figures S2-S3). +The log link models the score multiplicatively, and the Tweedie family accommodates the skew directly rather than leaving it in the residuals. +We compared this against a Gaussian and a Gamma family: the Gaussian left strongly skewed residuals, and the Gamma fitted comparably to the Tweedie but did not converge (Supplement). Exponentiating a partial effect gives a multiplicative ratio relative to the grand-mean WIS, where 1 indicates average performance (e.g. a partial effect of −0.1 corresponds to a ratio of `r round(exp(-0.1), 2)`, a WIS approximately `r round((1 - exp(-0.1)) * 100, digits=0)`% lower than average). We report the exponentiated ratio (with 95% confidence intervals) in the main text (Table 2), with the raw partial effects on the log scale in the Supplement. -After adjustment forfeatures of the forecast target , no single structural approach dominated. +After adjustment for features of the forecast target, no single structural approach dominated. Adjusted point estimates clustered around the grand mean, and confidence intervals overlapped throughout (Table 2). The largest shifts were among agent-based and human judgement models, which appeared better than average in unadjusted estimates but showed no difference from other model structures after adjustment (adjusted ratios within 1% of the grand mean). We noted that adjustment consistently shrank the standard error compared to univariate estimates, narrowing the intervals around these overlapping estimates (@fig-plot-coeffs). diff --git a/report/supplement.qmd b/report/supplement.qmd index 3cf0fc2..b7e2df5 100644 --- a/report/supplement.qmd +++ b/report/supplement.qmd @@ -315,7 +315,51 @@ m.data |> geom_density() ``` -We used a \[\] distribution [skew, kurtosis]. +We therefore compared three error families on this specification, holding the formula and data fixed and varying only the family. +The Gaussian family leaves the skew in the residuals rather than modelling it. +The Gamma and Tweedie families both accommodate a right-skewed positive outcome, and fit this data almost identically. +This is not a coincidence: `mgcv` restricts the Tweedie power parameter to $1 < p < 2$, and the estimate reaches 1.99, the top of that range. +The data are therefore asking for a variance function at least as steep as the Gamma ($p = 2$), and the fitted Tweedie is a Gamma in all but name. +The two differ only in numerical behaviour, and the Gamma fit did not converge on either scale, so we report the Tweedie fit. + +The choice of family matters a great deal on the log scale and very little on the natural scale. +On the log scale, moving from Gaussian to Tweedie reduces residual skew from 5.8 to 0.6 and kurtosis from 77 to 9. +On the natural scale, WIS is skewed beyond what any Tweedie can absorb, and the residual skew is essentially unchanged (4.59 against 4.59 under a Gamma); there the change resolves the convergence failure rather than improving the distributional fit. + +```{r} +#| label: supplement-family-comparison + +read_csv(here("output", "diagnostics", "fit-summary.csv"), + show_col_types = FALSE) |> + filter(grepl("^family-", spec_label), scale == "log") |> + # the no-offset arm has a different response vector, so its likelihood is not + # comparable with the others; it is reported separately below + filter(spec_label != "family-tweedie-nooffset") |> + transmute( + Family = family, + `Deviance explained` = round(as.numeric(dev_expl), 3), + AIC = round(as.numeric(aic)), + `Residual skew` = round(as.numeric(resid_skew), 2), + `Residual kurtosis` = round(as.numeric(resid_kurtosis), 1) + ) |> + kable( + caption = paste0( + "Comparison of error families for the joint model on the log scale. ", + "Skew and kurtosis are of the deviance residuals, which place each family ", + "on a comparable footing. Lower absolute skew and kurtosis indicate a ", + "distributional assumption better matched to the outcome. ", + "The Gamma fit did not converge, and is shown for comparison only; ", + "the Tweedie fit is the one reported throughout." + ), + align = c("l", "r", "r", "r", "r") + ) |> + kable_styling(full_width = FALSE) +``` + +A small share of forecasts (553, or 0.27%) score exactly zero: perfect predictions of zero-incidence targets, almost all deaths in the smallest countries. +`process-data.R` adds a constant of 1e-7 to every score so these remain representable on a log link. +Because that constant is arbitrary, we refitted the Tweedie model with the constant removed and the exact zeros retained, which a Tweedie family with power parameter between 1 and 2 admits directly. +This made almost no difference (residual skew 0.57 against 0.58, deviance explained 0.382 against 0.380), so we retained the constant. We used the `mgcv` package v1.9-4 using R 4.5, with the formula: @@ -454,8 +498,9 @@ We tried several alternative model specifications. #### Fitting -We tried a Gaussian distribution with a log link. -Because LWIS is strongly right-skewed, this resulted in high residual skewness (5.5). +We first used a Gaussian distribution with a log link. +Because LWIS is strongly right-skewed, this left high residual skewness (5.8) and kurtosis (77). We avoid using additional data transformations (e.g. an additional log transform) on the outcome (LWIS) as this would violate propriety of the score. +We instead addressed the skew through the error family, comparing Gaussian, Gamma and Tweedie families as reported under [Model fitting](#model-fitting) above, and adopting the Tweedie family for the reported analysis. The fitted relationships and partial effects remained stable across parameterisations. diff --git a/submission/Revision_reviews-response.md b/submission/Revision_reviews-response.md index ee81dfe..8de7b4a 100644 --- a/submission/Revision_reviews-response.md +++ b/submission/Revision_reviews-response.md @@ -13,7 +13,7 @@ | X | 1 | I think you need to more thoroughly justify why only evaluate with WIS as opposed to or in addition to eg CRPS | Accept \- Add justification in Methods: WIS is proper scoring rule for interval forecasts, matches the quantile forecast format teams submitted, and allows log transformation. | #107 | | | 1 | Presumably all forecast intervals are for negatively binomial distributed count data (cases or death)? Please explicitly state this | Partly accept \- Clarified that observations are of count data (incident cases/deaths) although quantile forecasts don't require teams to specify distributional assumptions. Observations are counts; forecasts are quantile predictions. | #108 | | | 1 | Why log link in the GAMM? Since WIS is already calculated on log scale, it's not intuitive to me why the modelled effects are expected to be multiplicative | Accept \- Add justification: log-link accounts for right-skew in WIS distribution, ensures positive predictions, and multiplicative effects are interpretable as proportional changes in score. Alternatively, could re-run with identity link as sensitivity check. | #109 | -| | 1 | In the QQ plot in your supp results it looks like dispersion of WIS score in the fitted GAMM is a bit off? | Accept \- Examine QQ plots carefully. If concerning: (1) try alternative link function or error distribution, (2) acknowledge limitation explicitly in Methods/Discussion, (3) note that point estimates may be robust even if uncertainty estimates are affected. | #110 | +| X | 1 | In the QQ plot in your supp results it looks like dispersion of WIS score in the fitted GAMM is a bit off? | Accept \- The Gaussian family left deviance residuals with skew 5.8 and kurtosis 77. We compared Gaussian, Gamma and Tweedie error families and adopted a Tweedie family with log link, reducing residual skew to 0.58 and kurtosis to 9.2 and raising deviance explained from 0.29 to 0.38. Gamma fit comparably but did not converge. Reported in Methods and in a new supplementary comparison table. | #110 | | | 1 | it is vague that by model structure you mean mechanistic vs statistical etc, and by forecast target you mean single vs multi country. I think you need to word this sentence to more accurately reflect the language used later on | Accept \- Revise abstract to specify 'model structure (agent-based, mechanistic, semi-mechanistic, statistical)' and 'geographic specificity (single versus multiple countries)'. | #111 | | | 1 | I think it's worth expanding this paragraph to also point out that different models are by design meant to be good at different things | Accept \- Add 2-3 sentences to paragraph around line 69 noting that models serve multiple purposes, forecasting is one use case, and structural differences may be optimized for different objectives. | #112 | | | 1 | I'm not fully convinced that target-specificity is necessarily linked to how many target countries the model predicts to | Partly accept \- Acknowledge this limitation in Discussion. Clarify that we use number of targets as a PROXY for specificity, but recognize it's imperfect. Note that we lack data on actual model customization practices. Soften claims about 'specificity' vs 'number of targets'. | #113 |