Overall Comparison of Hospitalization Forecasters

# params <- list(
#   forecaster_parameters = tar_read(forecaster_parameter_combinations),
#   forecasts = tar_read(joined_forecasts),
#   scores = tar_read(joined_scores),
#   truth_data = tar_read(hhs_evaluation_data),
#   disease = "flu"
# )

if (params$disease == "flu") {
  base_forecaster_name <- "FluSight-baseline"
  ensemble_forecaster_name <- "FluSight-ensemble"
} else {
  base_forecaster_name <- "COVIDhub-baseline"
  ensemble_forecaster_name <- "COVIDhub_CDC-ensemble"
}

# Select forecasters for this notebook
outside_forecasters <- c(base_forecaster_name, ensemble_forecaster_name)
id_to_forecaster <- params$forecaster_parameters %>%
  purrr::imap(\(x, idx) {
    x %>%
      mutate(family_name = idx) %>%
      select(id, forecaster, family_name)
  }) %>%
  bind_rows() %>%
  rename(forecaster_function = forecaster)

overall_rating <- params$scores %>%
  summarize(mean_score = mean(wis), .by = forecaster) %>%
  arrange(mean_score) %>%
  filter(mean_score > 20) %>%
  left_join(id_to_forecaster, by = join_by(forecaster == id)) %>%
  mutate(
    family_name = if_else(is.na(family_name), "external", family_name),
    forecaster_function = if_else(is.na(forecaster_function), "external", forecaster_function)
  )
best_in_class <- overall_rating %>%
  slice_min(mean_score, by = family_name) %>%
  pull(forecaster)
forecasters_to_keep <- c(best_in_class, outside_forecasters) %>% unique()
# Load forecasts and scores
scores <- params$scores %>%
  filter(forecaster %in% forecasters_to_keep)
forecast_dates <- scores %>%
  pull(forecast_date) %>%
  unique()
forecasters <- best_in_class
aheads <- scores %>%
  pull(ahead) %>%
  unique()

# Define aggregation functions
Mean <- function(x) mean(x, na.rm = TRUE)
GeoMean <- function(x, offset = 0) exp(Mean(log(x + offset)))

Notebook Information

Best examples from each family of forecasters.

Disease

covid

Forecast Dates

2023-11-11, 2023-11-18, 2023-11-25, 2023-12-02, 2023-12-09, 2023-12-16, 2023-12-23, 2023-12-30, 2024-01-06, 2024-01-13, 2024-01-20, 2024-01-27, 2024-02-03, 2024-02-10, 2024-02-17, 2024-02-24, 2024-03-02, 2024-03-09, 2024-03-16, 2024-03-23, 2024-03-30, 2024-04-06, 2024-04-13, 2024-04-20, 2024-04-27

Forecaster Parameter Mapping and Overall Scores

The table is sorted by ascending WIS and contains all the forecasters in this notebook.

# Display the table
if (params$disease == "flu") {
  ignore_keys <- c("forecaster", "keys_to_ignore", "pop_scaling")
} else {
  ignore_keys <- c("forecaster", "keys_to_ignore")
}

param_table <- params$forecaster_parameters %>%
  map(\(x) x %>% filter(id %in% forecasters_to_keep)) %>%
  bind_rows() %>%
  left_join(id_to_forecaster, by = join_by(id)) %>%
  select(-any_of(ignore_keys)) %>%
  full_join(
    scores %>%
      group_by(forecaster) %>%
      summarize(
        mean_ae = round(Mean(ae), 2),
        geomean_ae = round(GeoMean(ae), 2),
        mean_wis = round(Mean(wis), 2),
        geomean_wis = round(GeoMean(wis), 2),
        mean_coverage_90 = round(Mean(coverage_90), 2),
      ) %>%
      rename(id = forecaster)
  ) %>%
  arrange(mean_ae) %>%
  relocate(family = family_name, id, mean_ae, geomean_ae, mean_wis, geomean_wis, mean_coverage_90)
datatable(param_table)

Show entries

