From a6c4dfcbf8217db2b7f5cc8c347897ae243a6e3a Mon Sep 17 00:00:00 2001 From: kathsherratt Date: Thu, 9 Jul 2026 19:34:57 +0100 Subject: [PATCH 1/7] add note about quantile forecasts --- NEWS.md | 1 + vignettes/scoring-multivariate-forecasts.Rmd | 47 ++++++++++++++++++++ 2 files changed, 48 insertions(+) diff --git a/NEWS.md b/NEWS.md index e1ee239eb..20c31ed9a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,7 @@ - Fixed `summarise_scores()` producing a data.table with duplicate column names when the input `scores` object had no score columns (e.g. because every metric in `score()` warned and returned nothing). `summarise_scores()` now matches metric columns by exact name rather than regex partial match, and errors with a clear message when there is nothing to summarise (#1179). - Added internal S3 generic `get_forecast_type_ids()` so each forecast type declares the columns (beyond the forecast unit) that identify a unique row. `get_duplicate_forecasts()` now uses this instead of hard-coded column names (#888). - Removed the deprecated vignettes `Deprecated-functions` and `Deprecated-visualisations`. The code for removed functions (`plot_predictions()`, `make_NA()`, `plot_ranges()`, `plot_score_table()`, `merge_pred_and_obs()`) can still be found in the [git history](https://github.com/epiforecasts/scoringutils/tree/d0cd8e2/vignettes) (#1158). +- The vignette "Scoring multivariate forecasts" now explains why the energy score and the variogram score cannot be applied to quantile forecasts. Both scores require samples from the joint predictive distribution, which a set of marginal quantiles does not provide. Treating quantile levels as sample ids silently assumes perfect rank correlation between the targets being forecast jointly. # scoringutils 2.2.0 diff --git a/vignettes/scoring-multivariate-forecasts.Rmd b/vignettes/scoring-multivariate-forecasts.Rmd index ea53a1976..55b5476aa 100644 --- a/vignettes/scoring-multivariate-forecasts.Rmd +++ b/vignettes/scoring-multivariate-forecasts.Rmd @@ -128,4 +128,51 @@ example_mv_point <- as_forecast_multivariate_point( score(example_mv_point) ``` +## What about quantile forecasts? + +Both the energy score and the variogram score need samples drawn from the joint predictive distribution. +Each sample gives one value for every target at once, and it is the variation across samples that tells the score how the targets move together. +Quantile forecasts do not carry this information. +A set of quantiles describes each target on its own and says nothing about the relationship between them. + +It is tempting to line the quantile levels up across targets and treat them as samples, for instance by passing `sample_id = "quantile_level"` to `as_forecast_multivariate_sample()`. +`scoringutils` will warn you about the leftover `quantile_level` column, but it will still return a score. +That score is not the one you want. +Pairing the 5% quantile in one country with the 5% quantile in another assumes the two are perfectly rank correlated, so every quantile forecast ends up scored as though it had predicted perfect dependence, whatever it actually predicted. + +The two scores suffer differently. +The energy score still responds to the marginal distributions, so its values stay in a plausible range, but it loses almost all of its ability to separate forecasts with different correlation structures. +The variogram score fares worse, because dependence is the only thing it measures. +Replace the dependence with an assumption and what is left is close to a function of the observed values alone, which can rank models in the wrong order. + +If you have quantile forecasts and want to score them jointly, you have to supply the missing dependence structure. +You can ask forecasters for samples or trajectories instead of quantiles, which is the route several forecast hubs have taken. +Alternatively you can rebuild a joint distribution from the marginal quantiles using ensemble copula coupling (Schefzik et al., 2013) or the Schaake shuffle (Clark et al., 2004). +Both borrow a rank structure from elsewhere, usually a raw ensemble or the historical observations, so any score you compute afterwards partly reflects that choice. + +Multivariate point forecasts are a different matter. +A point forecast does specify a joint distribution, a degenerate one that places all its mass on a single vector, which is why the variogram score is well defined for it. + If, at any point, you want to score the same forecast using different groupings, you'd have create a new separate forecast object with a different grouping and score that new forecast object. + +## References + +Clark, M., Gangopadhyay, S., Hay, L., Rajagopalan, B. and Wilby, R. (2004). +The Schaake Shuffle: A Method for Reconstructing Space-Time Variability in Forecasted Precipitation and Temperature Fields. +*Journal of Hydrometeorology*, 5(1), 243-262. + + +Gneiting, T., Stanberry, L. I., Grimit, E. P., Held, L. and Johnson, N. A. (2008). +Assessing probabilistic forecasts of multivariate quantities, with an application to ensemble predictions of surface winds. +*TEST*, 17, 211-235. + + +Schefzik, R., Thorarinsdottir, T. L. and Gneiting, T. (2013). +Uncertainty Quantification in Complex Simulation Models Using Ensemble Copula Coupling. +*Statistical Science*, 28(4), 616-640. + + +Scheuerer, M. and Hamill, T. M. (2015). +Variogram-Based Proper Scoring Rules for Probabilistic Forecasts of Multivariate Quantities. +*Monthly Weather Review*, 143(4), 1321-1334. + From ba794586ad656303a4fd26d73975395e8ba0c608 Mon Sep 17 00:00:00 2001 From: kathsherratt Date: Fri, 17 Jul 2026 14:29:02 +0100 Subject: [PATCH 2/7] Update scoring-multivariate-forecasts.Rmd --- vignettes/scoring-multivariate-forecasts.Rmd | 117 +++++++++++++------ 1 file changed, 79 insertions(+), 38 deletions(-) diff --git a/vignettes/scoring-multivariate-forecasts.Rmd b/vignettes/scoring-multivariate-forecasts.Rmd index 55b5476aa..8ac92d00d 100644 --- a/vignettes/scoring-multivariate-forecasts.Rmd +++ b/vignettes/scoring-multivariate-forecasts.Rmd @@ -67,6 +67,8 @@ The corresponding multivariate forecast would similarly specify a predictive dis In the following, let's assume that our samples were draws from a multivariate distribution all along (we just treated them as independent for the univariate case). +> **Note on the demonstration.** The forecasts used in this vignette are for demonstration only and do not represent joint draws of a multivariate distribution. This means the dependence structure scored by multivariate scoring is more an artefact of how the data were prepared than a property of the original forecasts. (See `create-example-data.R` in the `inst/` folder for detail.) + To tell `scoringutils` that we want to treat these as a multivariate forecast, we need to specify the columns that are pooled together to form a single multivariate forecast. We do this via the `joint_across` argument. For example, if we want to pool forecasts across locations and treat them as a single multivariate forecast, we could set `joint_across = c("location", "location_name")` (in our example, the two columns contain essentially the same information - we therefore have to include both in `joint_across` (or could alternatively delete one of them)). ```{r} @@ -107,6 +109,28 @@ score( ) ``` +A set of multivariate targets can be pooled to account for different correlation structures. +If, at any point, you want to score the same forecast using different groupings, you'd have create a new separate forecast object with a different grouping and score that new forecast object. +For example, to pool across horizons, we bring `target_end_date` along with `horizon`, because the two both vary within a trajectory. + +```{r} +example_cases <- na.omit( + example_sample_continuous[ + target_type == "Cases" & + model == "EuroCOVIDhub-ensemble" + ] +) + +example_traj <- as_forecast_multivariate_sample( + data = example_cases, + joint_across = c("horizon", "target_end_date") +) + +head(score(example_traj), 3) +``` + +Each single score now covers one whole trajectory over three horizons. + ## Multivariate point forecasts If you have point forecasts rather than samples, you can score them using the variogram score via `as_forecast_multivariate_point()`. @@ -128,51 +152,68 @@ example_mv_point <- as_forecast_multivariate_point( score(example_mv_point) ``` -## What about quantile forecasts? +## Multivariate quantile forecasts -Both the energy score and the variogram score need samples drawn from the joint predictive distribution. -Each sample gives one value for every target at once, and it is the variation across samples that tells the score how the targets move together. -Quantile forecasts do not carry this information. -A set of quantiles describes each target on its own and says nothing about the relationship between them. +Both the energy score and the variogram score need samples drawn from the joint predictive distribution. We cannot use the same multivariate scoring for quantile forecasts as they do not describe the relationship between targets in this way. It is tempting to line the quantile levels up across targets and treat them as samples (for instance by passing `sample_id = "quantile_level"` to `as_forecast_multivariate_sample()`). `scoringutils` will warn you about the leftover `quantile_level` column, but it will still return a score. However that score is not the one you want, as in this case, every quantile forecast ends up scored as though it had predicted perfect dependence (regardless of what the underlying model actually predicted). -It is tempting to line the quantile levels up across targets and treat them as samples, for instance by passing `sample_id = "quantile_level"` to `as_forecast_multivariate_sample()`. -`scoringutils` will warn you about the leftover `quantile_level` column, but it will still return a score. -That score is not the one you want. -Pairing the 5% quantile in one country with the 5% quantile in another assumes the two are perfectly rank correlated, so every quantile forecast ends up scored as though it had predicted perfect dependence, whatever it actually predicted. +## Comparing the energy and variogram scores -The two scores suffer differently. -The energy score still responds to the marginal distributions, so its values stay in a plausible range, but it loses almost all of its ability to separate forecasts with different correlation structures. -The variogram score fares worse, because dependence is the only thing it measures. -Replace the dependence with an assumption and what is left is close to a function of the observed values alone, which can rank models in the wrong order. +The energy and variogram scores are complementary to each other. The energy score summarises overall accuracy across all the pooled targets at once, a multivariate generalisation of the CRPS. The variogram score instead looks only at the differences between the specific targets that were pooled together to compare the size of the observed difference against the size of predicted differences. It is pairwise on whatever dimension `joint_across` specified. -If you have quantile forecasts and want to score them jointly, you have to supply the missing dependence structure. -You can ask forecasters for samples or trajectories instead of quantiles, which is the route several forecast hubs have taken. -Alternatively you can rebuild a joint distribution from the marginal quantiles using ensemble copula coupling (Schefzik et al., 2013) or the Schaake shuffle (Clark et al., 2004). -Both borrow a rank structure from elsewhere, usually a raw ensemble or the historical observations, so any score you compute afterwards partly reflects that choice. +For example, when pooling across forecast horizon, the variogram compares the differences between weeks and is therefore sensitive to the shape of the trajectory. But because it uses only differences between targets, it is insensitive to any bias that shifts the overall trajectory away from the observed data. -Multivariate point forecasts are a different matter. -A point forecast does specify a joint distribution, a degenerate one that places all its mass on a single vector, which is why the variogram score is well defined for it. +We can see this comparison between the two scores, comparing two models' forecasts of 1-3 week ahead cases in Germany as above. The dashed line shows the observed trajectory, with cases falling after the first week. -If, at any point, you want to score the same forecast using different groupings, you'd have create a new separate forecast object with a different grouping and score that new forecast object. - -## References +```{r} +#| fig.width: 7 +#| fig.height: 3.4 +#| fig.alt: >- +#| Median forecasts with 50% intervals for two models over three weekly +#| horizons in Germany, alongside the observed trajectory as a dashed line. +#| The ensemble stays roughly flat near the first observed value while EpiNow2 +#| rises away from it, and the observed cases fall steeply. +library(ggplot2) +library(data.table) + +cases_de <- as.data.table(example_sample_continuous)[ + target_type == "Cases" & + location == "DE" & + forecast_date == "2021-05-03" & + model %in% c("EuroCOVIDhub-ensemble", "epiforecasts-EpiNow2") +] -Clark, M., Gangopadhyay, S., Hay, L., Rajagopalan, B. and Wilby, R. (2004). -The Schaake Shuffle: A Method for Reconstructing Space-Time Variability in Forecasted Precipitation and Temperature Fields. -*Journal of Hydrometeorology*, 5(1), 243-262. - +band <- cases_de[, .( + median = median(predicted), + lower = quantile(predicted, 0.25), + upper = quantile(predicted, 0.75) +), by = .(model, target_end_date)] + +observed <- unique(cases_de[, .(target_end_date, observed)]) + +ggplot(band, aes(target_end_date, median, colour = model, fill = model)) + + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, colour = NA) + + geom_line() + + geom_line( + data = observed, aes(target_end_date, observed), + inherit.aes = FALSE, linetype = "dashed" + ) + + geom_point( + data = observed, aes(target_end_date, observed), + inherit.aes = FALSE + ) + + labs(x = "Target date", y = "Weekly cases", colour = "Model", fill = "Model") + + theme_scoringutils() + + theme(legend.position = "bottom") +``` -Gneiting, T., Stanberry, L. I., Grimit, E. P., Held, L. and Johnson, N. A. (2008). -Assessing probabilistic forecasts of multivariate quantities, with an application to ensemble predictions of surface winds. -*TEST*, 17, 211-235. - +We pool the three horizons into a single multivariate forecast for each model, then score both. -Schefzik, R., Thorarinsdottir, T. L. and Gneiting, T. (2013). -Uncertainty Quantification in Complex Simulation Models Using Ensemble Copula Coupling. -*Statistical Science*, 28(4), 616-640. - +```{r} +cases_mv <- as_forecast_multivariate_sample( + data = na.omit(cases_de), + joint_across = c("horizon", "target_end_date") +) +score(cases_mv)[, c("model", "energy_score", "variogram_score")] +``` -Scheuerer, M. and Hamill, T. M. (2015). -Variogram-Based Proper Scoring Rules for Probabilistic Forecasts of Multivariate Quantities. -*Monthly Weather Review*, 143(4), 1321-1334. - +The two scores disagree about which model performed better. The energy score prefers the ensemble, whose level stays close to the observed cases, while the variogram score prefers EpiNow2. The observed cases move sharply from week to week, and EpiNow2 predicts changes of a similar size, whereas the near-flat ensemble predicts changes that are too small. Note that EpiNow2 wins on the variogram score even though it has the wrong direction of the change. The energy score rewards getting the overall level right; the variogram score rewards getting the size of the joint movement right. From a9e9fa77131a6fb9d499664b999287de615229d1 Mon Sep 17 00:00:00 2001 From: kathsherratt Date: Fri, 17 Jul 2026 14:29:04 +0100 Subject: [PATCH 3/7] Update NEWS.md --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 20c31ed9a..414c575f9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,7 @@ - Fixed `summarise_scores()` producing a data.table with duplicate column names when the input `scores` object had no score columns (e.g. because every metric in `score()` warned and returned nothing). `summarise_scores()` now matches metric columns by exact name rather than regex partial match, and errors with a clear message when there is nothing to summarise (#1179). - Added internal S3 generic `get_forecast_type_ids()` so each forecast type declares the columns (beyond the forecast unit) that identify a unique row. `get_duplicate_forecasts()` now uses this instead of hard-coded column names (#888). - Removed the deprecated vignettes `Deprecated-functions` and `Deprecated-visualisations`. The code for removed functions (`plot_predictions()`, `make_NA()`, `plot_ranges()`, `plot_score_table()`, `merge_pred_and_obs()`) can still be found in the [git history](https://github.com/epiforecasts/scoringutils/tree/d0cd8e2/vignettes) (#1158). -- The vignette "Scoring multivariate forecasts" now explains why the energy score and the variogram score cannot be applied to quantile forecasts. Both scores require samples from the joint predictive distribution, which a set of marginal quantiles does not provide. Treating quantile levels as sample ids silently assumes perfect rank correlation between the targets being forecast jointly. +- Added a more descriptive explanation of the use of energy and variogram scores in the vignette "Scoring multivariate forecasts", including an extended description of use for pooling over single-origin forecast horizon and a multi-model comparison. Added a note explaining where these scores cannot be applied to quantile forecasts. # scoringutils 2.2.0 From 1f6d33e63839dccaebb91bba40e1fc340a693d1c Mon Sep 17 00:00:00 2001 From: kathsherratt Date: Sat, 18 Jul 2026 12:48:15 +0100 Subject: [PATCH 4/7] Ignore .quarto directory in vignettes Co-Authored-By: Claude Fable 5 --- vignettes/.gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/vignettes/.gitignore b/vignettes/.gitignore index 097b24163..9e2bd63c1 100644 --- a/vignettes/.gitignore +++ b/vignettes/.gitignore @@ -1,2 +1,4 @@ *.html *.R + +/.quarto/ From f0a9ebc3bc8ab91210623333eec7826da6e4226e Mon Sep 17 00:00:00 2001 From: kathsherratt Date: Sat, 18 Jul 2026 12:48:46 +0100 Subject: [PATCH 5/7] Link NEWS item to PR #1193 Co-Authored-By: Claude Fable 5 --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 414c575f9..62f5ad3f9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,7 @@ - Fixed `summarise_scores()` producing a data.table with duplicate column names when the input `scores` object had no score columns (e.g. because every metric in `score()` warned and returned nothing). `summarise_scores()` now matches metric columns by exact name rather than regex partial match, and errors with a clear message when there is nothing to summarise (#1179). - Added internal S3 generic `get_forecast_type_ids()` so each forecast type declares the columns (beyond the forecast unit) that identify a unique row. `get_duplicate_forecasts()` now uses this instead of hard-coded column names (#888). - Removed the deprecated vignettes `Deprecated-functions` and `Deprecated-visualisations`. The code for removed functions (`plot_predictions()`, `make_NA()`, `plot_ranges()`, `plot_score_table()`, `merge_pred_and_obs()`) can still be found in the [git history](https://github.com/epiforecasts/scoringutils/tree/d0cd8e2/vignettes) (#1158). -- Added a more descriptive explanation of the use of energy and variogram scores in the vignette "Scoring multivariate forecasts", including an extended description of use for pooling over single-origin forecast horizon and a multi-model comparison. Added a note explaining where these scores cannot be applied to quantile forecasts. +- Added a more descriptive explanation of the use of energy and variogram scores in the vignette "Scoring multivariate forecasts", including an extended description of use for pooling over single-origin forecast horizon and a multi-model comparison. Added a note explaining where these scores cannot be applied to quantile forecasts (#1193). # scoringutils 2.2.0 From 5987b5cbf0cb228252a0aa94e327ac6de34e496b Mon Sep 17 00:00:00 2001 From: kathsherratt Date: Sat, 18 Jul 2026 12:50:21 +0100 Subject: [PATCH 6/7] Revert "Ignore .quarto directory in vignettes" This reverts commit 1f6d33e63839dccaebb91bba40e1fc340a693d1c. --- vignettes/.gitignore | 2 -- 1 file changed, 2 deletions(-) diff --git a/vignettes/.gitignore b/vignettes/.gitignore index 9e2bd63c1..097b24163 100644 --- a/vignettes/.gitignore +++ b/vignettes/.gitignore @@ -1,4 +1,2 @@ *.html *.R - -/.quarto/ From 1d9122542d1ea7b7036f58cb936f687e270e0662 Mon Sep 17 00:00:00 2001 From: kathsherratt Date: Sat, 18 Jul 2026 12:59:28 +0100 Subject: [PATCH 7/7] Suppress namespace_linter for purrr::partial in vignette purrr is in Imports but not attached in the vignette, so the qualified call is required for the example to run. Co-Authored-By: Claude Fable 5 --- vignettes/scoring-multivariate-forecasts.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/scoring-multivariate-forecasts.Rmd b/vignettes/scoring-multivariate-forecasts.Rmd index 8ac92d00d..f7e52ca0d 100644 --- a/vignettes/scoring-multivariate-forecasts.Rmd +++ b/vignettes/scoring-multivariate-forecasts.Rmd @@ -102,7 +102,7 @@ score( example_multiv, metrics = list( energy_score = energy_score_multivariate, - variogram_score = purrr::partial( + variogram_score = purrr::partial( # nolint: namespace_linter. variogram_score_multivariate, p = 1 ) )