## ----include=FALSE------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4,
  warning = FALSE,
  message = FALSE
)
oldpar <- par(no.readonly = TRUE)
oldopts <- options()
library(dplyr)
library(ggplot2)
library(kableExtra)

## ----setup--------------------------------------------------------------------
library(modelimportance)
library(dplyr)
library(ggplot2)

## ----forecast_data------------------------------------------------------------
# Specify forecasts to remove: MOBS-GLEAM_FLUH for location 25 on 2022-11-26,
# PSI-DICE for location 48 on 2022-12-10
forecast_to_remove <- tibble(
  model_id = c("MOBS-GLEAM_FLUH", "PSI-DICE"),
  location = c("25", "48"),
  target_end_date = as.Date(c("2022-11-26", "2022-12-10"))
)

# Filter out the specified forecasts from the original data
forecast_data <- forecast_data_raw |>
  anti_join(forecast_to_remove, by = c("model_id", "location", "target_end_date"))

# Display the forecast data
forecast_data

## ----target_data--------------------------------------------------------------
target_data <- target_data_raw |>
  dplyr::filter(
    target_end_date %in% unique(forecast_data$target_end_date),
    location %in% unique(forecast_data$location),
    target == "wk inc flu hosp"
  )

target_data

## ----fig-example-median-lomo, echo=FALSE--------------------------------------
target_data |>
  mutate(model_id = "Observed") |>
  rename(value = oracle_value) |>
  rbind(
    forecast_data |>
      select("target_end_date", "target", "location", "value", "model_id")
  ) |>
  ggplot(aes(x = target_end_date)) +
  geom_point(aes(y = value, color = model_id, shape = model_id), size = 2) +
  geom_point(
    data = target_data,
    aes(y = oracle_value, color = "Observed", shape = "Observed"),
    size = 3
  ) +
  facet_wrap(~location,
    scales = "free_y",
    labeller = labeller(location = function(x) paste0("Location: ", x))
  ) +
  scale_x_date(
    breaks = target_data$target_end_date,
    date_labels = "%Y-%m-%d",
    expand = expansion(add = c(5, 5))
  ) +
  scale_color_manual(
    name = "model_id",
    values = c(
      "Flusight-baseline" = "#F8766D",
      "MOBS-GLEAM_FLUH" = "#00BA38",
      "PSI-DICE" = "#619CFF",
      "Observed" = "black"
    ),
    limits = c(
      "Flusight-baseline",
      "MOBS-GLEAM_FLUH",
      "PSI-DICE",
      "Observed"
    )
  ) +
  scale_shape_manual(
    name = "model_id",
    values = c(
      "Flusight-baseline" = 16,
      "MOBS-GLEAM_FLUH" = 17,
      "PSI-DICE" = 15,
      "Observed" = 18
    ),
    limits = c(
      "Flusight-baseline",
      "MOBS-GLEAM_FLUH",
      "PSI-DICE",
      "Observed"
    )
  ) +
  labs(
    y = "Weekly Hospitalization",
    x = "Date",
    title = "Forecasts of incident deaths generated on November 19, 2022"
  )

## ----lomo---------------------------------------------------------------------
scores_lomo <- model_importance(
  forecast_data = forecast_data,
  oracle_output_data = target_data,
  ensemble_fun = "simple_ensemble",
  importance_algorithm = "lomo"
)

print(scores_lomo)

## -----------------------------------------------------------------------------
summary(scores_lomo)

## -----------------------------------------------------------------------------
s <- summary(scores_lomo)
s$all_tasks

## -----------------------------------------------------------------------------
s$model_summary

## -----------------------------------------------------------------------------
s$task_winners

## ----echo=FALSE---------------------------------------------------------------
ggplot(scores_lomo, aes(x = model_id, y = importance, fill = model_id)) +
  geom_col() +
  coord_flip() +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.25) +
  facet_grid(
    cols = vars(target, horizon, location, target_end_date),
    scales = "free_x"
  ) +
  labs(
    x = "Model ID", y = "Importance Score",
    title = "Model Importance by Task"
  )

## -----------------------------------------------------------------------------
aggregate(scores_lomo, by = "model_id", na_action = "drop", fun = mean)

## -----------------------------------------------------------------------------
aggregate(scores_lomo, by = "model_id", na_action = "worst", fun = mean)

## -----------------------------------------------------------------------------
aggregate(scores_lomo, by = "model_id", na_action = "average", fun = mean)

## ----code-model_importance-summary-lasomo-equal-------------------------------
scores_lasomo_eq <- model_importance(
  forecast_data = forecast_data,
  oracle_output_data = target_data,
  ensemble_fun = "simple_ensemble",
  importance_algorithm = "lasomo",
  subset_wt = "equal"
)
aggregate(scores_lasomo_eq, by = "model_id", na_action = "drop", fun = mean)

## ----code-model_importance-summary-lasomo-perm--------------------------------
scores_lasomo_perm <- model_importance(
  forecast_data = forecast_data,
  oracle_output_data = target_data,
  ensemble_fun = "simple_ensemble",
  importance_algorithm = "lasomo",
  subset_wt = "perm_based"
)
aggregate(scores_lasomo_perm, by = "model_id", na_action = "drop", fun = mean)

## ----include=FALSE------------------------------------------------------------
options(oldopts)
par(oldpar)