Search:

	family	id	mean_ae	geomean_ae	mean_wis	geomean_wis	mean_coverage_90	trainer	lags	pop_scaling	n_training	outcome	extra_sources	scale_method	seasonal_method	center_method	nonlin_method	model_used	filter_agg_level	drop_non_seasons	quantiles_by_geo	aheads	residual_tail	residual_center	forecaster_function

	family	id	mean_ae	geomean_ae	mean_wis	geomean_wis	mean_coverage_90	trainer	lags	pop_scaling	outcome	extra_sources	scale_method	seasonal_method	center_method	nonlin_method	model_used	filter_agg_level	drop_non_seasons	quantiles_by_geo	aheads	residual_tail	residual_center	forecaster_function
1	scaled_pop_exogenous	emasculative.jackal	97.56	34.49	60.36	25.64	0.81	quantreg	0,7,14,21,0,7,14	false	hhs	nssp	quantile											scaled_pop
2	climate_linear	unwasteful.hammerheadbird	113.98	44.18	72.86	34.11	0.74				hhs		quantile		median	quart_root	climate_linear	state	true	true	0,7,14,21,28	0.7	0.127	climate_linear_ensembled
3	scaled_pop_main	forgivable.invisiblerail	121.64	54.07	79.23	41.95	0.84	quantreg	0,7,14	false	hhs													scaled_pop
4		COVIDhub_CDC-ensemble	130.06	0	81.19	39.5	0.99
5	scaled_pop_season	unlunar.titi	131.31	48.65	85.42	37.94	0.81	quantreg	0,7,14,21	false	hhs			covid										scaled_pop_seasonal
6	scaled_pop_season	porcelain.americanbulldog	131.31	48.65	85.42	37.94	0.81	quantreg	0,7,14,21	false	hhs			window										scaled_pop_seasonal
7	scaled_pop_season	islandlike.cow	131.31	48.65	85.42	37.94	0.81	quantreg	0,7,14,21	false	hhs			covid,window										scaled_pop_seasonal
8	scaled_pop_season	cowardly.americanrobin	131.31	48.65	85.42	37.94	0.81	quantreg	0,7,14,21	false	hhs			climatological										scaled_pop_seasonal
9	scaled_pop_season	interchurch.opossum	131.31	48.65	85.42	37.94	0.81	quantreg	0,7,14,21	false	hhs			climatological,window										scaled_pop_seasonal
10	flatline_forecaster	specified.turkey	134.71	0	97.68	56.44	0.96				hhs													flatline_fc

Showing 1 to 10 of 11 entries

Previous1 2Next

Score Plots

The WIS plots are relative to the baseline. These plots are aggregated with the geometric mean.
The absolute error plots are also relative to the baseline and are aggregated with the arithmetic mean.
The dashed black line in all plots is the ensemble.

var <- "wis"
group_cols <- c("forecaster", "forecast_date", "ahead")

# Aggregate metric across groups
df <- scores %>%
  select(all_of(c(group_cols, var))) %>%
  drop_na(!!sym(var)) %>%
  summarize(!!var := GeoMean(!!sym(var)), .by = all_of(group_cols)) %>%
  filter(ahead >= 0)

# Make sure we don't divide by zero
if (
  df %>%
    filter(forecaster == base_forecaster_name & near(!!sym(var), 0)) %>%
    nrow() > 0
) {
  warning("scale_by_forecaster will divide by zero in column ", var)
}

# Normalize the metric by the baseline
normalized_df <- df %>%
  pivot_wider(names_from = forecaster, names_prefix = var, values_from = !!sym(var)) %>%
  mutate(across(starts_with(var), ~ .x / !!sym(paste0(var, base_forecaster_name)))) %>%
  pivot_longer(cols = starts_with(var), names_to = "forecaster", values_to = var) %>%
  mutate(forecaster = stringr::str_remove(forecaster, var)) %>%
  filter(forecaster != base_forecaster_name)

facets.label <- str_glue("{aheads} days ahead")
names(facets.label) <- aheads
subtitle <- sprintf(
  "Forecasts made over %s to %s",
  format(min(forecast_dates), "%B %d, %Y"),
  format(max(forecast_dates), "%B %d, %Y")
)
p <- ggplot(
  normalized_df,
  aes(x = forecast_date, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_hline(yintercept = 1, linetype = 1, color = "black") +
  facet_grid(rows = vars(ahead)) +
  facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead = facets.label)) +
  scale_y_log10(n.breaks = 10) +
  scale_color_discrete() +
  guides(color = guide_legend(ncol = 2)) +
  labs(title = subtitle, x = "Forecast Dates", y = "Geometric Mean WIS")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))

