DiagrammeR::grViz('
digraph processes {
graph [rankdir = TB, fontname = "inherit", bgcolor = "white"]
node [shape = box, style = filled, fillcolor = "white",
color = "black", fontcolor = "black", fontname = "inherit"]
edge [color = "black"]
subgraph cluster_FG {
label = "Forecast-generating process"; color = "black"; fontcolor = "black"
F1 [label = "Forecaster participation"]
F2 [label = "Model specification"]
F3 [label = "Model implementation"]
}
subgraph cluster_TG {
label = "Target-generating process"; color = "black"; fontcolor = "black"
T1 [label = "Epidemic dynamics"]
T2 [label = "Observation process"]
T3 [label = "Target specification"]
}
F [label = "Forecast"]
T [label = "Target"]
P [label = "Predictive performance"]
F1 -> F; F2 -> F; F3 -> F
F -> P
T1 -> T; T2 -> T; T3 -> T
T -> P
}
')Using model-based evaluation to interpret variation in infectious disease forecast performance
Using model-based evaluation to interpret variation in infectious disease forecast performance
Katharine Sherratt (1), Rok Grah (2), Bastian Prasse (2), Friederike Becker (3), Jamie McLean (1), Sam Abbott (1), Sebastian Funk (1)
- Centre for Mathematical Modelling of Infectious Diseases, London School of Hygiene & Tropical Medicine
- European Centre for Disease Prevention and Control
- Institute of Statistics, Karlsruhe Institute of Technology
Email: katharine.sherratt@.ac.uk
Under development
Abstract
Forecasters predicting infectious disease outbreaks have met with varying success. Some of this variation in performance comes from the method used to make a forecast, when different models are better or worse at prediction. The rest comes from the target being forecast, when some outbreaks are easier or harder to predict than others. However, when many forecasters each predict many different targets, it becomes difficult to trace the impact of these factors shaping performance. Here we use a regression model to separate the effect of the forecasting method, from the difficulty of the target, in forecast performance.
We evaluated forecasts of weekly COVID-19 cases and deaths over two years across 32 European countries, scoring them against observed data with the Weighted Interval Score (WIS). We expected a model’s structure to shape how well it predicted, so we classified 47 models by structure (agent-based, mechanistic, semi-mechanistic, statistical, or human judgement) and estimated how much structure alone affected performance. A generalised additive mixed model let us adjust for everything that makes a target easier or harder to predict: the outcome being forecast, its level and trend, the dominant variant, the country, the forecast horizon, and differences between individual models.
Once we accounted for the difficulty of the target, no single type of model performed best. Differences in European COVID-19 forecast performance were driven more by which targets were hard to predict than by which modelling approach a forecaster used.
This approach sits between informal and fully formal ways of handling bias in evaluation studies. As infectious disease forecasting grows, we encourage evaluators to choose from a wider range of study designs, matching the formality of the method to the question, so they can isolate the part of performance they actually want to measure.
Main text
Background
Forecasters attempting to predict infectious disease outbreaks have met with varying success. This variation might come from many contributing factors. One such set of varying factors is the predictive ability of different forecasting methods. Another set of varying factors is the predictive difficulty of different forecast targets. With increasing numbers of both forecasters and forecast targets, we have more potential to learn about which factors can be modified to improve forecast performance, and to set appropriate expectations for performance between different forecast settings. This requires careful accounting among the many factors contributing to performance. Here, we suggest one approach to this using a joint modelling framework, and apply this to illustrate varying performance of forecasts of Covid-19 in Europe.
Some of the variation in performance can be attributed to differences between infectious disease target estimands. We see better predictive performance with stable epidemic dynamics [1] in larger populations [2], and when the target is observed with timely, accurate surveillance [3] [4], or has informative leading indicators [5]. Meanwhile, performance is consistently better at shorter than longer horizons [6], possibly based on the disease-specific serial interval that links disease to observation [7], prompting suggestions that there are inherent limits to the predictability of infectious diseases based on the complexity of the disease-surveillance-response system [8]. This variability makes it difficult to make fair comparisons across multiple forecast targets.
At the same time, variation in performance between forecasters arises from individual differences in the forecast estimation process: the structural approach and implementation choices made in developing a forecasting model. Forecasters may use many methods, ranging from purely data-driven statistical methods, to mechanistic models that embed explicit hypotheses about biological or epidemiological drivers of transmission. Different structural approaches may be better suited to different purposes, for example scenario projection, policy analysis, or real-time nowcasting. This choice demands a trade-off among complementary strengths and weaknesses of different methods [9–11], any of which may impact predictive accuracy. Previous work has compared predictive performance between these approaches, and shown little evidence for substantial differences [3,12,13].
This creates two sets of contributing factors to any single instance of realised performance Figure 1. Identifying these influences is important in order to fairly evaluate between different forecasters, and to set appropriate expectations for performance between different forecast settings.
To evaluate the contribution of different factors to predictive performance, we need a comparable sample of many independent forecasting models across many forecast targets [14,15]. This comparability is achieved in participatory forecast challenges, which collate multiple forecasters’ predictions for pre-specified targets in a standardised format [15]. This offers the opportunity for consistent evaluation between models, with the potential to explicitly compare the value of different modelling approaches [15]. However, such analysis is challenging, as it is prone to selection bias and confounding when evaluating among forecasters with differential participation across forecast targets. Aggregated comparisons may conflate differences in the predictive performance of each modelling approach, with differences in forecast target difficulty.
Some evaluators facing this challenge have used a range of study designs to account for this in evaluation. This typically includes controlling for factors in the data-generating process underlying the target. These study designs include: restricting the sample using exclusion criteria; stratification; data transformation; and assessing performance relative to a baseline (Table 1). Each of these study design strategies may usefully simplify some aspect of the problem. However, evaluators facing large, open forecasting challenges must handle a high degree of variation across both many forecast targets and many forecasting models simultaneously. The noisiness of these mixed signals means that any of these strategies for simplification may not be sufficient to isolate a single component of forecast performance and trace its effect.
| Study design | Example usage | Considerations |
|---|---|---|
| Restriction | Data revised retrospectively when forecasts contributed prospectively; exclude any target with e.g. >10% difference in retrospectively observed counts [7] | Simple and easy to communicate; loses power by excluding some targets altogether, and penalises forecasters who may appropriately forecast the revised data |
| Stratification | Present forecast scores stratified by factors of interest (e.g. model structure) or factors assumed to confound performance (e.g. forecast horizon) [3] | Classical observational study design; increasingly sparse sample with multiple interacting strata |
| Transformation | Normalise (e.g. per 100,000 population) and/or transform (e.g. log) each forecast-observation pair [16] | Maintains comparability across targets of different magnitudes, and the log transform may be closer to reflecting underlying epidemic growth rate; moving away from natural error scale requires clear communication |
| Matching | Post-process the set of forecast scores to create a pairwise comparison between forecasters. This relative score may be aggregated to include any forecaster against any other, or restricted to matching forecasters per-target [17]. The pairwise comparison may be against a ‘baseline forecaster’ available for all targets [18] | Able to include all forecast predictions; however the appropriate choice of comparator or baseline is unclear, and gives no interpretation of performance on an absolute scale (when all models perform well or poorly) |
For example, we may wish to assess the performance of different model designs used by a set of independent forecasters participating in any number of different forecast targets. When assessing across multiple forecast targets, we would need to separate out the model structure from the adaptability of different model structures to produce predictions across the different forecast targets. Faced with multiple forecasting targets, a forecaster might produce forecasts for any given target given only updated data, as is typical of statistical time-series and many semi-mechanistic models. An alternative strategy is to structure a model to mechanistically capture the data-generating processes of a specific forecast target, for example emphasising particular epidemic dynamics or exploiting local data sources. This suggests the hypothesis that design strategy has a separate influence on predictive performance, confounding an estimate of the effect of model structure alone on predictive performance.
Here, we develop a framework to identify the partial effect of model structure on predictive performance, while jointly adjusting for factors contributing to forecast target difficulty. We use proper scoring rules to evaluate short-term forecasts of COVID-19 across Europe. We analyse performance by model structure, using a generalised additive mixed model to also account for variation in the target epidemiological outcome, trend and level of incident outcomes, dominant variant phase, country, forecast horizon, and random error of each individual model. Rather than stratifying scores across these sparse, interacting subgroups, our hierarchical approach pools information across sparsely populated groups through partial shrinkage of the random effects, and represents continuous sources of difficulty such as incidence as smooth functions. We include a directed acyclic graph in the Supplement to represent the structure of our assumptions, but note that our work is exploratory rather than formally causal. We find that, despite superficial differences between model structures, most of the variation in observed performance can be explained by factors associated with the forecast target rather than the forecasting method. We conclude by noting the increasing breadth and application of infectious disease forecasting efforts, and call for a wider range of evaluative study designs able to accurately and robustly capture their value.
Methods
Study design
Forecast performance
We consider forecasts of incident COVID-19 cases and deaths between 8 March 2021 to 10 March 2023 for 32 European countries.
We used forecast data collected via the public European COVID-19 Forecast Hub, described in [19]. The platform solicited real-time forecasts for between one and four weeks ahead for 32 European countries. Any forecaster was eligible to participate in the Hub, and there were no selection criteria.
We excluded forecasts that did not report the full set of 23 quantiles, to ensure fair comparison among probabilistic results, and excluded baseline and ensemble models created by the Hub team. Full study eligibility and a STROBE flow diagram are detailed in the Supplement.
We evaluated forecasts against observed data collated by Johns Hopkins University (JHU). This included cumulative daily reported cases or deaths [20]. We calculated incident weekly counts, and used the dataset available as of 10 March 2023, when JHU stopped collating data.
Prior to evaluation, we normalised both forecasts and observations to incidence per 100,000 population. We then log-transformed both forecasts and observations. This transform leads to more comparable scores across locations and times, and can be interpreted as evaluating the ability of models to predict the growth rate [16,21]. To avoid zeros, we added a small amount (1e-7, or less than 1/1000th of the smallest positive value in the data) to all values.
We evaluated forecasts against observed data using the Weighted Interval Score (WIS) [17]. The WIS is a proper scoring rule and a generalisation of the continuous ranked probability score (CRPS). This assigns each probabilistic forecast a score that reflects the forecast’s dispersion, overprediction, and underprediction.
Forecasting model structure
We designed a classification scheme for participating forecasters’ model structures. This included: “Agent-based” (with detailed individual-level interactions), “Mechanistic” (broadly based on compartmental models), “Semi-mechanistic” (statistical models with epidemiological elements, e.g. delay distributions and/or generation intervals), “Statistical” (context-agnostic time-series models), “Judgement” (ensembles of human-judgement based predictions).
We classified each model using self-reported metadata: an optional description of methods, and links to model code, website or publication (schema and example in the Supplement). Four members of the research team (KS, SF, FB, JM) classified models independently, with at least two rating each model. The majority vote determined the final classification, with ties resolved by an additional independent rater and all classifications reviewed by two investigators (SF, KS) for final revisions (see Supplement).
Forecast target difficulty
We included the epidemiological target (incident cases or deaths) as a fixed factor within a single model, rather than stratifying the analysis. With only two levels, a fixed factor gives one fully-identified, interpretable contrast (deaths versus cases on the log scale), whereas a random effect would require more levels to identify its variance. We then included a variety of covariates to account for forecast task difficulty.
We first included forecast horizon to account for the increasing difficulty of forecasting further into the future. We then accounted for variation in epidemic transmission dynamics underlying each forecast target. As the transmission state is unobserved, we adjusted for observed factors associated with it. Factors upstream of transmission included: country location, as a proxy for dynamics e.g. national policy or population structure influencing transmission; and dominant variant phase. We used publicly available sequence data to classify each target location-week according to the dominant circulating SARS-CoV-2 variant, producing a discrete sequence of phases (see Supplement) [22]. Factors downstream of transmission included: count of incident outcomes in the target week; and trend of incidence: “Stable”, “Decreasing”, or “Increasing”, based on the difference over a three-week moving average of incidence (with a change of +/-5% as “Stable”; see Supplement).
Task participation and residual confounding
We classified each model as participating in forecasting either a single country location or more than one location. We used a binary classification rather than a count to address data sparsity (few teams submitted for between 2 and 32 countries), and because forecasters could change the number of targets they submitted over time.
We additionally included a random effect for individual model identity, to separate the between-method-class contribution from variation between individual model implementations. We did not attempt to adjust for unobserved characteristics of the forecasting team, such as capacity, infrastructure, or experience, although these plausibly shape both model structure choices and realised forecast performance.
Analysis
We used a single hierarchical model structure fit across both epidemiological targets, with target included as a fixed factor. We treated the WIS of each forecast \(i\) as Gaussian with a log link, so that the linear predictor \(\eta_i\) acts multiplicatively on expected score:
\[ \mathrm{WIS}_i \sim \mathcal{N}(\mu_i,\, \sigma^2), \qquad \log \mu_i = \eta_i. \]
The linear predictor decomposed into a fixed intercept, a fixed effect for the epidemiological target, a sum of random-effect contributions from categorical covariates, and two smooth terms for continuous covariates:
\[ \eta_i = \beta_0 \;+\; \beta_{\text{target}}\,\mathbb{1}[\text{Deaths}_i] \;+\; \underbrace{\sum_{g \in \mathcal{G}} u^{g}_{j_g(i)}}_{\text{random effects}} \;+\; \underbrace{f_{\text{Inc}}(\mathrm{Incidence}_i) \;+\; f^{\text{Horizon}}_{k(i)}(\mathrm{Horizon}_i)}_{\text{smooths}}. \]
Here \(\beta_{\text{target}}\) is a fixed parametric coefficient (cases as the reference level), giving the deaths-versus-cases contrast on the log scale; it is not part of the zero-mean Gaussian random-effects set \(\mathcal{G}\).
Here \(j_g(i)\) denotes the level of grouping factor \(g\) observed for forecast \(i\). We assume random effects are drawn independently from a zero-mean Gaussian with its own variance:
\[ u^{g}_{j} \sim \mathcal{N}(0,\, \sigma^2_g), \qquad g \in \mathcal{G} = \{\text{Method},\, \text{CountryTargets},\, \text{Trend},\, \text{Variant},\, \text{Location},\, \text{Model}\}. \]
The two smooths captured continuous covariates. \(f_{\text{Inc}}\) was a penalised thin-plate regression spline over (log-transformed) incidence, and \(\{f^{\text{Horizon}}_k\}\) was a factor-smooth interaction (sum-to-zero constrained, with \(k=3\) basis functions per forecasting model), allowing the horizon effect to vary by model.
We used summary statistics to describe WIS by model structure and target strategy. We fitted all models using the mgcv package in R [23], using fast restricted maximum likelihood estimation (fREML). We checked model fit with visual inspection of residuals via Q-Q plots. To assess confounding, we fitted univariate models to contrast unadjusted and adjusted estimates of the partial effect of model structure on WIS.
All data and code are publicly available online [24].
Results
We evaluated the predictive performance of 47 forecasting models produced by 37 separate modelling teams participating in the European COVID-19 Forecast Hub. We identified both differential participation and varying performance among models, using the weighted interval score (WIS) as a standardised measure of performance. We then used a generalised additive mixed model (GAMM) to associate predictive performance with model structure. We accounted for confounding from forecasters’ selection of geographic targets (the number of countries a model targeted), and forecast target difficulty (factors associated with the epidemiological context).
Forecaster characteristics
| Models (%) | Models (%) | Single-country (%) | |
|---|---|---|---|
| Overall | 42 (100%) | 38 (100%) | NA |
| Method | |||
| Semi-mechanistic | 9 (21.4%) | 10 (26.3%) | 2/12 (17%) |
| Statistical | 11 (26.2%) | 7 (18.4%) | 4/12 (33%) |
| Mechanistic | 16 (38.1%) | 16 (42.1%) | 10/17 (59%) |
| Agent-based | 3 (7.1%) | 2 (5.3%) | 3/3 (100%) |
| Judgement | 3 (7.1%) | 3 (7.9%) | 0/3 (0%) |
| Geographic scope | |||
| Single-country | 19 (45.2%) | 14 (36.8%) | NA |
| Multi-country | 23 (54.8%) | 24 (63.2%) | NA |
We summarised 47 model characteristics by model structure and geographic target selection (Table 1). We included all participant forecasts for any target; full eligibility criteria and exclusion counts are shown in the STROBE flow diagram (Supplementary Figure S1). Over the study period, forecasters could participate in any of 26,624 available forecast targets (comprising 104 weeks, 32 countries, 4 horizons, and 2 epidemiological outcomes). Median forecaster participation was 3% (IQR: 2-12%, ranging from 0 to 96) across available targets, reflecting substantial variation in participation across countries and time.
We classified 17 models as using mechanistic structure, 12 as statistical, 12 as semi-mechanistic, 3 as agent-based and 3 that used human judgement forecasting (Table 1; Supplementary Table S2). Raters disagreed on 17 (36%) model classifications. The majority of 2/3 was used as the final classification, with additional manual review which in all cases retained the majority decision. Raters most often disagreed when classifying semi-mechanistic models, with 8 out of 12 receiving one or more votes for a statistical model structure. Under the final classification, approximately one-third of all forecast predictions were from mechanistic, semi-mechanistic, and statistical models each. Agent-based and judgement models provided fewer forecasts, representing <2% of forecasts.
We classified 19 models as forecasting for a single country and 28 as forecasting for multiple countries; the countries targeted are detailed in the Supplement. Among the 28 multi-country models, only 2 consistently forecast for the same number of countries throughout the study period.
Forecast performance
We descriptively summarised forecast performance across the study sample. Unadjusted performance varied across the entire cohort of models over space and time. We saw consistently worsening performance during summer 2021, coinciding with emergence of the Delta variant. By model structure, we noted apparent improved performance from agent-based and judgement models, compared to mechanistic and statistical approaches; and single-country models appeared to outperform multi-country models, particularly for death outcomes from spring 2022. However, model structure differed substantially between geographic targets and across features of the epidemiological target, making crude unadjusted comparisons unreliable.
Effect of model structure
We fitted a generalised additive mixed model to give adjusted estimates of the partial effect of model structure, while controlling for varying forecaster target selection and epidemic dynamics between targets. In this structure, partial effects are deviations from the grand mean WIS under a sum-to-zero constraint, so a negative value indicates better-than-average performance. The WIS was highly right-skewed, and we used a log link to account for this (diagnostics in Supplementary Figures S2-S3). The log link models the score multiplicatively but leaves some skew in the residuals. A sensitivity analysis modelling the log-transformed score directly substantially improves the residual distribution while preserving the direction of all effects and leaving the model-structure conclusions unchanged (Supplement). Exponentiating a partial effect gives a multiplicative ratio relative to the grand-mean WIS, where 1 indicates average performance (e.g. a partial effect of −0.1 corresponds to a ratio of 0.9, a WIS approximately 10% lower than average). We report the exponentiated ratio (with 95% confidence intervals) in the main text (Table 2), with the raw partial effects on the log scale in the Supplement.
After adjustment forfeatures of the forecast target , no single structural approach dominated. Adjusted point estimates clustered around the grand mean, and confidence intervals overlapped throughout (Table 2). The largest shifts were among agent-based and human judgement models, which appeared better than average in unadjusted estimates but showed no difference from other model structures after adjustment (adjusted ratios within 1% of the grand mean). We noted that adjustment consistently shrank the standard error compared to univariate estimates, narrowing the intervals around these overlapping estimates (Figure 3).
| Variable | Group | Unadjusted ratio (95% CI) | Adjusted ratio (95% CI) |
|---|---|---|---|
| Method | Agent-based | 0.87 (0.71, 1.06) | 1 (1, 1) |
| Judgement | 0.76 (0.63, 0.92) | 1 (1, 1) | |
| Mechanistic | 1.08 (0.9, 1.29) | 1 (1, 1) | |
| Semi-mechanistic | 1.17 (0.97, 1.4) | 1 (1, 1) | |
| Statistical | 1.21 (1.01, 1.45) | 1 (1, 1) |
Other drivers of forecast performance
Among our selection of covariates, we found performance was more driven by the target epidemiological features than by forecaster model features (with full reporting in Supplement). These features were modelled jointly across both epidemiological targets, so each estimate is a single effect shared by case and death forecasts. Holding all other features constant, the most predictable forecast targets appeared during periods of stability in the epidemic curve, with better predictive performance relative to periods of increasing or decreasing trends (ratio 0.5 (0.25-1)). Forecasting was more difficult during increasing trends (ratio 1.57 (0.78-3.13)), compared to stable or decreasing periods.
Predictability also varied between phases of variant dominance. We estimated better predictive performance during the earlier phases dominated by the Alpha and Delta variants (Alpha ratio: 0.73 (0.61-0.87); Delta ratio: 0.83 (0.69-0.99)). We estimated worse performance after Omicron variants became dominant (with worse performance in the first Omicron BA.1 variant, ratio: 1.32 (1.1-1.59)). Death forecasts achieved a lower WIS on the log-incidence scale than case forecasts, as a fixed contrast holding all other features constant (ratio: 0.27 (0.26-0.27)).
We identified residual confounding, with substantial unexplained variation between individual model performance. Figure 4 shows the adjusted partial effect for each individual model. This can be interpreted as variation in model-specific performance after accounting for all other factors considered here. This indicated that the variables we selected to account for variation between forecasters (i.e. model structure; one or multiple geographic targets) were not sufficient to fully explain differences in performance.
Sensitivity analysis
We repeated this analysis to look at performance on the natural scale (rather than log-transformed predictions and observations), presented in the Supplement. In developing this work, we considered several alternative model specifications. We selected an initial set of variables a priori, and changed this selection after updating our view of confounding factors, for example including an effect for dominant variant phase. We indicate these choices in the Supplement. Regardless of a priori variable selection, we did not find any meaningful change in partial effect sizes among alternate model specifications. Refitting with an identity link in place of the log link also left all partial effects unchanged to two decimal places.
Discussion
In this work, we aimed to assess the partial effect of model structure on the accuracy of forecasts of Covid-19 across European countries. We developed a model-based approach to adjust for the key sources of variation between different forecast targets. After adjustment, no single structural approach dominated. Adjusted estimates for each model structure were no different from the overall average, though substantial variation in individual model performance remained unexplained. We suggest that when evaluating across multiple forecast targets, it is necessary to account for the factors contributing to variation between targets, before considering factors within the forecasting process.
In this work we had limited ability to demonstrate the model-based approach used here, with a small and biased sample reducing any power to detect true differences in performance between model structures. While we used a large dataset of predictions over time, we had a much smaller effective sample size of 47 independent models. We relied on an opportunistic cohort of models from voluntary contributions to the European Hub. For example, we only observed 3 agent-based models, where all three represented here targeted a single country. This meant our results may be biased by unobserved characteristics that differentially affected participation, such as the resource intensity of model development, or the capacity to incorporate domain expertise. These factors may drive both the choice of model structure and be more useful as targets for identifying potential for modifying forecaster performance. In addition, we may have misclassified the exposure, with very little metadata to classify model structures. Independent raters disagreed on the classification of most semi-mechanistic models, and we did not track changes to each model’s individual methods over time. A more systematic sample of models would support a deeper investigation of factors driving model performance. Future work could approach this by combining data from across multiple similar Hub projects and a more detailed assessment of the characteristics of the model cohort.
Our findings are consistent with the hypothesis that variation in forecast performance is driven more by differences among features of the forecast target, than by differences among forecasting strategies. Previous work has also failed to differentiate performance among model structures, across a variety of outbreaks and epidemic dynamics [3,12,13]. Some prior moderate evidence suggests that statistical methods gave better performance than more mechanistic approaches [5,25,26]. On the other hand, bespoke mechanistic models are more easily tailored to a single target which may give a specific performance advantage for such targets [27,28]. A study of COVID-19 forecasts in Germany and Poland found the two best-performing models were those agent-based models specifically tuned to the national context [1]. In this study, we noted that the more mechanistic models produced a wide range of predictive performance. This may reflect a greater range of flexibility to incorporate forecaster experience [25], and understanding of the data generating processes that lead the forecast target [29]: for example, country-specific policies influencing local epidemic dynamics, or the specific way in which the public data reflected underlying transmission processes. While not considered here, we note the widely replicated finding that ensemble combinations across multiple models often produces the most reliable predictive performance (e.g. [6,7,30]).
Greater investment [31] and standardisation [32,33] in infectious disease forecasting has created a much broader comparable sample of forecast performance. While this may reduce the likelihood of random error, it cannot alone account for bias in comparisons of performance. Our analysis represents a middle ground in formally considering bias when designing evaluation studies. On one hand, this degree of formality may not be necessary when there are minimal sources of confounding (e.g. a single forecast target), or if the effective sample size is large enough to support a simpler stratification. In our setting, however, cross-classifying forecasts by the categorical covariates alone (epidemiological target, horizon, location, trend, variant phase, geographic scope, and model structure) yields over 46,000 strata; the resulting cells are sparse and highly uneven, and incidence, a continuous covariate, cannot be stratified without arbitrary binning. On the other hand, a more formal analysis could fully specify the causal estimand with a directed acyclic graph; use this to derive a minimal sufficient adjustment set of covariates, and report a quantitative bias analysis (such as an E-value) for unmeasured confounding. Further model comparison projects could enable this by encouraging methodological diversity among participants, and ensuring that detailed structured information on model metadata (such as methodology and model revisions) is collated alongside numerical predictions [14]. We would particularly encourage studies that combine information from multiple collaborative efforts to increase sample size of independent models. This could enable assessing more detailed components of the forecasting process, such as the amount of local context that modellers use in creating or adapting their forecast model, and similarly whether performance scales with efforts to revise, incorporate domain knowledge, or evaluate model performance.
Comparisons like the one we have conducted here support a more informed view of the role of modelling in pandemic preparedness and response. By accounting for forecast target difficulty, we adjusted for a lower bound of error set by the forecast target itself, and estimated the realistic room for improvement from improving forecasting methods. Our specific results suggest that focusing on improving understanding of the forecast target may yield more gains in forecast performance than advancing specific modelling techniques. As the field of infectious disease forecasting develops, we encourage evaluators to consider a wider range of study designs, adopting informal to formal methods where appropriate in order to isolate specific components of forecast performance.
References
Supporting Information
Supplementary Figure 1. Eligibility criteria for models contributing case and death forecasts to the European COVID-19 Forecast Hub, March 2021 - March 2023
Supplementary Table 1. Model characteristics contributing to the European COVID-19 Forecast Hub, by method used, number of countries targeted, and number of forecasts contributed
Supplementary Figure 2. Epidemic trends (cases)
Supplementary Figure 3. Epidemic trends (deaths)
Supplementary Figure 4. Model diagnostics (cases)
Supplementary Figure 5. Model diagnostics (deaths)
Data availability
The codebase for this paper is publicly available at the Github repository “epiforecasts/eval-by-method” and Zenodo with DOI: 10.5281/zenodo.14903162. Forecast and observed data were sourced from the European COVID-19 Forecast Hub, available to view online. All Hub data are now archived at the Github repository: “european-modelling-hubs/covid19-forecast-hub-europe_archive” and Zenodo with DOI: 10.5281/zenodo.13986751. Specific data used in this work are available in the Github repository for this paper.