## ----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(tidyr)
library(ggplot2)
library(kableExtra)
library(modelimportance)

## ----results='asis'-----------------------------------------------------------
forecast_data_example |>
  dplyr::filter(
    .data$reference_date == "2022-12-17",
    .data$horizon != 0
  ) |>
  head(10) |>
  knitr::kable(
               format = "html",
               caption = "Table 1: Example of the model output for incident
               influenza hospitalizations (top 10 rows) extracted from
               `forecast_data_example` bundled in the `modelimportance` package, 
               which is originally from `hubExamples` package.") |>
  kableExtra::kable_styling(
    font_size = 12,
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

## ----fig.pos="t"--------------------------------------------------------------
forecast_data_example |>
  filter(
    reference_date == "2022-12-17",
    horizon != 0,
    output_type_id %in% c(0.05, 0.5, 0.95)
  ) |>
  pivot_wider(names_from = output_type_id, values_from = value) |>
  rename(lower = "0.05", upper = "0.95", value = "0.5") |>
  ggplot(aes(x = target_end_date)) +
  facet_grid(~model_id) +
  geom_point(aes(y = value, color = "medians"), size = 2) +
  geom_line(aes(y = value, color = "medians"), linewidth = 1) +
  geom_ribbon(
    aes(
      ymin = lower, ymax = upper,
      fill = "#3388FF"
    ),
    alpha = 0.5
  ) +
  geom_point(
    data = target_data_example |> filter(target_end_date <= "2022-12-17"),
    aes(y = observation, group = 1, color = "obs")
  ) +
  geom_line(
    data = target_data_example |> filter(target_end_date <= "2022-12-17"),
    aes(y = observation, group = 1, color = "obs")
  ) +
  geom_point(
    data = target_data_example |> filter(target_end_date > "2022-12-17"),
    aes(y = observation, group = 1, color = "truth"),
    shape = 1, alpha = 1
  ) +
  geom_line(
    data = target_data_example |> filter(target_end_date > "2022-12-17"),
    aes(y = observation, group = 1, color = "truth"),
    alpha = 0.75
  ) +
  # coord_cartesian(ylim = c(0, 500)) +
  scale_x_date(breaks = target_data_example$target_end_date, date_labels = "%Y-%m-%d") +
  labs(
    y = "Weekly Hospitalization",
    x = "Date"
  ) +
  scale_color_manual(
    name = "",
    values = c(
      "medians" = "DodgerBlue",
      "obs" = "Black",
      "truth" = "Black"
    ),
    labels = c(
      "Forecast (Predictive median)",
      "Observed data before forecasting",
      "Eventually observed value"
    )
  ) +
  scale_fill_manual("",
    values = "#3388FF",
    labels = "90% Prediction interval"
  ) +
  theme(
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8),
    axis.text.x = element_text(size = 6, angle = 90),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 8),
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    legend.position = "bottom",
    legend.direction = "vertical",
    legend.box = "horizontal"
  )