var <- "wis"
group_cols <- c("forecaster", "ahead")

# Aggregate metric across groups
df <- scores %>%
  select(all_of(c(group_cols, var))) %>%
  drop_na(!!sym(var)) %>%
  summarize(!!var := GeoMean(!!sym(var)), .by = all_of(group_cols)) %>%
  filter(ahead >= 0)

# Make sure we don't divide by zero
if (df %>% filter(forecaster == base_forecaster_name & near(!!sym(var), 0)) %>% nrow() > 0) {
  warning("scale_by_forecaster will divide by zero in column ", var)
}

# Normalize the metric by the baseline
normalized_df <- df %>%
  pivot_wider(names_from = forecaster, names_prefix = var, values_from = !!sym(var)) %>%
  mutate(across(starts_with(var), ~ .x / !!sym(paste0(var, base_forecaster_name)))) %>%
  pivot_longer(cols = starts_with(var), names_to = "forecaster", values_to = var) %>%
  mutate(forecaster = stringr::str_remove(forecaster, var)) %>%
  filter(forecaster != base_forecaster_name)

subtitle <- sprintf(
  "Forecasts made over %s to %s",
  format(min(forecast_dates), "%B %d, %Y"),
  format(max(forecast_dates), "%B %d, %Y")
)
p <- ggplot(
  normalized_df,
  aes(x = ahead, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_hline(yintercept = 1, linetype = 1, color = "black") +
  scale_y_log10(n.breaks = 20) +
  scale_color_discrete() +
  guides(color = guide_legend(ncol = 2)) +
  labs(title = subtitle, x = "Days ahead", y = "Geometric Mean WIS")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))

var <- "ae"
group_cols <- c("forecaster", "forecast_date", "ahead")

# Aggregate metric across groups
df <- scores %>%
  select(all_of(c(group_cols, var))) %>%
  drop_na(!!sym(var)) %>%
  summarize(!!var := Mean(!!sym(var)), .by = all_of(group_cols)) %>%
  filter(ahead >= 0)

# Make sure we don't divide by zero
if (
  df %>%
    filter(forecaster == base_forecaster_name & near(!!sym(var), 0)) %>%
    nrow() > 0
) {
  warning("scale_by_forecaster will divide by zero in column ", var)
}

# Normalize the metric by the baseline
normalized_df <- df %>%
  pivot_wider(names_from = forecaster, names_prefix = var, values_from = !!sym(var)) %>%
  mutate(across(starts_with(var), ~ .x / !!sym(paste0(var, base_forecaster_name)))) %>%
  pivot_longer(cols = starts_with(var), names_to = "forecaster", values_to = var) %>%
  mutate(forecaster = stringr::str_remove(forecaster, var)) %>%
  filter(forecaster != base_forecaster_name)

facets.label <- str_glue("{aheads} days ahead")
names(facets.label) <- aheads
subtitle <- sprintf(
  "Forecasts made over %s to %s",
  format(min(forecast_dates), "%B %d, %Y"),
  format(max(forecast_dates), "%B %d, %Y")
)
p <- ggplot(
  normalized_df,
  aes(x = forecast_date, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_hline(yintercept = 1, linetype = 1, color = "black") +
  facet_grid(rows = vars(ahead)) +
  facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead = facets.label)) +
  scale_y_log10(n.breaks = 10) +
  scale_color_discrete() +
  guides(color = guide_legend(ncol = 2)) +
  labs(title = subtitle, x = "Forecast Dates", y = "Arithmetic Mean Absolute Error")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))

var <- "ae"
group_cols <- c("forecaster", "ahead")

# Aggregate metric across groups
df <- scores %>%
  select(all_of(c(group_cols, var))) %>%
  drop_na(!!sym(var)) %>%
  summarize(!!var := Mean(!!sym(var)), .by = all_of(group_cols)) %>%
  filter(ahead >= 0)

# Make sure we don't divide by zero
if (df %>% filter(forecaster == base_forecaster_name & near(!!sym(var), 0)) %>% nrow() > 0) {
  warning("scale_by_forecaster will divide by zero in column ", var)
}

