Local hub playground

Introduction

In this session we will try to take some of the lessons on building and evaluating forecasts in practice from the earlier sessions and and apply them. Parts of this session will be more open-ended than other sessions, enabling you to spend time doing some hands-on modeling of your own, or spending time reviewing some of the best practices and key concepts that we have covered in the course so far.

Slides

Why use a local hub?

Objectives

The aims of this session are

to introduce the use of modeling hubs and hubverse-style tools for local model development and collaborative modeling projects, and
to practice building forecasting models using real data.

Setup

Source file

The source file of this session is located at sessions/hub-playground.qmd.

Libraries used

In this session we will use the dplyr package for data wrangling, the ggplot2 and ggtime packages for plotting, the fable package for building forecasting models, and the epidatr package for downloading epidemiological surveillance data.

Additionally, we will use some hubverse packages such as hubEvals, hubUtils, hubData, hubVis, hubEnsembles packages for building ensembles.

library("nfidd.forecasting")
library("dplyr")
library("ggplot2")
library("epidatr")
library("fable")
library("ggtime")
library("hubUtils")
library("hubEvals")
library("hubVis")
library("hubData")
library("hubEnsembles")
theme_set(theme_bw())

Tip

The best way to interact with the material is via the Visual Editor of RStudio.

Initialisation

We set a random seed for reproducibility. Setting this ensures that you should get exactly the same results on your computer as we do. This is not strictly necessary but will help us talk about the models.

set.seed(406) # for Ted Williams

A flu forecasting hub for model dev

We have built a live modeling hub for all participants of the class to use as a “sandbox” for building models and running analyses.

Tip

We are referring to this hub as a “sandbox” hub because it’s not an active hub in use by decision-makers. It’s a sandbox where we can build models and try things out without worrying about breaking them.

We will also use the term local hub in this session. By this we mean any modeling hub that you set up on your local filesystem to do a modeling project. We’re going to set up a hub in the style of the hubverse, because you get a lot of added benefits from working with data in that format. But you could build a local hub however you want for your own modeling experiments.

If you complete this session and the next one, you will be able to see your forecasts alongside the other forecasts at the online dashboard for this hub.

You should download this hub to your local laptop before continuing this session. Here are the two ways you can download this hub:

Download the hub as a fixed dataset using the link to this zip file.
If you are familiar with using GitHub (or want to get more familiar with it), then we recommend cloning the repo directly to your laptop.
- Login to GitHub using your existing GitHub account (or create one!).
- Clone this repository by clicking on the green “<> Code” button on the main page of the Sandbox Hub repo.
- For the rest of the code in this session to work smoothly, this repository should live inside the sismid/ directory that stores the material for this course.

Tip

Whether you use method 1 or 2 listed above, you should end up with a directory called sismid/sismid-ili-forecasting-sandbox/ on your local machine.

Caution

Note that you can clone the repo using HTTPS, but if you are going to be using and pushing to GitHub a lot, then setting up SSH keys and cloning via SSH is usually worth it.

ILI Forecasting at HHS Region level

We have done some modeling on ILI data in the earlier sessions. Our sandbox hub is set up to accept forecasts for the US level and the 10 Health and Human Services (HHS) Regions.

A forecast will produce quantile predictions of 1- through 4-week ahead forecasts of the estimated percentage of outpatient doctors office visits that are for ILI. This target is given the label “ili perc” in the hub. Forecasts can be submitted for any week during the 2015/2016, 2016/2017, 2017/2018, 2018/2019 or 2019/2020 influenza seasons.

Details of this hub

There are six models that already have made forecasts that live in the hub: hist-avg, delphi-epicast, protea-cheetah, lanl-dbmplus, kot-kot, and neu-gleam.

We will assume that you have the hub repository downloaded to your local machine as described in the section above.

hub_path <- here::here("sismid-ili-forecasting-sandbox")
hub_con <- connect_hub(hub_path)

We can look at one forecast from the hub to get a sense of what a forecast should look like:

