# outside_forecaster_subset <- c("COVIDhub-baseline", "COVIDhub-trained_ensemble", "COVIDhub_CDC-ensemble")
# i <- 5
# forecaster_families <- names(tar_read(forecaster_parameter_combinations))
# forecaster_family <- forecaster_families[[i]]
# params_subset <- tar_read(forecaster_parameter_combinations)[[forecaster_family]]
# forecasts_subset <- tar_read(joined_forecasts) %>% filter(forecaster %in% c(params_subset$id, outside_forecaster_subset))
# scores_subset <- tar_read(joined_scores) %>% filter(forecaster %in% c(params_subset$id, outside_forecaster_subset))

# params <- list(
#   forecaster_family = forecaster_family,
#   forecaster_parameters = params_subset,
#   forecasts = forecasts_subset,
#   scores = scores_subset,
#   truth_data = tar_read(hhs_evaluation_data),
#   disease = "covid"
# )

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"
}

# Load scores and filter them, get global variables
scores <- params$scores
forecast_dates <- scores %>%
  pull(forecast_date) %>%
  unique()
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

Disease

flu

Forecaster Family

flatline

Target 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

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", "outcome")
}

param_table <- params$forecaster_parameters %>%
  select(-any_of(ignore_keys)) %>%
  mutate(
    across(matches("^n_training$"), ~ as.character(.x)),
    across(matches("^trainer$"), ~ unlist(.x)),
    across(matches("^seasonal_method$"), ~ Vectorize(function(x) paste0(x, collapse = ", "))(.x))
  ) %>%
  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(id, mean_ae, geomean_ae, mean_wis, geomean_wis, mean_coverage_90)
datatable(param_table)

Score Plots

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"))
0.40.60.81.01.21.41.61.80.40.60.81.01.21.41.61.8DecJanFebMarApr0.40.60.81.01.21.41.61.8DecJanFebMarApr
forecasterblitheful.mayflyFluSight-ensembleUMass-flusionForecasts made over November 11, 2023 to April 20, 2024Forecast DatesGeometric Mean WIS0 days ahead7 days ahead14 days ahead21 days ahead28 days ahead
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 := 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 %>% filter(forecaster != ensemble_forecaster_name),
  aes(x = forecast_date, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_line(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = forecast_date, y = !!sym(var)),
    color = "black", linetype = 2
  ) +
  geom_point(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = forecast_date, y = !!sym(var)),
    color = "black", shape = 21, fill = "white"
  ) +
  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 WIS")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))
0.30.40.50.60.70.80.91.00.30.40.50.60.70.80.91.0DecJanFebMarApr0.30.40.50.60.70.80.91.0DecJanFebMarApr
forecasterblitheful.mayflyUMass-flusionForecasts made over November 11, 2023 to April 20, 2024Forecast DatesArithmetic Mean WIS0 days ahead7 days ahead14 days ahead21 days ahead28 days ahead
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 %>% filter(forecaster != ensemble_forecaster_name),
  aes(x = ahead, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_line(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = ahead, y = !!sym(var)),
    color = "black", linetype = 2
  ) +
  geom_point(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = ahead, y = !!sym(var)),
    color = "black", shape = 21, fill = "white"
  ) +
  geom_hline(yintercept = 1, linetype = 1, color = "black") +
  scale_y_log10(n.breaks = 10) +
  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"))
010200.650.700.750.800.850.900.951.001.051.10
forecasterblitheful.mayflyUMass-flusionForecasts made over November 11, 2023 to April 20, 2024Days aheadGeometric Mean WIS
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 := 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 %>% filter(forecaster != ensemble_forecaster_name),
  aes(x = ahead, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_line(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = ahead, y = !!sym(var)),
    color = "black", linetype = 2
  ) +
  geom_point(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = ahead, y = !!sym(var)),
    color = "black", shape = 21, fill = "white"
  ) +
  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 = "Arithmetic Mean WIS")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
  layout(hoverlabel = list(bgcolor = "white"))
010200.540.570.600.630.660.690.720.750.780.810.840.870.900.930.960.991.021.051.081.11
forecasterblitheful.mayflyUMass-flusionForecasts made over November 11, 2023 to April 20, 2024Days aheadArithmetic Mean WIS
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 %>% filter(forecaster != ensemble_forecaster_name),
  aes(x = forecast_date, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_line(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = forecast_date, y = !!sym(var)),
    color = "black", linetype = 2
  ) +
  geom_point(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = forecast_date, y = !!sym(var)),
    color = "black", shape = 21, fill = "white"
  ) +
  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"))
0.30.40.50.60.70.81.02.00.30.40.50.60.70.81.02.0DecJanFebMarApr0.30.40.50.60.70.81.02.0DecJanFebMarApr
forecasterblitheful.mayflyUMass-flusionForecasts made over November 11, 2023 to April 20, 2024Forecast DatesArithmetic Mean Absolute Error0 days ahead7 days ahead14 days ahead21 days ahead28 days ahead
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 %>% filter(forecaster != ensemble_forecaster_name),
  aes(x = ahead, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_line(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = ahead, y = !!sym(var)),
    color = "black", linetype = 2
  ) +
  geom_point(
    data = normalized_df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = ahead, y = !!sym(var)),
    color = "black", shape = 21, fill = "white"
  ) +
  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"))
010200.600.620.640.660.680.700.720.740.760.780.800.820.840.860.880.900.920.940.960.981.001.021.04
forecasterblitheful.mayflyUMass-flusionForecasts made over November 11, 2023 to April 20, 2024Days aheadArithmetic Mean Absolute Error
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 %>% filter(forecaster != ensemble_forecaster_name),
  aes(x = forecast_date, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_line(
    data = df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = forecast_date, y = !!sym(var)),
    color = "black", linetype = 2
  ) +
  geom_point(
    data = df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = forecast_date, y = !!sym(var)),
    color = "black", shape = 21, fill = "white"
  ) +
  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"))
0.000.250.500.751.000.000.250.500.751.00DecJanFebMarApr0.000.250.500.751.00DecJanFebMarApr
forecasterblitheful.mayflyFluSight-baselineUMass-flusionForecasts made over November 11, 2023 to April 20, 2024Forecast DatesArithmetic Mean 90% Coverage0 days ahead7 days ahead14 days ahead21 days ahead28 days ahead
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 %>% filter(forecaster != ensemble_forecaster_name),
  aes(x = ahead, y = !!sym(var))
) +
  geom_line(aes(color = forecaster, group = forecaster)) +
  geom_point(aes(color = forecaster, group = forecaster)) +
  geom_line(
    data = df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = ahead, y = !!sym(var)),
    color = "black", linetype = 2
  ) +
  geom_point(
    data = df %>% filter(forecaster == ensemble_forecaster_name),
    aes(x = ahead, y = !!sym(var)),
    color = "black", shape = 21, fill = "white"
  ) +
  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"))
010200.60.70.80.9
forecasterblitheful.mayflyFluSight-baselineUMass-flusionForecasts made over November 11, 2023 to April 20, 2024Days aheadArithmetic Mean 90% Coverage

Fan plots

Fan plots showing the 90% 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(
    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))

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

# 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(
    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))

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

# 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(
    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))

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