## -----------------------------------------------------------------------------
data.frame(
  "Output Type" = c("mean", "median", "quantile", "pmf"),
  "Scoring Rule" = c("SE", "AE", "WIS", "Log Score"),
  Description = c(
    "Squared error (SE): the squared difference between the predicted value and the observed value",
    "Absolute error (AE): the absolute difference between the predicted value and the observed value",
    paste0(
      "Weighted interval score (WIS): a quantile-based approximation of the continuous ranked",
      " probability score (CRPS) for evaluating quantile forecasts"
    ),
    "Logarithm of the probability assigned to the true outcome (LogScore)"
  ),
  check.names = FALSE
) |>
  knitr::kable(
    format = "html",
    caption = paste0(
      "Table 2: Pairs of output types and their associated scoring rules ",
      "for evaluating prediction performance."
    )
  ) |>
  kableExtra::kable_styling(
    font_size = 13,
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

## -----------------------------------------------------------------------------
knitr::include_graphics("figure-concept.png")

## -----------------------------------------------------------------------------
knitr::include_graphics("algorithm-lomo.jpg")

## -----------------------------------------------------------------------------
knitr::include_graphics("algorithm-lasomo.jpg")

## ----fig.pos="t"--------------------------------------------------------------
data.frame(n = 2:10) |>
  mutate(
    w_eq = 1 / (2^(n - 1) - 1),
    w_perm_min = 1 / ((n - 1) * choose(n - 1, floor((n - 1) / 2))),
    w_perm_max = 1 / (n - 1)
  ) |>
  ggplot(aes(x = n)) +
  geom_line(aes(y = w_eq, color = "eq"), size = 0.75, linetype = "longdash") +
  geom_point(aes(y = w_eq, color = "eq"), size = 2, shape = 15) +
  geom_line(aes(y = w_perm_min, color = "perm_min"), size = 0.75) +
  geom_point(aes(y = w_perm_min, color = "perm_min"), size = 2, shape = 16) +
  geom_line(aes(y = w_perm_max, color = "perm_max"), size = 0.75) +
  geom_point(aes(y = w_perm_max, color = "perm_max"), size = 2, shape = 17) +
  scale_color_manual(
    name = "Weights",
    values = c(
      "eq" = "#F8766D",
      "perm_min" = "#619CFF",
      "perm_max" = "#00BA38"
    ),
    labels = c(
      eq = expression(w^{
        eq
      }),
      perm_min = expression(w^{
        perm - min
      }),
      perm_max = expression(w^{
        perm - max
      })
    )
  ) +
  scale_x_continuous(breaks = 2:10) +
  labs(
    x = "Number of models (n)",
    y = "Weight assigned to a subset",
    color = "Weighting Scheme"
  ) +
  theme(
    axis.title.x = element_text(size = 9),
    axis.title.y = element_text(size = 9),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    legend.title = element_text(size = 9),
    legend.position = "bottom",
    legend.text = element_text(size = 9),
    legend.spacing.x = unit(0.5, "mm"),
    legend.key.width = unit(1, "cm")
  )

## ----code-model_importance----------------------------------------------------
# model_importance(
#   forecast_data, oracle_output_data, ensemble_fun,
#   importance_algorithm, subset_wt, min_log_score,
#   ...
# )

## -----------------------------------------------------------------------------
tbl_argument1 <- data.frame(
  "Argument" = c(
    "`forecast_data`", "`oracle_output_data`",
    "`ensemble_fun`",
    "`importance_algorithm`", "`subset_wt`",
    "`min_log_score`", "`...`"
  ),
  "Description" = c(
    "Forecasts",
    "Ground truth data",
    "Ensemble method",
    "Algorithm to calculate importance",
    "Method for assigning weight to subsets when using LASOMO algorithm",
    "Minimum value to replace for log score",
    "Optional arguments for `'simple_ensemble'`"
  ),
  "Possible Values" = c(
    "Table of model output",
    "Table of oracle output",
    "`'simple_ensemble'`, `'linear_pool'`",
    "`'lomo', 'lasomo'`",
    "`'equal', 'perm_based'`",
    "Non-positive numeric",
    "Varies"
  ),
  Default = c(
    "N/A",
    "N/A",
    "`'simple_ensemble'`",
    "`'lomo'`",
    "`'equal'`",
    -10,
    "`agg_fun='mean'`"
  ),
  check.names = FALSE
)


tbl_argument1 |>
  knitr::kable(format = "html", escape = FALSE,
               caption = "Table 3: Description of the arguments for the
               `model_importance()` function, including their purpose, possible
               values, and default settings.") |>
  kableExtra::kable_styling(
    font_size = 13,
    bootstrap_options = c("striped", "condensed", "responsive"),
    full_width = FALSE
  )

## -----------------------------------------------------------------------------
tbl_argument2 <- data.frame(
  "Argument" = c(
    "`importance_scores`",
    "`by`",
    "`na_action`",
    "`fun`", "`...`"
  ),
  "Description" = c(
    "Model importance scores produced by `model_importance()`",
    "Grouping variable(s) for summarization",
    "Method to handle `NA` values",
    "Function to summarize importance scores",
    "Optional arguments for `\"fun\"`"
  ),
  "Possible Values" = c(
    "data frame",
    "grouping variable(s)",
    "`\"drop\", \"worst\", \"average\"`",
    "summary function",
    "depends on `fun`"
  ),
  Default = c(
    "N/A",
    "`\"model_id\"`",
    "`\"drop\"`",
    "`mean`",
    "N/A"
  ),
  check.names = FALSE
)

tbl_argument2 |>
  knitr::kable(
    format = "html", escape = FALSE,
    caption = paste0(
      "Table 4: Description of the arguments for the `aggregate()` function for ",
      "`model_imp_tbl` objects, including their purpose, possible values, and default settings."
    )
  ) |>
  kableExtra::kable_styling(
    font_size = 13,
    bootstrap_options = c("striped", "condensed", "responsive"),
    full_width = FALSE
  )

## ----code-forecast_data-------------------------------------------------------
# forecast_data

## ----forecast_data, tidy = FALSE----------------------------------------------
forecast_data <- readRDS("vignette-example-forecast_data.rds")
forecast_data |>
  knitr::kable(
    format = "html"
  ) |>
  kableExtra::kable_styling(
    font_size = 12,
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

## ----code-target_data---------------------------------------------------------
# target_data

## ----target_data--------------------------------------------------------------
target_data <- readRDS("vignette-example-oracle_data.rds")
target_data |>
  knitr::kable(
    format = "html"
  ) |>
  kableExtra::kable_styling(
    font_size = 12,
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

## ----fig.pos="t"--------------------------------------------------------------

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"
  ) +
  theme(
    axis.title.x = element_text(size = 9),
    axis.title.y = element_text(size = 9),
    axis.text.x = element_text(size = 9),
    axis.text.y = element_text(size = 9),
    strip.text = element_text(size = 9),
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.spacing.x = unit(0.25, "mm"),
    legend.key.width = unit(0.4, "cm")
  )

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

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

## -----------------------------------------------------------------------------
# print(
#   scores_lomo |>
#     mutate(importance = round(importance, 2)) |>
#     rename(
#       ref_date = reference_date, h = horizon,
#       loc = location, t_end_date = target_end_date,
#       o_type = output_type, imp = importance
#     )
# )

## -----------------------------------------------------------------------------
print(scores_lomo |>
  mutate(importance = round(importance, 2)) |>
  rename(
    ref_date = reference_date, h = horizon, loc = location,
    t_end_date = target_end_date, o_type = output_type, imp = importance
  ))

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

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

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

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

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

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

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

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

## -----------------------------------------------------------------------------
# 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"
#   ) +
#   scale_x_discrete(labels = function(x) gsub("[-_]", "-\n", x)) +
#   theme(
#     axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
#     panel.spacing.x = unit(0.5, "lines"), legend.position = "none"
#   )

## -----------------------------------------------------------------------------
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"
  ) +
  scale_x_discrete(labels = function(x) gsub("[-_]", "-\n", x)) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    panel.spacing.x = unit(0.5, "lines"),
    legend.position = "none"
  )

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

## -----------------------------------------------------------------------------
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 = "worst", fun = mean)

## -----------------------------------------------------------------------------
# aggregate(scores_lomo, by = "model_id", na_action = "average", 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)

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

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

## -----------------------------------------------------------------------------
knitr::include_graphics("runtime_comparison-fitted.png")

## -----------------------------------------------------------------------------
knitr::include_graphics("runtime_parallel-fitted.png")

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