hub_con |> 
  filter(model_id == "delphi-epicast",
         origin_date == "2015-11-14") |> 
  collect_hub()

# A tibble: 1,012 × 9
   model_id     origin_date location target horizon target_end_date output_type output_type_id value
 * <chr>        <date>      <chr>    <chr>    <int> <date>          <chr>                <dbl> <dbl>
 1 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.01  0.504
 2 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.025 0.512
 3 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.05  0.525
 4 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.1   0.551
 5 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.15  0.578
 6 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.2   0.604
 7 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.25  0.630
 8 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.3   0.657
 9 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.35  0.683
10 delphi-epic… 2015-11-14  HHS Reg… ili p…       1 2015-11-21      quantile             0.4   0.709
# ℹ 1,002 more rows

Building forecasts for multiple locations

Let’s start as we did in other sessions by downloading the data that we want to work with.

First, we define some location information:

locs <-  c("nat", 
           "hhs1", "hhs2", "hhs3", "hhs4", "hhs5", 
           "hhs6", "hhs7", "hhs8", "hhs9", "hhs10")
location_formal_names <- c("US National", paste("HHS Region", 1:10))
loc_df <- data.frame(region = locs, location = location_formal_names)

The following code will download data from epidatr, but we have also downloaded these data as a static data object for reproducibility.

flu_data_hhs <- pub_fluview(regions = locs, 
                        epiweeks = epirange(200335, 202035)) |> 
  ## aligning with the notation and dates for the hub
  ## in epidatr, weeks are labeled with the first day of the epiweek
  ## in the hub, weeks are labeled with the last day of the epiweek
  mutate(origin_date = epiweek + 6) |> 
  left_join(loc_df) |> 
  select(location, origin_date, wili) |>
  as_tsibble(index = origin_date, key = location)

The following code will load the saved data into the session:

data(flu_data_hhs)

Building a simple forecast for multiple locations

We’re going to start by building a simple auto-regressive forecast without putting too much thought into which exact model specification we choose. Let’s try, for starters, an ARIMA(2,1,0) with a fourth-root transformation. Recall, that this means we are using two auto-regressive terms with one set of differencing of the data.

fit_arima210 <- flu_data_hhs |>
  filter(origin_date <= "2015-10-24") |> 
  model(ARIMA(my_fourth_root(wili) ~ pdq(2,1,0)))
forecast(fit_arima210, h=4) |>
  autoplot(flu_data_hhs |> 
             filter(
               origin_date <= "2015-10-24",
               origin_date >= "2014-09-01"
             )) +
  facet_wrap(.~location, scales = "free_y") +
  labs(y = "% of visits due to ILI",
       x = NULL)

Note that if you set up your tsibble object correctly, with the "location" column as a key variable, then fable will fit and forecast 11 separate versions of the ARIMA(2,1,0) model for you directly, one for each of the locations in the dataset.

Each model can be inspected as before

fit_arima210 |> 
  filter(location == "HHS Region 1") |> 
  report(model)

Series: wili 
Model: ARIMA(2,1,0) 
Transformation: my_fourth_root(wili) 

Coefficients:
         ar1     ar2
      0.0302  0.0354
s.e.  0.0397  0.0403

sigma^2 estimated as 0.004699:  log likelihood=800.67
AIC=-1595.33   AICc=-1595.29   BIC=-1581.98

fit_arima210 |> 
  filter(location == "HHS Region 9") |> 
  report(model)

Series: wili 
Model: ARIMA(2,1,0) 
Transformation: my_fourth_root(wili) 

Coefficients:
          ar1      ar2
      -0.0938  -0.0118
s.e.   0.0402   0.0402

sigma^2 estimated as 0.004732:  log likelihood=798.43
AIC=-1590.86   AICc=-1590.82   BIC=-1577.51

Ideas for more complex models

Below, you will be encouraged to try fitting your own models. Here are some ideas of other models you could try:

