Evaluating real-world ensembles

Introduction

This session explicitly synthesizes content from two previous sessions, as we will be building and evaluating ensembles using real-world forecast data. In particular, we will explore weighting models in the ensemble.

Slides

• Real-world ensembles

Objectives

The aim of this session is to develop advanced skills in working with ensemble forecasts and an appreciation for the challenges improving on equal-weight ensembles by estimating model weights.

NoteSetup

Source file

The source file of this session is located at sessions/eval-real-world-ensembles.qmd.

Libraries used

In this session we will use the nfidd package to access some stored datasets, the dplyr package for data wrangling, the ggplot2 library for plotting, and the following hubverse packages: hubData, hubUtils, hubEvals and hubEnsembles.

library("nfidd.forecasting")
library("dplyr")
library("ggplot2")
library("hubData")
library("hubUtils")
library("hubEvals")
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(78) # for Greg Maddux

Downloading the forecast data

We will start by accessing forecast data from the COVID-19 Forecast Hub as we did in the earlier session.

Accessing forecast data from the cloud

The following code can be used to retrieve COVID-19 forecasts from the cloud. To make these results more reproducible and to provide a snapshot of the data, we also provide the forecasts as a data object that can be loaded directly from the nfidd package, although the following code could be run or adapted for filtering different subsets of data (see “Going Further” section.

hub_path_cloud <- hubData::s3_bucket("s3://covid19-forecast-hub")
hub_con_cloud <- hubData::connect_hub(hub_path_cloud, skip_checks = TRUE)

data(covid_locations)

covid_forecasts <- hub_con_cloud |> 
  filter(
    reference_date <= as.Date("2025-07-12"),
    output_type == "quantile",
    target == "wk inc covid hosp",
    horizon >= 0
    ) |> 
  collect() |> 
  left_join(covid_locations)

This line of code will load all the forecasts saved in the package into your R session.

data(covid_forecasts)

Accessing target data from the cloud

Here is code to query the time-series data and oracle output data from the cloud-based hub.

covid_time_series <- connect_target_timeseries(hub_path_cloud) |> 
  filter(target == "wk inc covid hosp") |> 
  collect() |> 
  left_join(covid_locations)

hub_path_cloud <- hubData::s3_bucket("s3://covid19-forecast-hub")
covid_oracle_output <- connect_target_oracle_output(hub_path_cloud) |> 
  collect()
covid_oracle_output

# A tibble: 18,270 × 6
   location target_end_date target       output_type oracle_value output_type_id
   <chr>    <date>          <chr>        <chr>              <dbl> <chr>         
 1 US       2024-11-09      wk inc covi… quantile            7988 <NA>          
 2 01       2024-11-09      wk inc covi… quantile             115 <NA>          
 3 02       2024-11-09      wk inc covi… quantile              14 <NA>          
 4 04       2024-11-09      wk inc covi… quantile             532 <NA>          
 5 05       2024-11-09      wk inc covi… quantile              60 <NA>          
 6 06       2024-11-09      wk inc covi… quantile             643 <NA>          
 7 08       2024-11-09      wk inc covi… quantile             145 <NA>          
 8 09       2024-11-09      wk inc covi… quantile              82 <NA>          
 9 10       2024-11-09      wk inc covi… quantile              20 <NA>          
10 11       2024-11-09      wk inc covi… quantile               9 <NA>          
# ℹ 18,260 more rows

However, to ensure reproducibility (since this is a live dataset), we have downloaded this object for you already as of July 9, 2025, and it is available to load in from the course R package directly.

data(covid_time_series)

Ensembling forecasts

Recall that in the earlier session we analyzed and compared individual forecasts from the COVID-19 forecast hub, including the CovidHub-ensemble model put together by the hub administrators. In this session, we will build several new ensemble forecasts and compare these to the CovidHub-ensemble model. Based on the last results from the section above, we saw that there was some fairly consistent performance from a few of the models. This is one precondition for weighted ensembles to have a chance at performing well.

We will use ensembling methods and tools presented in the Forecast Ensembles session. Here are the different approaches we will use:

• lop_unweighted_all: an unweighted linear opinion pool (LOP) of all available forecasts each week.
• lop_unweighted_select: an unweighted LOP of the most commonly submitted forecasts.
• lop_inv_wis_weighted_select: a weighted LOP of models that submitted more than 60% of the time.
• median_inv_wis_weighted_select: a weighted median ensemble of models that submitted more than 60% of the time.

Note that the existing CovidHub-ensemble forecast is basically median-unweighted-all in the above nomenclature. It would be hard to do any weighted-all ensembles since when models are missing a lot of forecasts it might be hard to estimate an accurate weight for them.

Finding frequently submitted models

Let’s start by identifying the subset of models that submitted more than 60% of the time.

## 33 dates, 4 horizons, 53 locations, 60%
threshold_60pct <- 33 * 4 * 53 * 0.6
model_subset <- covid_forecasts |> 
  filter(output_type_id == 0.5) |> 
  group_by(model_id) |> 
  summarize(targets = n()) |> 
  filter(targets > threshold_60pct) |> 
  pull(model_id)
model_subset

 [1] "CEPH-Rtrend_covid"        "CMU-TimeSeries"          
 [3] "CMU-climate_baseline"     "CovidHub-baseline"       
 [5] "CovidHub-ensemble"        "MOBS-GLEAM_COVID"        
 [7] "NEU_ISI-AdaptiveEnsemble" "OHT_JHU-nbxd"            
 [9] "UM-DeepOutbreak"          "UMass-ar6_pooled"        
[11] "UMass-gbqr"              

Build unweighted LOP ensembles

Let’s start with the “easy” part, building the unweighted ensembles.

TipLOP details

Note that when we build LOPs with quantile based forecasts, the algorithm needs to approximate a full distribution (the tails are not well-defined in a quantile representation). There are some assumptions made about how this is approximated. For the really gory details, you can check out the hubEnsembles documentation

Note that trying to build all of the forecasts in one call to linear_pool() exceeded the memory on your instructors laptop, so we will do a purrr::map_dfr() “loop” to do this one reference_date at a time. (This code takes a few minutes to run.)

# Get a list of unique reference dates, removing the last one
reference_dates <- covid_forecasts |> 
  filter(reference_date != as.Date("2025-07-12")) |> 
  pull(reference_date) |> 
  unique()

# Map over each date and apply linear_pool separately
lop_unweighted_all <- purrr::map_dfr(
  reference_dates, 
  function(date) {
    covid_forecasts |> 
      filter(model_id != "CovidHub-ensemble", 
             reference_date == date) |> 
      hubEnsembles::linear_pool(model_id = "lop_unweighted_all")
})


lop_unweighted_select <- purrr::map_dfr(
  reference_dates, 
  function(date) {
    covid_forecasts |> 
      filter(model_id %in% model_subset,
             model_id != "CovidHub-ensemble", 
             reference_date == date) |> 
      hubEnsembles::linear_pool(model_id = "lop_unweighted_select")
})

Build weighted ensemble

Inverse WIS weighting is a simple method that weights the forecasts by the inverse of their WIS over some period. Note that identifying what this period should be in order to produce the best forecasts is not straightforward as predictive performance may vary over time if, say, some models are better at making predictions during different parts of an outbreak. For example, some models might be good at predicting exponential growth at the beginning of an outbreak but not as good at seeing the slowdown.

The main benefit of WIS weighting over other methods is that it is simple to understand and implement. However, it does not optimise the weights directly to produce the best forecasts. It relies on the hope that giving more weight to better performing models yields a better ensemble.

We will use the WIS scores that we have computed above to generate weights for our ensembles. The weights for one reference_date will be based on the inverse of the cumulative average WIS achieved by the model up to the week prior to the current reference_date.

Caution

Note that since our scores are computed based on data available at the end of the season, this is kind of “cheating”, since we’re using information that would not have been available to us in real time. If we were to do this entirely correctly, we would, for each week, compute the scores based on the available data at that time. But, for instructional purposes, this will suffice.

Let’s start by calculating WIS values by model and target_end_date. This will allow us to see what scores are available for a model’s forecasts that it had for dates prior to the current forecast date.

scores_by_target_date <- covid_forecasts |> 
  filter(model_id %in% model_subset,
         location != "US") |> 
  score_model_out(
    covid_oracle_output,
    metrics = "wis",
    by = c("model_id", "target_end_date")
  )
scores_by_target_date

# A tibble: 396 × 3
   model_id          target_end_date   wis
   <chr>             <date>          <dbl>
 1 CEPH-Rtrend_covid 2024-11-23       27.1
 2 CEPH-Rtrend_covid 2024-11-30       29.0
 3 CEPH-Rtrend_covid 2024-12-07       27.7
 4 CEPH-Rtrend_covid 2024-12-14       33.1
 5 CEPH-Rtrend_covid 2024-12-21       38.8
 6 CEPH-Rtrend_covid 2024-12-28       64.2
 7 CEPH-Rtrend_covid 2025-01-04       92.7
 8 CEPH-Rtrend_covid 2025-01-11       53.2
 9 CEPH-Rtrend_covid 2025-01-18       63.8
10 CEPH-Rtrend_covid 2025-01-25       87.7
# ℹ 386 more rows

Now, we will compute a cumulative average WIS up to a given date for all models except the ensemble, which we are not going to include in our ensemble.

cum_wis <- scores_by_target_date |>
  filter(model_id != "CovidHub-ensemble") |> 
  arrange(model_id, target_end_date) |> 
  group_by(model_id) |> 
  mutate(
    cumulative_mean_wis = cummean(wis), ## cumulative mean WIS including current week
    cumulative_mean_wis_lag1 = lag(cumulative_mean_wis)
    ) |> 
  ungroup()
cum_wis

# A tibble: 359 × 5
   model_id     target_end_date   wis cumulative_mean_wis cumulative_mean_wis_…¹
   <chr>        <date>          <dbl>               <dbl>                  <dbl>
 1 CEPH-Rtrend… 2024-11-23       27.1                27.1                   NA  
 2 CEPH-Rtrend… 2024-11-30       29.0                28.1                   27.1
 3 CEPH-Rtrend… 2024-12-07       27.7                28.0                   28.1
 4 CEPH-Rtrend… 2024-12-14       33.1                29.2                   28.0
 5 CEPH-Rtrend… 2024-12-21       38.8                31.1                   29.2
 6 CEPH-Rtrend… 2024-12-28       64.2                36.7                   31.1
 7 CEPH-Rtrend… 2025-01-04       92.7                44.7                   36.7
 8 CEPH-Rtrend… 2025-01-11       53.2                45.7                   44.7
 9 CEPH-Rtrend… 2025-01-18       63.8                47.7                   45.7
10 CEPH-Rtrend… 2025-01-25       87.7                51.7                   47.7
# ℹ 349 more rows
# ℹ abbreviated name: ¹​cumulative_mean_wis_lag1

Tip

Note that we also compute a cumulative_mean_wis_lag1, which is the mean up through the previous week. This is to inject a bit of realism: how could we know in real-time how our forecast did for the current week?

So the cumulative_mean_wis_lag1 is the cumulative mean WIS that we could have observed at the target_end_date. (This all still assumes that whatever data we observed at the target_end_date was complete and never revised.)

We will proceed by computing the inverse mean WIS and using those to derive the weights.

weights_per_model <- cum_wis |> 
  mutate(inv_wis = 1 / cumulative_mean_wis_lag1) |>
  group_by(target_end_date) |> 
  mutate(inv_wis_total_by_date = sum(inv_wis, na.rm = TRUE)) |>
  ungroup() |> 
  mutate(weight = inv_wis / inv_wis_total_by_date) |> ## this normalises the weights to sum to 1
  select(model_id, target_end_date, weight)

And we will put in place a few small fixes to the weights:

• assign equal weights for the first two weeks. (There are missing forecasts here, and we haven’t learned much about the model performance yet anyways.)
• change target_end_date to reference_date because now we want the weights computed up to that target_end_date to be applied on that reference_date.
• set to 0 a missing weight for the NEU_ISI-AdaptiveEnsemble.

## assign equal weights for the first week
first_week_idx <- weights_per_model$target_end_date == as.Date("2024-11-23")
weights_per_model[first_week_idx, "weight"] <- 1/sum(first_week_idx)

second_week_idx <- weights_per_model$target_end_date == as.Date("2024-11-30")
weights_per_model[second_week_idx, "weight"] <- 1/sum(second_week_idx)

weights_per_model <- rename(weights_per_model, reference_date = target_end_date)

weights_per_model$weight[is.na(weights_per_model$weight)] <- 0

As always, it is good to plot your data as a sanity check. Here are the model weights plotted over time. For each week (x-axis), the weights assigned to each model are shown in a stacked bar plot, with the height of the bar equal to the weight of the model.

plot_weights <- ggplot(weights_per_model) +
  geom_col(aes(x = reference_date, y = weight, fill = model_id))
plotly::ggplotly(plot_weights)

TipTake 5 minutes

Looking at the above plot, what is the “story” that the weights tell when plotted over time? Which models did the ensemble like and dislike over time and how did it change? (Remember, in the first two weeks, all models received equal weight.)

NoteSolution

Starting in the third week, three models were assigned less weight than others: CMU-climate-baseline, CMU-TimeSeries and COVIDHub-baseline. Over time, CMU-TimeSeries and COVIDHub-baseline earned back some “trust” with higher weights, but CMU-climate-baseline stayed with a low weight.

In general, the weights did not change much after the first few weeks (except in the last week when UM-DeepOutbreak wasn’t submitted.)

Now that we have the weights and we’ve checked them, we can build some ensembles. First, we will build the median ensembles, both weighted and unweighted.

## build the weighted ensemble
median_inv_wis_weighted_select <- covid_forecasts |> 
  filter(model_id %in% model_subset,
         model_id != "CovidHub-ensemble") |> 
  simple_ensemble(weights = weights_per_model, 
                  agg_fun = median,
                  model_id = "median_inv_wis_weighted_select")

## build unweighted select median ensemble
median_unweighted_select <- covid_forecasts |> 
  filter(model_id %in% model_subset,
         model_id != "CovidHub-ensemble") |> 
  simple_ensemble(agg_fun = median,
                  model_id = "median_unweighted_select")

And then we will build the inverse WIS weighted LOP ensemble.

Caution

It takes several minutes for the LOP ensembles to be built if you are running the code on your machine.

lop_inv_wis_weighted_select <- purrr::map_dfr(
  reference_dates, 
  function(date) {
    covid_forecasts |> 
      filter(model_id %in% model_subset,
             model_id != "CovidHub-ensemble", 
             reference_date == date) |> 
      linear_pool(model_id = "lop_inv_wis_weighted_select",
                  weights = weights_per_model)
})

Visually compare ensembles

As always, it is really important to look at the forecasts to develop some intuition about the similarities and differences between the models. The following code subsets the output to look just at the ensembles (including the original CovidHub-ensemble) and then plots forecasts for three locations (GA, MA and CA) at one date. Feel free to modify the code to look at other locations or dates.

covid_forecasts_w_ensembles <- bind_rows(
  covid_forecasts,
  lop_unweighted_all,
  lop_unweighted_select,
  lop_inv_wis_weighted_select,
  median_inv_wis_weighted_select,
  median_unweighted_select
)

model_subset_ens <- c("CovidHub-ensemble", 
                      "lop_unweighted_all",
                      "lop_unweighted_select",
                      "median_inv_wis_weighted_select",
                      "median_unweighted_select",
                      "lop_inv_wis_weighted_select")

covid_forecasts_w_ensembles |> 
  filter(
    reference_date %in% as.Date("2025-01-18"),
    abbreviation %in% c("GA", "MA", "CA"),
    model_id %in% model_subset_ens
  ) |> 
  hubVis::plot_step_ahead_model_output(
    target_data = covid_time_series |> 
      filter(as_of == as.Date("2025-03-05"),
             abbreviation %in% c("GA", "MA", "CA")),
    use_median_as_point = TRUE,
    facet = "abbreviation",
    facet_nrow = 3,
    facet_scales = "free_y",
    intervals = 0.95,
    x_col_name = "target_end_date",
    x_target_col_name = "date",
    pal_color = "Set3", 
    title = "Weekly hospitalizations due to COVID-19 (data and forecasts)"
  )

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

Warning: ! `output_type_id` column must be a numeric. Converting to numeric.

TipTake 5 minutes

Look at the forecasts above, and possibly explore other dates and/or locations. Are there noticeable differences between the point estimates or between the prediction interval widths? If so, what are they?

NoteSolution

In the plot shown above, there are not noticeable differences between point estimates. All of the LOP models have wider PIs compared to the other ensembles.

Compare all ensembles

We will compare the performance of all ensembles just at the state level based on quantitative scoring rules, since the scores by reference date were only computed at the state level.

covid_forecasts_w_ensembles |> 
  filter(model_id %in% model_subset_ens,
         location != "US") |> 
  score_model_out(
    covid_oracle_output
  ) |> 
  arrange(wis) |> 
  knitr::kable(digits = 2)

model_id	wis	overprediction	underprediction	dispersion	bias	interval_coverage_50	interval_coverage_90	ae_median
median_unweighted_select	26.57	3.67	11.21	11.69	-0.13	0.59	0.92	40.07
CovidHub-ensemble	26.58	3.74	12.20	10.64	-0.18	0.53	0.90	40.40
median_inv_wis_weighted_select	27.61	4.07	12.04	11.50	-0.14	0.57	0.90	41.51
lop_unweighted_all	28.10	2.75	9.68	15.67	-0.12	0.66	0.99	40.40
lop_unweighted_select	28.29	2.71	9.78	15.81	-0.13	0.67	0.99	40.69
lop_inv_wis_weighted_select	28.30	2.73	11.02	14.54	-0.15	0.64	0.97	41.44

Remember, in the above table, the CovidHub-ensemble is kind of like median-unweighted-all because it uses a median ensemble of all models available in every week without weights. Overall, the ensembles have fairly similar average WIS values across all dates and states. Using a median unweighted ensemble of frequently submitting models marginally outperformed the actual CovidHub-ensemble.

TipTake 5 minutes

Why do you think weighting does not improve the forecasts more?

Seeing these patterns, and the plots of the weights above, is there anything you might change if you were given the opportunity to compute the weights in a different way?

Write down your thoughts about these questions and compare answers with a neighbor.

NoteSolution

Why aren’t weighted ensembles better?

When model performance varies a lot from one week to the next, it may be hard to estimate/anticipate/predict which model should get more weight in a given week.

What might your consider changing about the weights?

There isn’t one right answer here, but some things that you might consider to do differently are:

1. Figuring out a way to encourage the model to stay closer to equal weights early on. The weights are being perhaps overly influenced by small sample sizes there.
2. Thinking about using a rolling window of mean WIS to compute weights rather than overall cumulative. This could be good, but also assumes that recent performance will predict future performance, which is not always true!

It does not seem like weighting adds a lot of value to these forecasts, although perhaps with additional data on model performance the weights could be improved.

We can also look at how the ensemble performance varied over time.

ensemble_scores_by_reference_date <- covid_forecasts_w_ensembles |> 
  filter(model_id %in% model_subset_ens,
         location != "US") |> 
  score_model_out(
    covid_oracle_output,
    metrics = "wis",
    by = c("model_id", "reference_date")
  )

p <- ggplot(ensemble_scores_by_reference_date) +
  geom_line(aes(x = reference_date, y = wis, color = model_id))
plotly::ggplotly(p)

It is interesting to note how while the LOP ensembles did slightly worse overall, they actually had slightly more accurate forecasts near the peak of the season, especially in the weeks where the scores were the highest. Overall, though, weighting does not seem to have a large impact on the forecast accuracy.

Going further

Challenge

• This session only includes data from the year leading up to July 2025. Adjust the code to download data from the past year and re-run analyses on the most recent year of forecasts. Consider adding an ensemble where you manually pick a few models and assign them “human judgment” weights. Compare the results from those in the session above.

Methods in practice

Here are a set of papers that talk about attempts to build ensembles from a set of models:

• An effort using FluSight models that generated multiple years of prior forecasts for training the weights. (Reich et al. 2019)
• Building ensembles to predict dengue in Vietnam.(Colón-González et al. 2021)
• Predicting health care demand in France during COVID-19. (Paireau et al. 2022)
• Predicting COVID-19 in Europe using ensembles. (Sherratt et al. 2023)
• Evaluating ensembling strategies for COVID-19 forecasts in the US. (Ray et al. 2023)
• Some work on adaptively estimating weights as functions of other variables. (Ray and Reich 2018; Yamana et al. 2016; McAndrew and Reich 2021)

Wrap up

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

References

Colón-González, Felipe J., Leonardo Soares Bastos, Barbara Hofmann, et al. 2021. “Probabilistic Seasonal Dengue Forecasting in Vietnam: A Modelling Study Using Superensembles.” PLOS Medicine 18 (3): e1003542. https://doi.org/10.1371/journal.pmed.1003542.

McAndrew, Thomas, and Nicholas G. Reich. 2021. “Adaptively Stacking Ensembles for Influenza Forecasting.” Statistics in Medicine 40 (30): 6931–52. https://doi.org/10.1002/sim.9219.

Paireau, Juliette, Alessio Andronico, Nathanaël Hozé, et al. 2022. “An Ensemble Model Based on Early Predictors to Forecast COVID-19 Health Care Demand in France.” Proceedings of the National Academy of Sciences 119 (18): e2103302119. https://doi.org/10.1073/pnas.2103302119.

Ray, Evan L., Logan C. Brooks, Jacob Bien, et al. 2023. “Comparing Trained and Untrained Probabilistic Ensemble Forecasts of COVID-19 Cases and Deaths in the United States.” International Journal of Forecasting 39 (3): 1366–83. https://doi.org/10.1016/j.ijforecast.2022.06.005.

Ray, Evan L., and Nicholas G. Reich. 2018. “Prediction of Infectious Disease Epidemics via Weighted Density Ensembles.” PLOS Computational Biology 14 (2): e1005910. https://doi.org/10.1371/journal.pcbi.1005910.

Reich, Nicholas G., Craig J. McGowan, Teresa K. Yamana, et al. 2019. “Accuracy of Real-Time Multi-Model Ensemble Forecasts for Seasonal Influenza in the U.S.” PLOS Computational Biology 15 (11): e1007486. https://doi.org/10.1371/journal.pcbi.1007486.

Sherratt, Katharine, Hugo Gruson, Rok Grah, et al. 2023. “Predictive Performance of Multi-Model Ensemble Forecasts of COVID-19 Across European Nations.” eLife 12 (April): e81916. https://doi.org/10.7554/eLife.81916.

Yamana, Teresa K., Sasikiran Kandula, and Jeffrey Shaman. 2016. “Superensemble Forecasts of Dengue Outbreaks.” Journal of The Royal Society Interface 13 (123): 20160410. https://doi.org/10.1098/rsif.2016.0410.