# Normalize the metric by the baseline
normalized_df <- df %>%
  pivot_wider(names_from = forecaster, names_prefix = var, values_from = !!sym(var)) %>%
  mutate(across(starts_with(var), ~ .x / !!sym(paste0(var, base_forecaster_name)))) %>%
  pivot_longer(cols = starts_with(var), names_to = "forecaster", values_to = var) %>%
  mutate(forecaster = stringr::str_remove(forecaster, var)) %>%
  filter(forecaster != base_forecaster_name)

subtitle <- sprintf(
  "Forecasts made over %s to %s",
  format(min(forecast_dates), "%B %d, %Y"),
  format(max(forecast_dates), "%B %d, %Y")
)
p <- ggplot(
  normalized_df,
  aes(x = ahead, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_hline(yintercept = 1, linetype = 1, color = "black") +
  scale_y_log10(n.breaks = 25) +
  scale_color_discrete() +
  guides(color = guide_legend(ncol = 2)) +
  labs(title = subtitle, x = "Days ahead", y = "Arithmetic Mean Absolute Error")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))

var <- "coverage_90"
group_cols <- c("forecaster", "forecast_date", "ahead")

# Aggregate metric across groups
df <- scores %>%
  drop_na(!!sym(var)) %>%
  summarize(!!var := Mean(!!sym(var)), .by = all_of(group_cols)) %>%
  filter(ahead >= 0)

facets.label <- str_glue("{aheads} days ahead")
names(facets.label) <- aheads
subtitle <- sprintf(
  "Forecasts made over %s to %s",
  format(min(forecast_dates), "%B %d, %Y"),
  format(max(forecast_dates), "%B %d, %Y")
)
p <- ggplot(
  df,
  aes(x = forecast_date, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_hline(yintercept = .8, linetype = 1, color = "black") +
  facet_grid(rows = vars(ahead)) +
  facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead = facets.label)) +
  scale_color_discrete() +
  guides(color = guide_legend(ncol = 2)) +
  labs(title = subtitle, x = "Forecast Dates", y = "Arithmetic Mean 90% Coverage")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))

var <- "coverage_90"
id_cols <- c("forecaster", "ahead")

# Aggregate metric across groups
df <- scores %>%
  select(all_of(c(id_cols, var))) %>%
  drop_na(!!sym(var)) %>%
  summarize(!!var := Mean(!!sym(var)), .by = all_of(id_cols)) %>%
  filter(ahead >= 0)

subtitle <- sprintf(
  "Forecasts made over %s to %s",
  format(min(forecast_dates), "%B %d, %Y"),
  format(max(forecast_dates), "%B %d, %Y")
)
p <- ggplot(
  df,
  aes(x = ahead, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_hline(yintercept = .8, linetype = 1, color = "black") +
  scale_color_discrete() +
  guides(color = guide_legend(ncol = 2)) +
  labs(title = subtitle, x = "Days ahead", y = "Arithmetic Mean 90% Coverage")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))

Fan plots

Fan plots showing the 80% prediction intervals for the forecasts made by the CMU forecasters and the outside forecasters. The black line is the truth data.

plot_dates <- forecast_dates[seq(1, length(forecast_dates), by = 4)]

# We plot a subset of the dates and geos for the fan plot
geo_vals <- c("ca", "fl", "pa", "tx", "ny")
forecast_subset <- params$forecasts %>%
  filter(
    forecaster %in% forecasters_to_keep,
    geo_value %in% geo_vals,
    forecast_date %in% plot_dates
  ) %>%
  mutate(quantile = as.character(quantile)) %>%
  pivot_wider(names_from = "quantile", values_from = "prediction") %>%
  mutate(ahead = as.numeric(target_end_date - forecast_date))

