Question: REF-op (prod’s tracker) and sqrt constant 0.018 (E16) differ in several settings. How do WIS and coverage change as one setting at a time is changed from one to the other?

Setup:

  • Disease: flu. Base forecaster: the clean windowed_seasonal replay (53 locations, the 55 hub rounds). Horizons h0–h3.
  • Learning truth: exact (re-run at every round on the data published by then; E00), one week of extra delay.
  • Steps 0–3 run with the 2023-24 season present (the HHS-stitched backfill from E13); steps 4–6 and the references are cold.
  • Scores over the two live seasons, each on its own rounds; states only (US excluded).

Each row changes one setting from the row above:

step change config after the change
0 REF-op — sqrt, adaptive mult 0.03 floor 1e-3, window 20, carry; slow term (mult 0.003) warm-started from the 2023-24 burn-in; leak 0.1 on the fast term
1 no warm start slow_init = NULL the slow term starts at zero
2 no slow term lr_slow = NULL single fast term with the leak
3 no leak fast_decay = 0 sqrt, adaptive, mult 0.03, window 20, carry, with the 2023-24 season present
4 no burn-in season drop the 2023-24 rounds as 3, cold: the eta pool starts empty
5 per-100k units scales = population / 100k as 4, on the square root of admissions per 100k
6 sqrt constant 0.018 (E16) lr = 0.018 sqrt per 100k, constant 0.018

The 2023-24 season does two things: with slow_init it sets the slow term’s starting value (step 1 removes this), and its residuals fill the adaptive learning rate’s window, so eta is sized from the first live round (step 4 removes this). Burn-in outcomes never step the offsets. Step 5 changes only units, so that step 6 changes only the learning rate.

Outside the chain: a branch, sqrt 0.018 warm, which is step 6 with the 2023-24 season present and slow_init = "burn_in_quantile" (the warm start of step 0, here on the single term); and REF-op cold (prod’s covid config: step 0 without the warm start or the 2023-24 season), so the state scatters below can pair warm with warm and cold with cold.

Steps 0, 3, 4 and 6 are the cached runs of E13 (REF-op warm), E02 (reference + burn-in), E01/E02 (reference) and E16 (sqrt 0.018), so they are not independent checks of those notebooks.

WIS reduction % = 100 × (base − calibrated) / base, so positive = calibrated WIS lower than base. L1 coverage bias: lower is better. “Month-avg” coverage bias is computed within each calendar month and then averaged.

source(here::here("scripts/calibration/calibration_ws_replay.R"))
wi <- ws_inputs("flu")
fc_burn_in <- bind_rows(ws_burn_in_forecasts("flu"), wi$fc)
# Each config is written as the experiment that first ran it writes it, so cached runs are reused.
refop_warm <- ws_configs(wi$rate_scales, burn_in = TRUE)[["REF-op warm"]]
sqrt_adaptive <- list(
  ref = "paper", transform = "sqrt", lr_args = list(mult = 0.03, floor = 1e-3),
  burn_in_seasons = character(0), slow_init = NULL, lr_window = 20, season_policy = "carry"
)
const_cfg <- ws_configs(wi$rate_scales, burn_in = FALSE)[["E11 constant 0.1"]]
STEPS <- list(
  `0 REF-op` = list(cfg = refop_warm, burn_in = TRUE),
  `1 no warm start` = list(cfg = list(ref = "op", burn_in_seasons = "2023-2024", slow_init = NULL), burn_in = TRUE),
  `2 no slow term` = list(cfg = list(ref = "op", burn_in_seasons = "2023-2024", slow_init = NULL, lr_slow = NULL), burn_in = TRUE),
  `3 no leak` = list(cfg = utils::modifyList(sqrt_adaptive, list(burn_in_seasons = "2023-2024")), burn_in = TRUE),
  `4 no burn-in season` = list(cfg = sqrt_adaptive, burn_in = FALSE),
  `5 per-100k units` = list(cfg = c(sqrt_adaptive, list(scales = wi$rate_scales)), burn_in = FALSE),
  `6 sqrt constant 0.018 (E16)` = list(cfg = utils::modifyList(const_cfg, list(lr = 0.018, transform = "sqrt")), burn_in = FALSE)
)
# Outside the chain: step 6 with a warm start, and REF-op without one (E16's cached run).
REFS <- list(
  `branch: sqrt 0.018 warm` = list(cfg = utils::modifyList(STEPS[[7]]$cfg, list(burn_in_seasons = "2023-2024", slow_init = "burn_in_quantile")), burn_in = TRUE),
  `REF-op cold` = list(cfg = ws_configs(wi$rate_scales, burn_in = FALSE)[["REF-op cold"]], burn_in = FALSE)
)
run <- function(spec) {
  fc <- if (spec$burn_in) fc_burn_in else wi$fc
  do.call(cal_run_cached, c(list(fc, wi$truth), spec$cfg, list(learn = "exact", vintages = wi$vintages, workers = params$workers)))
}
cals <- purrr::map(c(STEPS, REFS), run)
SEASONS <- c(LIVE, "both")
STEP_COLS <- setNames(c("#000000", "#E69F00", "#56B4E9", "#009E73", "#CC79A7", "#0072B2", "#D55E00", "#AA4499", "#999999"), names(c(STEPS, REFS)))
step_m <- purrr::imap(cals, function(cal, nm) {
  fc <- cal$forecasts %>% mutate(season = season_of(reference_date)) %>% filter(location != "US")
  ws_scores(fc) %>%
    select(season, horizon, which, wis, cal_err) %>%
    tidyr::pivot_wider(names_from = which, values_from = c(wis, cal_err)) %>%
    mutate(step = nm, wis_pct = 100 * (wis_base - wis_cal) / wis_base)
}) %>%
  bind_rows() %>%
  mutate(
    step = factor(step, levels = names(c(STEPS, REFS))),
    season = factor(season, levels = SEASONS, labels = c(LIVE, "both seasons")),
    h = factor(paste0("h", horizon), levels = paste0("h", 0:3))
  )