use the fable::VAR() model to build a multivariate vector auto-regressive model.
use the fable::NNETAR() model to build a neural network model.
add forecasts from the renewal model (this may take some time to run)
add forecasts from other variants of the fable::ARIMA() models
try out other specifications of transformations (e.g. a log transformation instead of a fourth-root) or the number of fourier terms.
incorporate another model of your choice.
consider working with a generative AI coding tool to build a model. If you are familiar with tools like Claude, Codex or other AI-based coding assistants, they might be very effective and taking an idea of a model that you’d like to translate into a forecast model and building a prototype quickly that would work in the framework of this course.
build some ensembles of all (or a subset) of the models you created.

Designing forecast validation and testing

There are five seasons worth of weeks that we can make forecasts for in our sandbox hub. We will design an validation and testing experiment where we will use the first two seasons (2015/2016 and 2016/2017) as our “validation set”. We can fit as many models as we want to on the weeks during these two seasons. However, to show that we have a modeling process that can create a reliable forecasting model, we want to pass only 1-2 models on to the test phase to measure their performance. We want to pass along only models that we think will perform well in the test phase.

The figure below shows the validation and testing split. The first two validation seasons are used to tune the models and select ones that you like best. The next three testing seasons are used to measure performance of any selected models that are passed to that phase. Note that we intentionally didn’t show the data from the three testing seasons so as not to bias you beforehand by getting a glimpse of the dynamics on those seasons.

And the figure below here shows the expanding window cross-validation where for each week that a forecast is made (rows of the figure), the blue squares indicate weeks whose data are used in training and the red squares indicate the weeks for which forecasts are made. Note that the training data actually extends back to 2003, but we are only showing some of those weeks for easier visualization.

Real-time data and reporting backfill

When these seasons were forecast in real time, forecasters did not have the finalized wILI values plotted above. ILINet is revised for weeks after it is first reported, as additional outpatient-provider reports arrive — a phenomenon known as backfill. The course package includes flu_data_hhs_versions, which records every reported version of each week’s wILI, with the as_of column giving the date each version was published.

data(flu_data_hhs_versions)

flu_data_hhs_versions |>
  filter(location == "US National",
         origin_date %in% as.Date(c("2016-01-02", "2016-01-30", "2016-02-27"))) |>
  ggplot(aes(as_of, wili, colour = factor(origin_date))) +
  geom_line() +
  geom_point(size = 1) +
  labs(x = "data version (as_of date)",
       y = "% of visits due to ILI",
       colour = "week ending",
       title = "Reporting backfill: revisions to US National wILI")

Each line shows how the estimate for a single week changed as later releases revised it. The leftmost point is the first reported value — what a forecaster targeting the following weeks would actually have seen — and it often drifts substantially before settling.

Another way to see backfill is to overlay several versions of the recent series — how the last weeks of data looked at each successive weekly release around the 2015/2016 peak.

win_start <- as.Date("2015-12-01")

focus_versions <- flu_data_hhs_versions |>
  filter(location == "US National",
         as_of >= as.Date("2015-12-01"), as_of <= as.Date("2016-03-12")) |>
  distinct(as_of) |>
  arrange(as_of) |>
  pull(as_of)

cols <- scales::seq_gradient_pal("blue", "red", "Lab")(seq(0,1,length.out=length(focus_versions)))

flu_data_hhs_versions |>
  filter(location == "US National",
         as_of %in% focus_versions,
         origin_date >= win_start, origin_date <= as_of) |>
  ggplot(aes(origin_date, wili, colour = factor(as_of), group = as_of)) +
  geom_line() +
  geom_point(size = 1) +
  labs(x = NULL, y = "% of visits due to ILI",
       colour = "data version\n(as_of)",
       title = "US National wILI near the 2015/2016 peak, by data version") +
  scale_color_manual(values = cols)

