Question: does calibration change WIS and coverage the same way for the submitted ensemble as for one fixed forecaster?

Setup. Two forecasters, everything else shared:

  • A: the submitted flu ensemble, read from the hub as in E00 and E03–E05.
  • B: the clean windowed_seasonal replay (E12–E17), with the 2023-24 HHS backfill (E13) as its burn-in season.
rows identical (round, location, horizon, level) rows for A and B: hub submission rounds, submissions excluded from the hub reads dropped from both
seasons 2023-24 burn-in (27 rounds), 2024-25 (27) and 2025-26 (28) scored
locations all 53 learned; scores are states only
horizons h0–h3 (the ensemble’s h−1 is a different model)
learning truth NHSN as published each week, re-run from scratch each week, one week of extra revision delay
scoring truth the store’s NHSN truth

Trackers (section 2), all sqrt scale, adaptive learning rate (mult 0.03, floor 0.001, window 20, carry):

  • single: REF-paper at sqrt.
  • leaky: single with the fast offset shrinking 10% per round.
  • operating: REF-op (leaky plus a slow offset and a warm start from the burn-in season’s errors).

This rebuilds the original E07 (a script since deleted), which learned from finalized data, used the replay with prod’s data substitutions and included h−1.

“WIS reduction %” is 100 × (WIS_base − WIS_cal) / WIS_base, so positive = calibrated WIS lower than base. “L1 coverage bias” (lower is better) is the mean over the 23 levels of |empirical coverage − level|.

source(here::here("scripts/calibration/calibration_ws_replay.R"))
wi <- ws_inputs("flu")
hub_fc <- suppressMessages(hub_read_forecasts(hub_dir = here::here(params$hub_dir), target = HUB_FLU_TARGET)) %>%
  filter(horizon %in% 0:3)
ws_fc <- bind_rows(ws_burn_in_forecasts("flu"), wi$fc)
# Keep only (round, location, horizon, level) rows both forecasters have.
key <- c("reference_date", "location", "horizon", "level")
common <- inner_join(distinct(hub_fc, across(all_of(key))), distinct(ws_fc, across(all_of(key))), by = key)
FCS <- list(
  `A ensemble` = semi_join(hub_fc, common, by = key),
  `B windowed_seasonal` = semi_join(ws_fc, common, by = key)
)
sqrt_single <- list(ref = "paper", transform = "sqrt", lr_args = list(mult = 0.03, floor = 1e-3), burn_in_seasons = "2023-2024")
TRACKERS <- list(
  single = sqrt_single,
  leaky = c(sqrt_single, list(fast_decay = 0.1)),
  operating = list(ref = "op")
)
grid <- tidyr::expand_grid(forecaster = names(FCS), tracker = names(TRACKERS))
cals <- purrr::pmap(grid, function(forecaster, tracker) {
  do.call(cal_run_cached, c(list(FCS[[forecaster]], wi$truth), TRACKERS[[tracker]], list(learn = "exact", vintages = wi$vintages, workers = params$workers)))
})
purrr::imap(FCS, \(fc, nm) fc %>% mutate(season = season_of(reference_date)) %>% count(season, name = nm)) %>%
  purrr::reduce(full_join, by = "season") %>%
  left_join(FCS[[1]] %>% mutate(season = season_of(reference_date)) %>% group_by(season) %>% summarize(rounds = n_distinct(reference_date)), by = "season") %>%
  knitr::kable()
season A ensemble B windowed_seasonal rounds
2023-2024 123372 123372 27
2024-2025 131376 131376 27
2025-2026 136528 136528 28 
GEOS <- c(states = "states only (US excluded)")
SEASONS <- c(LIVE, "both")
score <- function(cal) {
  fc <- cal$forecasts %>% mutate(season = season_of(reference_date))
  purrr::imap(WS_GEOS[names(GEOS)], function(keep, g) ws_scores(fc %>% filter(keep(location))) %>% mutate(geos = g)) %>%
    bind_rows() %>%
    select(geos, season, horizon, which, wis, cal_err) %>%
    tidyr::pivot_wider(names_from = which, values_from = c(wis, cal_err)) %>%
    mutate(wis_pct = 100 * (wis_base - wis_cal) / wis_base)
}
m <- purrr::map2(cals, seq_len(nrow(grid)), \(cal, i) score(cal) %>% mutate(grid[i, ])) %>%
  bind_rows() %>%
  mutate(
    forecaster = factor(forecaster, levels = names(FCS)), tracker = factor(tracker, levels = names(TRACKERS)),
    geos = factor(geos, levels = names(GEOS), labels = GEOS),
    season = factor(season, levels = SEASONS, labels = c(LIVE, "both seasons")),
    h = factor(paste0("h", horizon), levels = paste0("h", 0:3))
  )