mv <- cal_month_view(cals, seasons = LIVE)
# Five configs, so the month plots stay readable; the tables below them cover every config.
mv_plot_cfgs <- c("0 REF-op", "1 no warm start", "4 no burn-in season", "6 sqrt constant 0.018 (E16)", "branch: sqrt 0.018 warm")
mv_cols <- STEP_COLS

Month view: WIS, coverage, median error

Scores by calendar month of the reference date, one figure per season, states only (US excluded). Rows are the three metrics, columns the horizons:

  • WIS reduction %: 100 × (WIS_base − WIS_cal) / WIS_base within the month, so positive = calibrated WIS lower than base. Grey bars are the month’s share of the season’s base WIS, so a large swing in a month with a small bar costs little overall.
  • L1 coverage bias reduction %: 100 × (bias_base − bias_cal) / bias_base, so positive = calibrated coverage bias lower than base. The bias is computed within the month: the mean over the 23 levels of |share of outcomes at or below the quantile − level|. A season-pooled bias lets under-prediction on the ramp cancel over-prediction after the peak; the monthly one does not. When the base’s bias is already small, a small absolute change shows as a large percentage; the month-average column of the coverage table below gives the absolute values.
  • Median AE reduction %: the same as WIS reduction % for the absolute error of the 0.5 quantile, so positive = calibrated median closer to the outcome.

A month holds one to five rounds of about 52 states, and outcomes in the same round move together, so a single month’s coverage bias is noisy and sits above zero even for a calibrated forecaster. That’s why the plot shows it relative to the base.

mv_plot_month <- mv$month
if (!is.null(get0("mv_plot_cfgs"))) {
  mv_plot_month <- mv_plot_month %>% filter(config %in% mv_plot_cfgs) %>% mutate(config = droplevels(config))
}
for (s in sort(unique(mv_plot_month$season))) {
  cat("\n\n", cal_h(2), " ", s, "\n\n", sep = "")
  print(cal_month_plot(mv_plot_month, s, colours = get0("mv_cols"), linetype = get0("mv_linetype"), colour_by = get0("mv_colour_by") %||% "config"))
}

2024-2025

2025-2026

cat("\n\n", cal_h(2), " By horizon\n\n", sep = "")

By horizon

cat("WIS reduction % (positive = calibrated WIS lower than base), states only; `both` pools the seasons.\n\n")

WIS reduction % (positive = calibrated WIS lower than base), states only; both pools the seasons.

