In epiforecasts/europe-covid-forecast: Forecasting Covid-19 related metrics for the ECDC Forecast Hub

library(scoringutils)
library(ggplot2)
library(dplyr)
library(DT)
library(data.table)
library(here)
library(knitr)
library(covidHubUtils)
library(lubridate)
library(purrr)
options(knitr.duplicate.label = "allow")
knitr::opts_chunk$set(echo = FALSE, 
                      message = FALSE, 
                      warning = FALSE)

include_ranking <- TRUE
include_forecast_plot <- TRUE
include_avail_forecasts <- FALSE
include_ranks_over_time <- FALSE

include_countries <- include_forecast_plot || 
  include_forecast_plot || 
  include_avail_forecasts || 
  include_ranks_over_time

data <- params$data

horizons <- data %>% 
  filter(!is.na(true_value), !is.na(horizon)) %>% 
  .$horizon %>% 
  unique

horizons <- horizons[as.integer(horizons) < 5]
horizons <- horizons[order(as.integer(horizons))]

summarise_by = c("model", "target_variable")

target_variables <- c(Cases = "inc case", Deaths = "inc death")

{.unlisted .unnumbered}

---

out <- paste0("\n\n## Forecast visualisation (", params$location_name,
              ") {.tabset .tabset_fade}\n\n")

out <- c(out, "Forecast visualisations. The date of the tab marks the date on",
         "which a forecast was made.\n\n")

r paste(if (include_forecast_plot) knit(text = out), collapse = '\n\n')

out <- NULL
forecast_dates <-
  data$forecast_date[!is.na(data$forecast_date)] %>%
  unique() %>%
  sort() %>%
  as.character() %>%
  rev()
loc <- params$location_name
out <- c(out, knit_child(here::here("reports", "evaluation",
                                    "template-plot-forecasts.Rmd")))

r paste(if (include_forecast_plot) knit(text = out), collapse = '\n\n')

{.unlisted .unnumbered}

---

out <- paste0("\n\n## Forecast scores (", loc, ") {.tabset .tabset-fade}\n\n")

out <- c(out, "Scores separated by target and forecast horizon.\n\n")

r paste(if (include_ranking) knit(text = out, quiet = TRUE), collapse = '\n\n')

out <- NULL
for (variable in names(target_variables)) {
  out <- c(out, paste("\n\n###", variable, " {.tabset .tabset-fade}\n\n"))
  for (this_horizon in horizons) {
    filter_list = list(paste0("target_variable == '",
                              target_variables[[variable]], "'"),
                     "!(is.na(quantile))",
                     paste0("horizon == ", this_horizon),
                     paste0("location_name == '", loc, "'"))
    horizon_string <-
      paste0(this_horizon, " week",
             if_else(this_horizon > 1, "s", ""), " ahead horizon")
    out <- c(out, paste("\n\n####", horizon_string, "\n\n"))
    out <- c(out, knit_child(here::here("reports", "evaluation",
                                        "template-ranking-table.Rmd"),
                            quiet = TRUE))
  }
}

r paste(if (include_ranking) knit(text = out, quiet = TRUE), collapse = '\n\n')

{.unlisted .unnumbered}

---

Evaluation metrics

• Relative skill is a metric based on the weighted interval score (WIS) that is using a 'pairwise comparison tournament'. All pairs of forecasters are compared against each other in terms of the weighted interval score. The mean score of both models based on the set of common targets for which both models have made a prediction are calculated to obtain mean score ratios. The relative skill is the geometric mean of these mean score ratios. Smaller values are better and a value smaller than one means that the model beats the average forecasting model.
• The weighted interval score is a proper scoring rule (meaning you can't cheat it) suited to scoring forecasts in an interval format. It has three components: sharpness, underprediction and overprediction. Sharpness is the width of the prediction interval. Over- and underprediction only come into play if the prediction interval does not cover the true value. They are the absolute value of the difference between the upper or lower bound of the prediction interval (depending on whether the forecast is too high or too low).
• coverage deviation is the average difference between nominal and empirical interval coverage. Say the 50 percent prediction interval covers only 20 percent of all true values, then the coverage deviation is 0.5 - 0.2 = -0.3. The coverage deviation value in the table is calculated by averaging over the coverage deviation calculated for all possible prediction intervals. If the value is negative you have covered less then you should. If it is positve, then the forecasts could be a little more confident.
• bias is a measure between -1 and 1 that expresses the tendency to underpredict (-1) or overpredict (1). In contrast to the over- and underprediction components of the WIS it is bound between -1 and 1 and cannot go to infinity. It is therefore less susceptible to outliers.
• aem is the absolute error of the median forecasts. A high aem means the median forecasts tend to be far away from the true values.

{.unlisted .unnumbered}

---

out <- paste0("\n\n## Scores over time (", loc, ") {.tabset .tabset-fade}\n\n")
out <- c(out, "Visualisation of the weighted interval score over time. In",
         "addition, the components of the interval score, sharpness (how",
         "narrow are forecasts - smaller is better), and penalties for",
         "underprediction and overprediction are shown. Scores are again",
         "separated by forecast horizon\n\n")

r paste(knit(text = out), collapse = '\n\n')

out <- NULL

for (this_horizon in horizons) {
  horizon_string <-
    paste0(this_horizon, " week", if_else(this_horizon > 1, "s", ""),
           " ahead horizon {.tabset .tabset-fade}\n\n ")
    out <- c(out, paste("\n\n###", horizon_string, "\n\n"))
    out <- c(out, knit_child(here::here("reports", "evaluation",
                                        "template-scores-and-truth-time.Rmd"),
                             quiet = TRUE))
  }

r paste(knit(text = out), collapse = '\n\n')

If you want to learn more about a model, you can go the the 'data-processed'-folder of the European Forecast Hub github repository, select a model and access the metadata file with further information provided by the model authors.

epiforecasts/europe-covid-forecast documentation built on Jan. 15, 2025, 8:57 p.m.