Each line is the series as it stood at one release date, showing only the most recent weeks. The right-most point of each line is a week’s first-reported value; following that same week across the later releases shows how it was subsequently revised. A forecaster working in late January had only one of the purple lines to go on — the later weeks had not been reported yet, and even the reported ones would still later be corrected.

To make our retrospective experiment realistic, the cross-validation datasets used below are built from these vintage snapshots: for each forecast date, a split contains only the data that had actually been reported by then.

Validation-phase forecasts

We’ll follow the same time-series cross-validation recipe that we used in the session on forecast evaluation to run multiple forecasts at once for our validation phase.

You could adapt the code below to run validation-phase forecasts.

Here are all the valid origin_dates that we could submit forecasts for, read directly from the hub:

origin_dates <- hub_path |>
  read_config("tasks") |>
  get_round_ids()
origin_dates

  [1] "2015-10-17" "2015-10-24" "2015-10-31" "2015-11-07" "2015-11-14" "2015-11-21" "2015-11-28"
  [8] "2015-12-05" "2015-12-12" "2015-12-19" "2015-12-26" "2016-01-02" "2016-01-09" "2016-01-16"
 [15] "2016-01-23" "2016-01-30" "2016-02-06" "2016-02-13" "2016-02-20" "2016-02-27" "2016-03-05"
 [22] "2016-03-12" "2016-03-19" "2016-03-26" "2016-04-02" "2016-04-09" "2016-04-16" "2016-04-23"
 [29] "2016-04-30" "2016-05-07" "2016-10-22" "2016-10-29" "2016-11-05" "2016-11-12" "2016-11-19"
 [36] "2016-11-26" "2016-12-03" "2016-12-10" "2016-12-17" "2016-12-24" "2016-12-31" "2017-01-07"
 [43] "2017-01-14" "2017-01-21" "2017-01-28" "2017-02-04" "2017-02-11" "2017-02-18" "2017-02-25"
 [50] "2017-03-04" "2017-03-11" "2017-03-18" "2017-03-25" "2017-04-01" "2017-04-08" "2017-04-15"
 [57] "2017-04-22" "2017-04-29" "2017-05-06" "2017-10-21" "2017-10-28" "2017-11-04" "2017-11-11"
 [64] "2017-11-18" "2017-11-25" "2017-12-02" "2017-12-09" "2017-12-16" "2017-12-23" "2017-12-30"
 [71] "2018-01-06" "2018-01-13" "2018-01-20" "2018-01-27" "2018-02-03" "2018-02-10" "2018-02-17"
 [78] "2018-02-24" "2018-03-03" "2018-03-10" "2018-03-17" "2018-03-24" "2018-03-31" "2018-04-07"
 [85] "2018-04-14" "2018-04-21" "2018-04-28" "2018-05-05" "2018-10-13" "2018-10-20" "2018-10-27"
 [92] "2018-11-03" "2018-11-10" "2018-11-17" "2018-11-24" "2018-12-01" "2018-12-08" "2018-12-15"
 [99] "2018-12-22" "2018-12-29" "2019-01-05" "2019-01-12" "2019-01-19" "2019-01-26" "2019-02-02"
[106] "2019-02-09" "2019-02-16" "2019-02-23" "2019-03-02" "2019-03-09" "2019-03-16" "2019-03-23"
[113] "2019-03-30" "2019-04-06" "2019-04-13" "2019-04-20" "2019-04-27" "2019-05-04" "2019-10-12"
[120] "2019-10-19" "2019-10-26" "2019-11-02" "2019-11-09" "2019-11-16" "2019-11-23" "2019-11-30"
[127] "2019-12-07" "2019-12-14" "2019-12-21" "2019-12-28" "2020-01-04" "2020-01-11" "2020-01-18"
[134] "2020-01-25" "2020-02-01" "2020-02-08" "2020-02-15" "2020-02-22" "2020-02-29"

