library(aws.s3)
library(data.table)
library(dplyr)
library(DT)
library(ggplot2)
library(plotly)
library(readr)
library(stringr)
library(tidyr)
# params <- list(forecaster_set = "scaled_pop_season")

# Load the table of parameter combinations
s3load("flu_2023_forecaster_parameter_combinations.rds", bucket = "forecasting-team-data")

# Select forecasters for this notebook
cmu_forecasters <- forecaster_parameter_combinations_[[params$forecaster_set]]$id
base_forecaster_name <- "FluSight-baseline"
ensemble_forecaster_name <- "FluSight-ensemble"
outside_forecasters <- c(base_forecaster_name, ensemble_forecaster_name)

# Load scores and filter them, get global variables
s3load(object = "flu_2023_joined_scores.rds", bucket = "forecasting-team-data")
scores <- flu_2023_joined_scores %>%
    filter(forecaster %in% c(cmu_forecasters, outside_forecasters))
forecast_dates <- scores %>%
    pull(forecast_date) %>%
    unique()
forecasters <- c(cmu_forecasters, outside_forecasters)
aheads <- scores %>%
    pull(ahead) %>%
    unique()
forecaster_family <- params$forecaster_set

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

Notebook Information

Forecaster Family

no_recent_quant

Target Dates

2023-10-14, 2023-10-21, 2023-10-28, 2023-11-04, 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, 2023-10-07

Forecaster Parameter Mapping and Overall Scores

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

# Display the table
param_table <- forecaster_parameter_combinations_[[params$forecaster_set]] %>%
    select(-any_of(c("forecaster", "keys_to_ignore", "pop_scaling"))) %>%
    {
        if ("n_training" %in% colnames(.)) {
            (.) %>% mutate(n_training = as.character(n_training))
        } else {
            .
        }
    } %>%
    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_80 = round(Mean(coverage_80), 2),
            ) %>%
            rename(id = forecaster)
    ) %>%
    arrange(mean_ae) %>%
    relocate(id, mean_ae, geomean_ae, mean_wis, geomean_wis, mean_coverage_80)
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 %>% 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 = "Geometric Mean WIS")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
    layout(hoverlabel = list(bgcolor = "white"))
3e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e-011e+003e+001e+013e+011e+023e+021e+033e+031e+04OctJanApr3e-011e+003e+001e+013e+011e+023e+021e+033e+031e+04
forecastercrazed.whitepelicanfungal.bantengsoutheastern.amphibianwised.terrapinForecasts made over October 07, 2023 to April 20, 2024Forecast DatesGeometric Mean WIS0 days ahead7 days ahead14 days ahead21 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"))
3e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e+041e+053e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e+041e+053e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e+041e+05OctJanApr3e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e+041e+05
forecastercrazed.whitepelicanfungal.bantengsoutheastern.amphibianwised.terrapinForecasts made over October 07, 2023 to April 20, 2024Forecast DatesArithmetic Mean WIS0 days ahead7 days ahead14 days ahead21 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"))
0510152013103010030010003000
forecastercrazed.whitepelicanfungal.bantengsoutheastern.amphibianwised.terrapinForecasts made over October 07, 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"))
051015205e-011e+002e+003e+005e+001e+012e+013e+015e+011e+022e+023e+025e+021e+032e+033e+035e+031e+042e+043e+04
forecastercrazed.whitepelicanfungal.bantengsoutheastern.amphibianwised.terrapinForecasts made over October 07, 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"))
3e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e+041e+053e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e+041e+053e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e+041e+05OctJanApr3e-011e+003e+001e+013e+011e+023e+021e+033e+031e+043e+041e+05
forecastercrazed.whitepelicanfungal.bantengsoutheastern.amphibianwised.terrapinForecasts made over October 07, 2023 to April 20, 2024Forecast DatesArithmetic Mean Absolute Error0 days ahead7 days ahead14 days ahead21 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"))
051015205e-017e-011e+002e+003e+005e+007e+001e+012e+013e+015e+017e+011e+022e+023e+025e+027e+021e+032e+033e+035e+037e+031e+042e+043e+045e+04
forecastercrazed.whitepelicanfungal.bantengsoutheastern.amphibianwised.terrapinForecasts made over October 07, 2023 to April 20, 2024Days aheadArithmetic Mean Absolute Error
var <- "coverage_80"
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 80% Coverage")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
    layout(hoverlabel = list(bgcolor = "white"))
0.000.250.500.751.000.000.250.500.751.000.000.250.500.751.00OctJanApr0.000.250.500.751.00
forecastercrazed.whitepelicanFluSight-baselinefungal.bantengsoutheastern.amphibianwised.terrapinForecasts made over October 07, 2023 to April 20, 2024Forecast DatesArithmetic Mean 80% Coverage0 days ahead7 days ahead14 days ahead21 days ahead
var <- "coverage_80"
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 80% Coverage")

ggplotly(p, tooltip = "text", height = 800, width = 1000) %>%
    layout(hoverlabel = list(bgcolor = "white"))
051015200.00.20.40.60.8
forecastercrazed.whitepelicanFluSight-baselinefungal.bantengsoutheastern.amphibianwised.terrapinForecasts made over October 07, 2023 to April 20, 2024Days aheadArithmetic Mean 80% Coverage

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.

# Load the forecasts and the truth data
s3load(object = "flu_2023_joined_forecasts.rds", bucket = "forecasting-team-data")
s3load(object = "flu_2023_truth_data.rds", bucket = "forecasting-team-data")

# We plot a subset of the dates and geos for the fan plot
plot_dates <- seq.Date(as.Date("2023-10-07"), by = "4 weeks", length.out = 8)
geo_vals <- c("ca", "fl", "pa", "tx")
forecast_subset <- flu_2023_joined_forecasts %>%
    filter(
        forecaster %in% c(cmu_forecasters, outside_forecasters),
        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 = flu_2023_truth_data %>% filter(geo_value %in% geo_vals),
        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_y") +
    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"))
02505007501000125015001750200022502500275030003250350002004006008001000120014001600180020002200050000010000001500000200000025000003000000350000040000000500000100000015000002000000250000030000003500000400000005000001000000150000020000002500000300000035000004000000OctJanApr05000001000000150000020000002500000300000035000004000000OctJanAprOctJanAprOctJanApr
Monthly Forecasts and Truth DataReference DateWeekly Sums of HospitalizationscaflpatxFluSight-ensembleFluSight-baselinesoutheastern.amphibianwised.terrapincrazed.whitepelicanfungal.banteng