FC_COLS <- c(`A ensemble` = "#0072B2", `B windowed_seasonal` = "#D55E00")

1. The two bases

Base WIS

Uncalibrated WIS (2 × mean pinball, in admissions) on identical rows, and the ratio B/A (below 1: windowed_seasonal is better).

m %>%
  filter(tracker == "single") %>%
  select(season, h, forecaster, wis_base) %>%
  tidyr::pivot_wider(names_from = forecaster, values_from = wis_base) %>%
  transmute(season, h, A = .data[[names(FCS)[1]]], B = .data[[names(FCS)[2]]], `B/A` = B / A) %>%
  arrange(season, h) %>%
  knitr::kable(digits = c(0, 0, 1, 1, 3), caption = "Base WIS, A ensemble vs B windowed_seasonal, states only")
Base WIS, A ensemble vs B windowed_seasonal, states only
season h A B B/A
2024-2025 h0 64.6 60.4 0.934
2024-2025 h1 89.2 88.3 0.991
2024-2025 h2 110.5 112.0 1.013
2024-2025 h3 136.7 138.2 1.011
2025-2026 h0 28.6 38.0 1.328
2025-2026 h1 49.3 60.0 1.218
2025-2026 h2 67.2 74.6 1.111
2025-2026 h3 78.8 82.4 1.046
both seasons h0 46.3 49.0 1.058
both seasons h1 68.9 73.9 1.073
both seasons h2 88.5 92.9 1.051
both seasons h3 107.2 109.8 1.024

Base bias by month

Share of outcomes above the base median and L1 coverage bias, by calendar month of the reference date, states only.

base_fc <- purrr::map(which(grid$tracker == "single"), \(i) cals[[i]]$forecasts %>% mutate(forecaster = grid$forecaster[i])) %>%
  bind_rows() %>%
  mutate(season = season_of(reference_date)) %>%
  filter(season %in% LIVE, !is.na(truth), !is.na(value_base))
month_base <- function(d) {
  d %>%
    group_by(across(c(any_of(c("forecaster", "season", "horizon", "month")), level))) %>%
    summarize(gap = mean(truth <= value_base) - first(level), above = mean(truth > value_base), n = n(), .groups = "drop") %>%
    group_by(across(any_of(c("forecaster", "season", "horizon", "month")))) %>%
    summarize(above = above[level == 0.5], cal_err = mean(abs(gap)), n = first(n), .groups = "drop")
}
base_fc_states <- base_fc %>%
  filter(location != "US") %>%
  mutate(month = factor(format(reference_date, "%b"), levels = month.abb[c(8:12, 1:7)]))
mb <- month_base(base_fc_states) %>% mutate(h = factor(paste0("h", horizon), levels = paste0("h", 0:3)))
mb %>%
  tidyr::pivot_longer(c(above, cal_err), names_to = "metric") %>%
  mutate(metric = factor(metric, levels = c("above", "cal_err"), labels = c("share above median", "L1 coverage bias (lower is better)"))) %>%
  ggplot(aes(month, value, colour = forecaster, group = forecaster)) +
  geom_hline(data = data.frame(metric = factor("share above median", levels = c("share above median", "L1 coverage bias (lower is better)")), y = 0.5), aes(yintercept = y), colour = "grey50") +
  geom_line() +
  geom_point(size = 1.1) +
  facet_grid(metric + season ~ h, scales = "free", space = "free_x", labeller = labeller(metric = label_wrap_gen(20))) +
  scale_colour_manual(values = FC_COLS, name = NULL) +
  labs(x = "month of reference date", y = NULL) +
  theme(legend.position = "bottom")

2. What calibration does to each

WIS reduction % (positive = calibrated WIS lower than base) and L1 coverage bias (lower is better), states only.