We provide a ready-made time-series cross-validation dataset for each season in the course package (flu_data_hhs_tscv_season1 through flu_data_hhs_tscv_season5), one row per location, week, and split. Each .split is one forecast date, and its wili values are the vintage values available in real time at that date — so the same week can take different values in different splits, exactly as backfill produced them. The datasets are constructed like a tsibble::stretch_tsibble() expanding window, but with the finalized values swapped for the real-time ones.

Make season 1 forecasts

The cross-validation dataset for the first season (2015/2016) is flu_data_hhs_tscv_season1:

data(flu_data_hhs_tscv_season1)
flu_data_hhs_tscv_season1

# A tsibble: 214,005 x 4 [7D]
# Key:       location, .split [330]
   location     origin_date   wili .split
   <chr>        <date>       <dbl>  <int>
 1 HHS Region 1 2003-08-30  0.559       1
 2 HHS Region 1 2003-09-06  0.423       1
 3 HHS Region 1 2003-09-13  0.0677      1
 4 HHS Region 1 2003-09-20  0.0694      1
 5 HHS Region 1 2003-09-27  0.594       1
 6 HHS Region 1 2003-10-04  0.572       1
 7 HHS Region 1 2003-10-11  0.449       1
 8 HHS Region 1 2003-10-18  0.614       1
 9 HHS Region 1 2003-10-25  0.661       1
10 HHS Region 1 2003-11-01  0.761       1
# ℹ 213,995 more rows

And now we will run all of the forecasts and make 1 through 4 week forecasts, including the generate() function to generate predictive samples and then extracting the quantiles. This code also reformats the fable-style forecasts into the required hubverse structure for this hub.

quantile_levels <- c(0.01, 0.025, seq(0.05, 0.95, 0.05), 0.975, 0.99)

cv_forecasts_season1 <- 
  flu_data_hhs_tscv_season1 |> 
  model(
    arima210 = ARIMA(my_fourth_root(wili) ~ pdq(2,1,0))
  ) |> 
  generate(h = 4, times = 100, bootstrap = TRUE) |> 
  ## the following 3 lines of code ensure that there is a horizon variable in the forecast data
  group_by(.split, .rep, location, .model) |> 
  mutate(horizon = row_number()) |> 
  ungroup() |> 
  as_tibble() |> 
  ## make hubverse-friendly names
  rename(
    target_end_date = origin_date,
    value = .sim,
    model_id = .model
  ) |> 
  left_join(loc_df) |> 
  mutate(origin_date = target_end_date - horizon * 7L) |> 
  ## compute the quantiles
  group_by(model_id, location, origin_date, horizon, target_end_date) |> 
  reframe(tibble::enframe(quantile(value, quantile_levels), "quantile", "value")) |> 
  mutate(output_type_id = as.numeric(stringr::str_remove(quantile, "%"))/100, .keep = "unused",
         target = "ili perc",
         output_type = "quantile",
         model_id = "sismid-arima210")

Joining with `by = join_by(location)`

Tip

Note that you need to choose a model_id for your model. It is required by our hub to have the format of [team abbreviation]-[model abbreviation]. When you are making a local hub, you should pick a team abbreviation that is something short and specific to you. For now, we’re going to use sismid as the team abbreviation. (See last line of code above for where I create the name.)

The above table has 30360 rows, which makes sense because there are

30 origin dates
11 locations
4 horizons
23 quantiles

And \(30 \cdot 11 \cdot 4 \cdot 23 = 30,360\).

Let’s plot one set of forecasts as a visual sanity check to make sure our reformatting has worked. The coloured observations are the data as it was available on the forecast date (2015-12-19) — what the model was actually working from — while the grey line shows the finalized data for the whole season, i.e. what eventually happened. Where the two diverge in the weeks just before the forecast, you are seeing reporting backfill.

## observed data as it was available on the 2015-12-19 forecast date (vintage)
obs_asof <- flu_data_hhs_tscv_season1 |>
  tibble::as_tibble() |>
  group_by(.split) |>
  filter(max(origin_date) == as.Date("2015-12-19")) |>
  ungroup() |>
  filter(origin_date >= as.Date("2015-10-01")) |>
  rename(observation = wili)

