Model WIS with a Tweedie family (#159)

CLAUDE.md

@@ -123,9 +123,9 @@ source(here("R", "process-data.R"))
 #    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"),
-          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")

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`

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,37 +72,42 @@ 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 |>
     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))
 
-  # 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,15 +157,18 @@ 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
       )
   }
 
   # 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)
   }

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)
+}

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

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:
 
 $$

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).