p <- ggplot(
  data = forecast_subset,
  aes(x = target_end_date, group = forecast_date)
) +
  geom_ribbon(aes(ymin = `0.1`, ymax = `0.9`, fill = forecast_date), alpha = 0.3) +
  geom_line(aes(y = `0.5`, color = forecast_date), linetype = 2L) +
  geom_point(aes(y = `0.5`, color = forecast_date), size = 0.75) +
  geom_line(
    data = params$truth_data %>%
      filter(
        geo_value %in% geo_vals,
        target_end_date >= min(plot_dates),
        target_end_date <= max(plot_dates) + 5 * 7
      ),
    aes(x = target_end_date, y = true_value),
    inherit.aes = FALSE, na.rm = TRUE,
    color = "black", linetype = 1
  ) +
  scale_y_continuous(n.breaks = 15) +
  facet_grid(factor(forecaster, levels = param_table$id) ~ geo_value, scales = "free") +
  labs(x = "Reference Date", y = "Weekly Sums of Hospitalizations", title = "Monthly Forecasts and Truth Data") +
  theme(legend.position = "none")

ggplotly(p, tooltip = "text", height = 300 * length(param_table$id), width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))

# We plot a subset of the dates and geos for the fan plot
geo_vals <- c("mn", "co", "ky", "md")
forecast_subset <- params$forecasts %>%
  filter(
    forecaster %in% forecasters_to_keep,
    geo_value %in% geo_vals,
    forecast_date %in% plot_dates
  ) %>%
  mutate(quantile = as.character(quantile)) %>%
  pivot_wider(names_from = "quantile", values_from = "prediction") %>%
  mutate(ahead = as.numeric(target_end_date - forecast_date))

p <- ggplot(
  data = forecast_subset,
  aes(x = target_end_date, group = forecast_date)
) +
  geom_ribbon(aes(ymin = `0.1`, ymax = `0.9`, fill = forecast_date), alpha = 0.3) +
  geom_line(aes(y = `0.5`, color = forecast_date), linetype = 2L) +
  geom_point(aes(y = `0.5`, color = forecast_date), size = 0.75) +
  geom_line(
    data = params$truth_data %>%
      filter(
        geo_value %in% geo_vals,
        target_end_date >= min(plot_dates),
        target_end_date <= max(plot_dates) + 5 * 7
      ),
    aes(x = target_end_date, y = true_value),
    inherit.aes = FALSE, na.rm = TRUE,
    color = "black", linetype = 1
  ) +
  scale_y_continuous(n.breaks = 15) +
  facet_grid(factor(forecaster, levels = param_table$id) ~ geo_value, scales = "free") +
  labs(x = "Reference Date", y = "Weekly Sums of Hospitalizations", title = "Monthly Forecasts and Truth Data") +
  theme(legend.position = "none")

ggplotly(p, tooltip = "text", height = 300 * length(param_table$id), width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))

# We plot a subset of the dates and geos for the fan plot
geo_vals <- c("nv", "vt", "de", "ri", "nm")
forecast_subset <- params$forecasts %>%
  filter(
    forecaster %in% forecasters_to_keep,
    geo_value %in% geo_vals,
    forecast_date %in% plot_dates
  ) %>%
  mutate(quantile = as.character(quantile)) %>%
  pivot_wider(names_from = "quantile", values_from = "prediction") %>%
  mutate(ahead = as.numeric(target_end_date - forecast_date))

p <- ggplot(
  data = forecast_subset,
  aes(x = target_end_date, group = forecast_date)
) +
  geom_ribbon(aes(ymin = `0.1`, ymax = `0.9`, fill = forecast_date), alpha = 0.3) +
  geom_line(aes(y = `0.5`, color = forecast_date), linetype = 2L) +
  geom_point(aes(y = `0.5`, color = forecast_date), size = 0.75) +
  geom_line(
    data = params$truth_data %>%
      filter(
        geo_value %in% geo_vals,
        target_end_date >= min(plot_dates),
        target_end_date <= max(plot_dates) + 5 * 7
      ),
    aes(x = target_end_date, y = true_value),
    inherit.aes = FALSE, na.rm = TRUE,
    color = "black", linetype = 1
  ) +
  scale_y_continuous(n.breaks = 15) +
  facet_grid(factor(forecaster, levels = param_table$id) ~ geo_value, scales = "free") +
  labs(x = "Reference Date", y = "Weekly Sums of Hospitalizations", title = "Monthly Forecasts and Truth Data") +
  theme(legend.position = "none")

ggplotly(p, tooltip = "text", height = 300 * length(param_table$id), width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))