## finalized data for the whole season, shown in grey for context
finalized_season <- flu_data_hhs |>
  tibble::as_tibble() |>
  filter(origin_date >= as.Date("2015-10-01"), origin_date <= as.Date("2016-05-15"))

cv_forecasts_season1 |>
  filter(origin_date == "2015-12-19") |>
  plot_step_ahead_model_output(
    obs_asof,
    x_target_col_name = "origin_date",
    x_col_name = "target_end_date",
    use_median_as_point = TRUE,
    facet = "location",
    facet_scales = "free_y",
    facet_nrow = 3,
    interactive = FALSE
  ) +
  geom_line(
    data = finalized_season,
    aes(x = origin_date, y = wili),
    colour = "grey60", linewidth = 0.4, inherit.aes = FALSE
  ) +
  theme(legend.position = "bottom")

Warning: ! `model_out_tbl` must be a `model_out_tbl`. Class applied by default

Make season 2 forecasts

The cross-validation dataset for the second season (2016/2017) is flu_data_hhs_tscv_season2:

data(flu_data_hhs_tscv_season2)
flu_data_hhs_tscv_season2

# A tsibble: 223,619 x 4 [7D]
# Key:       location, .split [319]
   location     origin_date   wili .split
   <chr>        <date>       <dbl>  <int>
 1 HHS Region 1 2003-08-30  0.559       1
 2 HHS Region 1 2003-09-06  0.423       1
 3 HHS Region 1 2003-09-13  0.0677      1
 4 HHS Region 1 2003-09-20  0.0694      1
 5 HHS Region 1 2003-09-27  0.594       1
 6 HHS Region 1 2003-10-04  0.572       1
 7 HHS Region 1 2003-10-11  0.449       1
 8 HHS Region 1 2003-10-18  0.614       1
 9 HHS Region 1 2003-10-25  0.661       1
10 HHS Region 1 2003-11-01  0.761       1
# ℹ 223,609 more rows

Generate and reformat season 2 forecasts:

cv_forecasts_season2 <- 
  flu_data_hhs_tscv_season2 |> 
  model(
    arima210 = ARIMA(my_fourth_root(wili) ~ pdq(2,1,0))
  ) |> 
  generate(h = 4, times = 100, bootstrap = TRUE) |> 
  ## the following 3 lines of code ensure that there is a horizon variable in the forecast data
  group_by(.split, .rep, location, .model) |> 
  mutate(horizon = row_number()) |> 
  ungroup() |> 
  as_tibble() |> 
  ## make hubverse-friendly names
  rename(
    target_end_date = origin_date,
    value = .sim,
    model_id = .model
  ) |> 
  left_join(loc_df) |> 
  mutate(origin_date = target_end_date - horizon * 7L) |> 
  ## compute the quantiles
  group_by(model_id, location, origin_date, horizon, target_end_date) |> 
  reframe(tibble::enframe(quantile(value, quantile_levels), "quantile", "value")) |> 
  mutate(output_type_id = as.numeric(stringr::str_remove(quantile, "%"))/100, .keep = "unused",
         target = "ili perc",
         output_type = "quantile",
         model_id = "sismid-arima210")

Joining with `by = join_by(location)`

Save the forecasts locally

Now let’s save the validation forecasts locally so that we can incorporate these forecasts with the existing forecasts in the hub.

First, we need to create a model metadata file (formatted as a YAML file) and store it in the hub. This can be something very minimal, but there are a few required pieces. Here is some simple code to write out a minimal model metadata file.

this_model_id <- "sismid-arima210"

metadata_filepath <- file.path(
  hub_path,
  "model-metadata", 
  paste0(this_model_id, ".yml"))

my_text <- c("team_abbr: \"sismid\"", 
             "model_abbr: \"arima210\"", 
             "designated_model: true")

writeLines(my_text, metadata_filepath)

And now we can write the files out, one for each origin_date.