print(knitr::kable(cal_month_wis_table(mv)))
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
0 REF-op +3.2 +1.9 +1.6 +1.6 +2.2 -0.1 -0.0 +0.2 +2.8 +1.1 +0.9 +1.1
1 no warm start +2.0 +1.0 +0.3 -0.4 +1.7 +0.1 +0.4 +0.9 +1.9 +0.6 +0.3 +0.1
2 no slow term +1.8 +0.9 +0.3 -0.2 +1.4 +0.1 +0.4 +0.9 +1.7 +0.6 +0.4 +0.2
3 no leak +1.9 +0.4 -0.6 -1.5 +2.3 -0.7 -0.7 -0.2 +2.0 -0.0 -0.6 -1.0
4 no burn-in season +1.8 +0.0 -1.4 -2.6 +2.2 -0.6 -0.7 -0.2 +1.9 -0.2 -1.1 -1.7
5 per-100k units +1.8 +0.0 -1.4 -2.6 +2.2 -0.6 -0.7 -0.2 +1.9 -0.2 -1.1 -1.7
6 sqrt constant 0.018 (E16) +1.7 +0.6 +0.3 +0.1 +3.0 +1.0 +1.1 +1.5 +2.2 +0.8 +0.6 +0.6
branch: sqrt 0.018 warm +3.4 +2.1 +2.2 +3.3 +3.4 +0.9 +0.8 +1.2 +3.4 +1.6 +1.6 +2.5
REF-op cold +2.2 +0.9 -0.0 -0.8 +1.8 +0.1 +0.3 +0.9 +2.0 +0.6 +0.1 -0.2 
print(knitr::kable(cal_month_cov_table(mv), 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 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
0 REF-op 0.085 / 0.115 0.079 / 0.133 0.077 / 0.175 0.066 / 0.208 0.035 / 0.078 0.032 / 0.111 0.027 / 0.129 0.029 / 0.146 0.060 / 0.095 0.049 / 0.121 0.043 / 0.150 0.038 / 0.175
1 no warm start 0.109 / 0.122 0.110 / 0.139 0.124 / 0.176 0.132 / 0.213 0.059 / 0.078 0.052 / 0.101 0.051 / 0.122 0.049 / 0.128 0.084 / 0.099 0.080 / 0.119 0.087 / 0.147 0.090 / 0.167
2 no slow term 0.113 / 0.124 0.115 / 0.140 0.128 / 0.175 0.140 / 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.101 / 0.170
3 no leak 0.092 / 0.116 0.085 / 0.137 0.090 / 0.173 0.091 / 0.209 0.011 / 0.078 0.010 / 0.105 0.012 / 0.128 0.014 / 0.144 0.050 / 0.096 0.045 / 0.120 0.044 / 0.149 0.040 / 0.174
4 no burn-in season 0.085 / 0.115 0.077 / 0.138 0.078 / 0.174 0.075 / 0.213 0.008 / 0.078 0.010 / 0.105 0.014 / 0.130 0.016 / 0.143 0.045 / 0.095 0.039 / 0.120 0.037 / 0.150 0.032 / 0.176
5 per-100k units 0.085 / 0.115 0.077 / 0.138 0.078 / 0.174 0.075 / 0.213 0.008 / 0.078 0.010 / 0.105 0.014 / 0.130 0.016 / 0.143 0.045 / 0.095 0.039 / 0.120 0.037 / 0.150 0.032 / 0.176
6 sqrt constant 0.018 (E16) 0.099 / 0.115 0.104 / 0.136 0.120 / 0.170 0.137 / 0.204 0.021 / 0.075 0.023 / 0.098 0.024 / 0.119 0.022 / 0.138 0.059 / 0.094 0.062 / 0.116 0.071 / 0.143 0.078 / 0.169
branch: sqrt 0.018 warm 0.076 / 0.115 0.070 / 0.131 0.065 / 0.174 0.057 / 0.218 0.018 / 0.083 0.027 / 0.122 0.035 / 0.156 0.051 / 0.175 0.040 / 0.098 0.034 / 0.126 0.032 / 0.164 0.040 / 0.195
REF-op cold 0.106 / 0.121 0.107 / 0.140 0.117 / 0.175 0.124 / 0.213 0.057 / 0.078 0.050 / 0.100 0.049 / 0.121 0.047 / 0.127 0.081 / 0.098 0.077 / 0.119 0.082 / 0.146 0.085 / 0.167

What each step changes

The change each step makes to the step before it, both seasons pooled, states only. WIS: change in WIS reduction % (positive = the step lowers WIS). Coverage: change in L1 coverage bias, season-pooled and month-average (negative = the step improves coverage).

wis_both <- step_m %>%
  filter(season == "both seasons", step %in% names(STEPS)) %>%
  transmute(step, h, v = wis_pct)
cov_both <- mv$cov %>%
  filter(season == "both") %>%
  transmute(step = config, h = factor(paste0("h", horizon), levels = paste0("h", 0:3)), pooled = pooled_cal, mavg = mavg_cal)
delta <- wis_both %>%
  left_join(cov_both, by = c("step", "h")) %>%
  arrange(h, step) %>%
  group_by(h) %>%
  mutate(d_wis = v - lag(v), d_pooled = pooled - lag(pooled), d_mavg = mavg - lag(mavg)) %>%
  ungroup() %>%
  filter(step != names(STEPS)[1])
knitr::kable(
  delta %>%
    transmute(step, h, cell = sprintf("%+.1f / %+.3f / %+.3f", d_wis, d_pooled, d_mavg)) %>%
    tidyr::pivot_wider(names_from = h, values_from = cell) %>%
    arrange(step),
  caption = "Both seasons, states only: change from the previous step in WIS reduction % / season-pooled L1 coverage bias / month-avg L1 coverage bias"
)
Both seasons, states only: change from the previous step in WIS reduction % / season-pooled L1 coverage bias / month-avg L1 coverage bias
step h0 h1 h2 h3
1 no warm start -0.9 / +0.024 / +0.003 -0.5 / +0.032 / -0.002 -0.6 / +0.044 / -0.003 -1.0 / +0.052 / -0.008
2 no slow term -0.2 / +0.006 / +0.002 -0.0 / +0.006 / +0.002 +0.0 / +0.008 / +0.002 +0.1 / +0.011 / +0.003
3 no leak +0.4 / -0.040 / -0.006 -0.6 / -0.042 / -0.001 -1.0 / -0.051 / -0.000 -1.2 / -0.060 / +0.004
4 no burn-in season -0.1 / -0.005 / -0.000 -0.2 / -0.006 / +0.000 -0.5 / -0.007 / +0.001 -0.7 / -0.009 / +0.001
5 per-100k units +0.0 / +0.000 / +0.000 +0.0 / +0.000 / +0.000 +0.0 / +0.000 / +0.000 +0.0 / +0.000 / +0.000
6 sqrt constant 0.018 (E16) +0.3 / +0.014 / -0.001 +1.0 / +0.023 / -0.005 +1.7 / +0.034 / -0.007 +2.3 / +0.047 / -0.007

State by state: REF-op against sqrt 0.018, matched on the warm start

Per state, pooled over h0–h3: WIS reduction % of REF-op against sqrt 0.018 (above the line, REF-op is better). Three pairs: the chain’s ends (REF-op warm vs sqrt 0.018 cold), then warm vs warm and cold vs cold. Cold configs play the base through Nov–Dec 2024 (35–40% of that season’s base WIS), so only the matched pairs compare the trackers rather than the start.

PAIRS <- list(
  `chain ends: REF-op warm vs sqrt 0.018 cold` = c("0 REF-op", "6 sqrt constant 0.018 (E16)"),
  `warm vs warm` = c("0 REF-op", "branch: sqrt 0.018 warm"),
  `cold vs cold` = c("REF-op cold", "6 sqrt constant 0.018 (E16)")
)
state_pct <- cal_state_wis(cals[unique(unlist(PAIRS))]) %>%
  group_by(config, season, location) %>%
  summarize(pct = 100 * (sum(wis_base) - sum(wis_cal)) / sum(wis_base), .groups = "drop")
pair_pts <- purrr::imap(PAIRS, function(p, nm) {
  state_pct %>%
    filter(config %in% p) %>%
    mutate(which = ifelse(config == p[1], "refop", "sqrt")) %>%
    select(-config) %>%
    tidyr::pivot_wider(names_from = which, values_from = pct) %>%
    mutate(pair = nm)
}) %>%
  bind_rows() %>%
  mutate(pair = factor(pair, levels = names(PAIRS)))
ggplot(pair_pts, aes(sqrt, refop)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey50") +
  geom_hline(yintercept = 0, colour = "grey85") +
  geom_vline(xintercept = 0, colour = "grey85") +
  geom_point(alpha = 0.8) +
  facet_grid(season ~ pair) +
  labs(x = "sqrt 0.018: WIS reduction %, h0–h3", y = "REF-op: WIS reduction %, h0–h3")

States where REF-op has the lower WIS, out of the states scored, per pair and season.

knitr::kable(
  pair_pts %>%
    group_by(pair, season) %>%
    summarize(refop_better = sum(refop > sqrt), states = n(), .groups = "drop") %>%
    transmute(pair, season, cell = sprintf("%d / %d", refop_better, states)) %>%
    tidyr::pivot_wider(names_from = season, values_from = cell),
  caption = "States (of those scored) where REF-op beats sqrt 0.018, h0–h3 pooled"
)
States (of those scored) where REF-op beats sqrt 0.018, h0–h3 pooled
pair 2024-2025 2025-2026
chain ends: REF-op warm vs sqrt 0.018 cold 45 / 52 13 / 52
warm vs warm 23 / 52 15 / 52
cold vs cold 17 / 52 15 / 52