mv_trk <- cal_month_view(setNames(cals, paste(grid$forecaster, grid$tracker, sep = ", ")), seasons = LIVE)
knitr::kable(cal_month_wis_table(mv_trk))
config 2024-2025 h0 2024-2025 h1 2024-2025 h2 2024-2025 h3 2025-2026 h0 2025-2026 h1 2025-2026 h2 2025-2026 h3 both h0 both h1 both h2 both h3
A ensemble, single +2.2 +1.1 -0.2 -1.5 +2.5 -2.0 -3.1 -2.4 +2.3 -0.0 -1.3 -1.8
A ensemble, leaky +1.7 +1.0 +0.1 -0.5 +3.2 +0.4 -0.1 +0.2 +2.2 +0.8 -0.0 -0.3
A ensemble, operating +3.7 +2.9 +2.8 +4.3 +4.3 +0.1 -0.7 -0.9 +3.9 +1.9 +1.5 +2.3
B windowed_seasonal, single +1.9 +0.4 -0.6 -1.5 +2.2 -0.7 -0.7 -0.2 +2.0 -0.0 -0.6 -1.0
B windowed_seasonal, leaky +1.9 +1.0 +0.4 -0.2 +1.4 +0.1 +0.4 +0.9 +1.7 +0.6 +0.4 +0.2
B windowed_seasonal, operating +3.2 +1.9 +1.4 +1.7 +2.2 -0.1 -0.0 +0.3 +2.8 +1.1 +0.8 +1.2 
knitr::kable(cal_month_cov_table(mv_trk), caption = CAL_MONTH_COV_CAPTION)
L1 coverage bias (lower is better), states only, as season-pooled / month-avg. Season-pooled scores all of a season’s rounds together. Month-avg is the mean over calendar months of the bias computed within each month, every month weighted equally (‘both’: every season-month weighted equally).
config 2024-2025 h0 2024-2025 h1 2024-2025 h2 2024-2025 h3 2025-2026 h0 2025-2026 h1 2025-2026 h2 2025-2026 h3 both h0 both h1 both h2 both h3
base (A ensemble, single) 0.131 / 0.146 0.135 / 0.152 0.146 / 0.169 0.159 / 0.190 0.078 / 0.089 0.055 / 0.098 0.050 / 0.108 0.046 / 0.117 0.104 / 0.116 0.094 / 0.123 0.097 / 0.137 0.101 / 0.151
base (B windowed_seasonal, single) 0.127 / 0.130 0.129 / 0.144 0.148 / 0.179 0.169 / 0.211 0.090 / 0.097 0.081 / 0.112 0.084 / 0.130 0.086 / 0.146 0.108 / 0.112 0.105 / 0.127 0.115 / 0.153 0.127 / 0.177
A ensemble, single 0.091 / 0.117 0.088 / 0.121 0.088 / 0.143 0.081 / 0.181 0.034 / 0.087 0.038 / 0.114 0.034 / 0.118 0.030 / 0.123 0.031 / 0.101 0.026 / 0.117 0.027 / 0.130 0.025 / 0.150
A ensemble, leaky 0.112 / 0.132 0.116 / 0.139 0.125 / 0.159 0.130 / 0.180 0.052 / 0.076 0.035 / 0.095 0.034 / 0.108 0.033 / 0.110 0.081 / 0.102 0.074 / 0.115 0.078 / 0.132 0.081 / 0.143
A ensemble, operating 0.084 / 0.116 0.081 / 0.121 0.085 / 0.138 0.072 / 0.163 0.012 / 0.067 0.008 / 0.100 0.018 / 0.111 0.030 / 0.116 0.047 / 0.090 0.036 / 0.110 0.039 / 0.124 0.035 / 0.138
B windowed_seasonal, single 0.090 / 0.116 0.084 / 0.137 0.088 / 0.173 0.089 / 0.209 0.010 / 0.078 0.010 / 0.105 0.012 / 0.128 0.013 / 0.142 0.049 / 0.096 0.044 / 0.120 0.043 / 0.149 0.040 / 0.173
B windowed_seasonal, leaky 0.112 / 0.124 0.114 / 0.140 0.128 / 0.175 0.139 / 0.212 0.068 / 0.081 0.060 / 0.104 0.063 / 0.126 0.063 / 0.134 0.090 / 0.101 0.087 / 0.121 0.095 / 0.149 0.100 / 0.171
B windowed_seasonal, operating 0.087 / 0.116 0.081 / 0.134 0.082 / 0.175 0.066 / 0.208 0.037 / 0.078 0.033 / 0.111 0.030 / 0.128 0.027 / 0.145 0.061 / 0.096 0.051 / 0.121 0.047 / 0.150 0.037 / 0.174