# Group the forecasts by task id variables
groups <- bind_rows(cv_forecasts_season1, cv_forecasts_season2) |> 
  group_by(model_id, target, origin_date) |> 
  group_split()

# Save each group as a separate CSV
for (i in seq_along(groups)) {
  group_df <- groups[[i]]
  this_model_id <- group_df$model_id[1]
  this_origin_date <- group_df$origin_date[1]
  
  ## remove model_id from saved data, as it is implied from filepath
  group_df <- select(group_df, -model_id)
  
  ## path to the file from the working directory of the instructional repo
  model_folder <- file.path(
    hub_path,
    "model-output", 
    this_model_id)

  ## just the filename, no path
  filename <- paste0(this_origin_date, "-", this_model_id, ".csv")
  
  ## path to the file
  results_path <- file.path(
    model_folder, 
    filename)
  
  ## if this model's model-out directory doesn't exist yet, make it
  if (!file.exists(model_folder)) {
    dir.create(model_folder, recursive = TRUE)
  }
  
  write.csv(group_df, file = results_path, row.names = FALSE)
  
  ### if you run into errors, the code below could help trouble-shoot validations
  # hubValidations::validate_submission(
  #   hub_path = hub_path,
  #   file_path = file.path(this_model_id, filename)
  # )
}

Local hub evaluation

Now that you have your model’s forecasts in the local hub, you should be able to run an evaluation to compare the model’s forecasts to the model in the hub.

We first collect the forecasts from the hub (including our new model)

validation_origin_dates <- origin_dates[which(as.Date(origin_dates) <= as.Date("2017-05-06"))]

new_hub_con <- connect_hub(hub_path)

validation_forecasts <- new_hub_con |> 
  filter(origin_date %in% validation_origin_dates) |> 
  collect_hub()

and then compare them to the oracle output.

oracle_output <- connect_target_oracle_output(hub_path) |> 
  collect()

hubEvals::score_model_out(
  validation_forecasts,
  oracle_output
) |> 
  arrange(wis) |> 
  knitr::kable(digits = 2)

model_id	wis	overprediction	underprediction	dispersion	bias	interval_coverage_50	interval_coverage_90	ae_median
delphi-epicast	0.31	0.09	0.08	0.14	0.06	0.47	0.90	0.43
sismid-arima210	0.37	0.13	0.11	0.14	0.02	0.44	0.89	0.58
hist-avg	0.45	0.05	0.16	0.23	-0.14	0.55	0.94	0.68
kot-kot	2.49	1.43	0.00	1.06	0.71	0.13	0.97	4.46

Here we can see that our new sismid-arima210 model has lower probabilistic accuracy (higher WIS) than the delphi-epicast model, and has better accuracy than the hist-avg model. Also, it shows less bias than hist-avg or delphi-epicast, although a bit more dispersion (it is less sharp) than the delphi-epicast model.

Tip

Our cross-validation datasets use the vintage data available at each forecast date — the same real-time data (including backfill) that the pre-existing hub models had to work with. That makes this a fairer comparison than if we had trained on finalized data. The scoring truth (the oracle output) is, by contrast, each season’s final settled value (the data as of the following summer): we judge every model against what actually happened, but without letting later cross-season revisions leak in.

Going further

Challenge

If you complete this session, consider generating more model forecasts for the first two seasons. You could also move on to the next test-phase local hub session.

Wrap up

Review what you’ve learned in this session with the learning objectives

Introduction

Slides

Objectives

Source file

Libraries used

Initialisation

A flu forecasting hub for model dev

ILI Forecasting at HHS Region level

Details of this hub

Building forecasts for multiple locations

Building a simple forecast for multiple locations

Ideas for more complex models

Designing forecast validation and testing

Real-time data and reporting backfill

Validation-phase forecasts

Make season 1 forecasts

Make season 2 forecasts

Save the forecasts locally

Local hub evaluation

Going further

Challenge

Wrap up

References