Summary of included team forecasts

A total of 14,960,171 forecasts were submitted by 67 teams throughout the analysis period (see S1 Appendix for submission patterns over time). Because forecasts were submitted at multiple geographic scales, we stratified analyses for 1) national forecasts, 2) state (including all 50 states), territory (US Virgin Islands and Puerto Rico), and DC forecasts), 3) county-level forecasts (including all 3,143 counties and county equivalents), split into five equal sized groups based on county population size.

We first assessed the locations, horizons, and time periods forecasted by each team to ensure that forecasts included all required quantiles and horizons and to limit comparisons to teams with substantial overlapping spatiotemporal coverage. Briefly, teams were included if they submitted the full range of required quantiles, included at least 50 of states/territories/DC or 75% of counties, and produced forecasts at least four weeks into the future for at least 50% of the time points in the study period. At the national level, 22 sets of team forecasts met these criteria (5,136 forecasts across dates and forecast horizons), 23 sets of team forecasts met the state/territory level criteria (280,132 forecasts across jurisdictions, dates, and forecast horizons), and 15 sets of team forecasts met the county-level criteria (9,415,460 forecasts across counties, dates, and forecast horizons). Overall, 64.8% of all submitted forecasts were included in the analysis (9,700,728 forecasts). Of the included forecasts, 11 sets of team forecasts met the inclusion criteria for analyzing submissions across all geospatial scales (8,125,220 forecasts for specific locations, dates and forecast horizons).

Each team included in the analysis submitted forecasts that were generated from unique model structures, data inputs, and assumptions (S1 Appendix). Two naïve models (the COVIDhub-baseline and CEID-Walk) and four ensemble models (the COVIDhub-4_week_ensemble, the COVIDhub-trained_ensemble, LNQ-ens1, and UVA-Ensemble), which combined multiple forecasts into one, were included in the 26 models evaluated (see Table A in S1 Appendix). The COVIDhub-baseline model projects the number of reported cases in the most recent week as the median predicted value for the next 4 weeks. CEID-Walk is a random walk model with a simple method for removing outliers. A total of seven models included data on COVID-19 hospitalizations, 12 models incorporated demographic data, and seven models used mobility data. Of the 26 evaluated models, three (COVIDAnalytics-DELPHI, CU-select, and UCLA-SuEIR) assumed that social distancing and other behavioral patterns changed during the prediction period.

The evaluation period consisted of 1–4 week ahead forecasts submitted in the 73 weeks from July 28, 2020, through December 21, 2021. Multiple phases of the US epidemic were included: the late summer 2020 increase in several locations, a large late-fall/early-winter surge in 2020/2021, the rise and fall of the Delta variant in the summer and fall of 2021, and the early phase of the Omicron variant’s dominance in winter 2021 (Fig 1A). Performance of the national ensemble forecasts varied over this period (Fig 1B). For some forecasts, the median predictions were close to the cases eventually reported, and most reported numbers fell within the 95% PIs. However, forecasts made at other times, such as January 2021 or December 2021, diverged widely from the reported data. At those times, the forecasts missed substantial decreases and increases, respectively, with reported cases falling within the 95% prediction interval for only 1-week ahead forecasts.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

Weekly incident reported COVID-19 cases per 100K population, nationally (in black) and per state/territory/DC (in gray), over time in panel A. Panel B shows a subset of COVIDhub-4_week_ensemble forecasts (in green) over time, with the median predictions represented as lines and points and the 95% prediction intervals in bands. Reported incident cases (counts per week) are shown in gray. In both plots, the black, dashed vertical line shows the date that public communication of the case forecasts was paused.

https://doi.org/10.1371/journal.pcbi.1011200.g001

Aggregate performance

We evaluated aggregate forecast performance with two metrics: Weighted Interval Score (WIS), a proper score considering both precision and accuracy, and prediction interval coverage, an indicator of forecast uncertainty. Lower WIS values reflect forecasts with probability mass closer to observed values. We assessed scaled pairwise WIS relative to the baseline model (referred to throughout as relative WIS, or rWIS) for national and state/territory/DC forecasts (Fig 2). A rWIS less than one indicates performance that is better than the baseline model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Percent of weeks with complete submissions for all sets of team forecasts, scaled, pairwise relative Weighted Interval Score (rWIS; see Methods for description), observed 95% prediction interval coverage, by geographical scale of submitted forecasts.

Teams are sorted by increasing state/territory/DC rWIS values.

https://doi.org/10.1371/journal.pcbi.1011200.g002

Overall, seven of 22 team’s forecast models outperformed the COVIDhub-baseline model at the state/territory/DC level (i.e., had rWIS values less than 1.0), and 11 outperformed the baseline model at the national level. Six of these teams outperformed the baseline model at both scales: LNQ-ens1, COVIDhub-4_week_ensemble, USC-SI_kJalpha, LANL-GrowthRate, Microsoft-DeepSTIA, and CU-select.

Prediction interval coverage at the 95% level should be close to 95% for well calibrated forecasts. However, it was lower for most sets of team forecasts, with only one (LNQ-ens1) having coverage of at least 90% at all scales, while others were as low as 23%. PI coverage at 50% and 80% levels were also well below nominal levels for most sets of team forecasts, including the COVIDhub-4_week_ensemble (Fig 3). For the 50% prediction interval, no sets of team forecasts had coverage better than 36% at any scale. Only two sets of team forecasts had better coverage than the COVIDhub-4_week_ensemble for the geographic scales in which they submitted forecasts: LNQ-ens1 (all scales) and JHU_UNC_GAS-StatMechPool (state/territory/DC and large county levels). We also found that calibration, WIS, and prediction interval coverage were all worse at 4-week horizons compared to 1-week horizons (S6 Appendix).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

Expected and observed coverage rates for central 50%, 80% and 95% prediction intervals aggregated over time and horizon for national forecasts (panel A), state/territory/DC forecasts (panel B), and the largest county forecasts (panel C). The dashed line represents optimal expected coverage. Team forecasts that had closer to nominal coverage than the COVIDhub-4_week_ensemble model at all three coverage levels are labeled on the right side of the plots.

https://doi.org/10.1371/journal.pcbi.1011200.g003

Forecast skill also showed distinct patterns across jurisdictional scales, with rWIS decreasing for larger jurisdiction scales (e.g., national vs. state/territory) or population sizes (e.g., larger counties vs. smaller counties, Fig 4) for most sets of team forecasts. In contrast to this general trend, for three sets of team forecasts, that pattern was inverted, one team had no distinct pattern, and the COVIDhub-4_week_ensemble had markedly consistent rWIS across all scales. Consistent with the aggregate findings, both LNQ-ens1 and COVIDhub-4_week_ensemble had rWIS lower than 1.0 at all scales, while LANL-GrowthRate had rWIS greater than 1.0 for smaller counties.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Scaled, pairwise relative Weighted Interval Score (rWIS) (see Methods for description) by spatial scale for sets of team forecasts that submitted forecasts for the US nation, states/territories/DC, and all US counties.

WIS is averaged across all horizons. The COVIDhub-baseline model has, by definition, a rWIS of 1 (horizontal dashed line). Teams are ordered by increasing state/territory/DC rWIS with the most accurate model on the left. Points for each team are staggered horizontally to show overlapping WIS values.

https://doi.org/10.1371/journal.pcbi.1011200.g004

Performance across jurisdictions

There was additional variability in forecast skill between jurisdictions. Only two team forecasts (LNQ-ens1 and COVIDhub-4_week_ensemble) performed as well as or better than the baseline for all included states and territories (Fig 5). Variation was higher between team forecasts than between specific jurisdictions, but the baseline model tended to outperform more models in some jurisdictions (e.g., the baseline was better in Colorado, Kansas, Puerto Rico) than in others (e.g., the baseline was worse in Mississippi, South Carolina, West Virginia). Spatial correlation is intrinsic to COVID-19 spread making correlation in the forecasts likely regardless of forecast skill. We found some evidence of spatial correlation in rWIS for many team forecasts (ensemble Moran’s I: 0.36, p-value = 0.001, Fig A in S4 Appendix).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Scaled, pairwise relative Weighted Interval Score (rWIS; see Methods for description) by location for national and state/territory/DC forecasts, averaged across all horizons through the entire analysis period.

National estimates are displayed first, followed by jurisdictions in alphabetical order. Team forecasts are ordered by increasing average state/territory/DC rWIS.

https://doi.org/10.1371/journal.pcbi.1011200.g005

Performance over time

WIS also varied over time (Fig 6). For example, all models had relatively high WIS in December 2020-January 2021 and low WIS in June 2021. Prediction interval coverage also varied between teams and over time, with most team forecasts exhibiting times of low coverage. Across most time points, the baseline model outperformed many team forecasts, including the COVIDhub-4_week_ensemble, though the ensemble more often outperformed the baseline in both WIS and prediction interval coverage at the national, state/territory, and large county scales. Increased WIS and decreased prediction interval coverage generally occurred with increasing case counts, such as in the fall of 2020 and summer of 2021. The worst performance was in the early Omicron wave in the winter of 2021. For the last set of ensemble forecasts posted by the CDC in December 2021 (https://www.cdc.gov/coronavirus/2019-ncov/science/forecasting/forecasts-cases.html), the WIS reached the highest level ever for all scales and the reported case numbers were outside the 95% prediction interval for most locations at every forecast horizon.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Forecast accuracy over time, aggregated by geographic units, forecast horizon, and prediction date.

Panels A-C show average Weighted Interval Score (WIS); panels D-F show 95% prediction interval coverage. The black, dashed vertical line in all panels shows the date that public communication of the case forecasts was paused. The black, dashed horizontal line in panels D-F shows nominal 95% prediction interval coverage. National level forecasts are presented in A and D, state/territory/DC forecasts in B and E and large county forecasts in C and F.

https://doi.org/10.1371/journal.pcbi.1011200.g006

To further investigate these temporal patterns in performance, we first classified each forecasted week as increasing, peak, decreasing, or nadir based on the estimated time-varying reproduction number for that given week and jurisdiction (Fig A in S5 Appendix). We then fitted Gaussian generalized estimating equations (GEE) models for each set of team forecasts, using a normalized, log transformed WIS value per forecast time and location as the model outcome. The regression models were adjusted for each prediction horizon and included a natural spline with two degrees of freedom for the time/state reported case counts to adjust for intrinsic increases in WIS due to higher values in reported cases (see S7 Appendix). In agreement with the aggregated results (Fig 2), we found that the expected WIS at the mean number of case counts across all jurisdictions was lower than the baseline for the better performing models (6 team forecasts and the ensemble) and higher than the baseline for others (8 team forecasts).

Forecast skill and coverage also varied across epidemic phases (Fig 7B and 7C). Compared to the baseline model across all phases, overall WIS for most models was better in nadir and peak phases and worse in increasing and decreasing phases. Likewise, 95% prediction interval coverage was highest in the nadir phase for nearly all teams while coverage in other phases was mostly lower than 95 percent. LNQ-ens1 and the COVIDhub ensemble had better WIS values than the baseline model in all epidemic phases between August 1, 2020, and January 15, 2022, with LNQ-ens1 also exhibiting close to nominal coverage across all phases.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7.

Estimated marginal mean Weighted Interval Score (WIS) and 95% confidence intervals for mean cases from team-specific GEE models for all 51 jurisdictions (Panel A). The 95% confidence intervals for the COVIDhub-baseline model are shown in dashed red vertical lines. Panel B presents each team’s estimated marginal mean WIS per phase, scaled to the COVIDhub-baseline model’s estimated marginal mean WIS for all epidemic phases. Teams with higher estimated marginal mean WIS values (i.e., greater than 1.0) are presented in shades of orange while teams with lower estimated marginal mean WIS (i.e., less than 1.0) are shown in shades of green. Forecasts for a team in a particular phase are marked with an asterisk (*) if the 80% confidence interval of the expected WIS outcome (normalized and on the log scale) was estimated by a model to be lower than the average expected WIS of the COVIDhub-baseline model across all phases. Panel C shows each team’s mean 95% prediction interval coverage in each epidemic phase.

https://doi.org/10.1371/journal.pcbi.1011200.g007

We classified each forecast as increasing, decreasing, or stable or uncertain based on the 50% prediction interval relative to the most recent observed value. If the forecasts were able to correctly predict the direction of the epidemic phase, we would expect a high percentage of forecasts with an increasing trajectory to occur in an increasing epidemic phase, and likewise, a high percent of stable/uncertain trajectory forecasts in the peak and nadir phases and a high percent of forecasts with a decreasing trajectory in the decreasing epidemic phase. Increasing forecasts were most commonly made for the peak phase (predicting continued increases rather than a peak) followed by the increasing phase (Fig 8A), though overall the most common forecasts were for stable or uncertainty trajectories (Fig 8B). While the percent of decreasing forecasts for the decreasing epidemic phase was high, forecasts often predicted declines for the nadir as well and some teams also frequently predicted decreases for the increasing phase (Fig 8C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8.

Percent of forecasts with predicted increasing trajectory per epidemic phase (Panel A), predicted stable or uncertain trajectory per epidemic phase (Panel B), and predicted decreasing trajectory per epidemic phase (Panel C). In each plot, epidemic phase labels are in bold when they correspond with the predicted direction of the forecast.

https://doi.org/10.1371/journal.pcbi.1011200.g008

To examine whether our results were affected by reporting anomalies, we also conducted sensitivity analyses for data revisions, when data were revised at a later date, and for outlier data points, when reported cases were outside of weekly expected ranges (see S2 Appendix). We first identified weeks in which revised case counts as of April 2, 2022, differed from the case counts initially reported for that week, excluded them from the dataset, and reran the GEE models. With this partial dataset, the results were essentially unchanged. Next, we identified outliers as reported case counts outside of the expected range by at least two of the three following algorithms: a rolling median, a seasonal trend decomposition, and a seasonal trend decomposition without a seasonality term, each method over a 21-day window. Approximately 3% of weeks (686 of 27,489 total week-location combinations in the analysis period) had at least one day of reported cases identified as an outlier. We then excluded the weeks with outliers and the week following an outlier and reran the GEE models. This sensitivity analysis had comparable results to the models with the full data (see Fig C in S2 Appendix, Panel 1).

Inclusion criteria

Teams were included in the evaluation when they submitted forecasts with a complete set of quantiles for each 1- through 4-week ahead target predictions. Additionally, teams must have met the following inclusion criteria:

1. had predictions for at least 50 locations (states, territories, or DC) for the state, territory, and DC level analyses; and for at least 75% of counties included in each population size quintile per submission week for the county-level analyses;
2. had submissions for at least 50% of the weeks included in the analysis period per location forecasted.

Application of these inclusion criteria provided a more comparable set of forecasts for scoring with the attempt to reduce biased scores if teams only forecasted for a limited number of locations or number weeks. Teams meeting these inclusion criteria, and their submissions over time, are depicted in Figs A and B in S1 Appendix.

Ground truth

Forecasts were evaluated against the reported COVID-19 case reports collated by the Johns Hopkins Center for Systems Science and Engineering (CSSE) [48]. To calculate weekly incident reported cases, we subtracted the cumulative count for each Saturday from the cumulative count for the next Saturday, such that each incident weekly count reflects the number of cases reported from Sunday through Saturday in a given week. We aggregated reported counts from smaller geographic units into their larger unit. For example, counts in a given state are the aggregate of the county-level reported counts and national counts are the sum of all states, territories, and DC.

CSSE reports data in real-time. Thus, data may be revised if the reporting health system makes public updates to its surveillance data. At times, such revisions may result in negative daily counts or in increases to case counts if the date of cases is shifted from one day to another or the definition of a reportable case is changed. We examined the percent change between the first reported cases in each state, DC, and territory per date relative to the counts in the surveillance file from April 2, 2022. We also assessed the influence of revised data on the final model outcomes (see S2 Appendix) and the presence of negative case counts in the timeseries. Less than 1 percent of time points in the analysis period had negative daily case counts in the largest US counties. Negative counts were observed at the state/territory level only twice: in Missouri during the week of April 17, 2021, and Virgin Islands during the week ending October 10, 2020. The state of Florida reported 0 cases on November 27, 2021. We excluded all weeks and locations with negative counts as well as the week with 0 incidence in Florida in our primary analyses.

Additionally, we also examined whether a reported case count was an outlier in the case trend for each state. Anomalies in case data trends have not been uncommon throughout the pandemic, as reporting entities have uploaded large batches of surveillance data on a single day. To assess whether cases were outside of the expected range of reported cases over time, we applied three outlier detection algorithms, each with a 21-day window: a rolling median, a seasonal trend decomposition, and a seasonal trend decomposition without a seasonality term. We then classified a given count as an outlier if it was detected as such by at least two of the three algorithms. Using these data, we ran several sensitivity analyses to assess the likely impact of anomalous data points on model performance. Sensitivity analyses examining the robustness of our findings to reporting anomalies are presented in S2 Appendix.

Additional information about the CSSE data, and revisions to the dataset, is publicly available on a GitHub repository:

https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data.

Forecast locations

Forecasts for incident cases were submitted for the national level, 50 states, 5 territories (American Samoa, Guam, the Northern Mariana Islands, Puerto Rico, and the US Virgin Islands), the DC, and 3,142 US counties. We excluded two relatively new (2019) counties in Alaska (Federal Information Processing Standard code 02063 and 02066) because they were not included by most teams. Because fewer teams submitted forecasts for American Samoa, Guam, the Northern Mariana Islands, we excluded these territories from the analysis. Some teams treated DC as both a county and a jurisdiction, so we excluded DC from the county forecasts. In addition, because county population size and transmission are correlated and case counts and forecast performance are also correlated, we grouped counties into 5 quintiles based on their population sizes, with cut points at 8,908; 18,662; 36,742; and 93,230 people; most analyses used forecasts from the quintile with the largest population size (n = 628). We hypothesized that small counties would be more likely to have better forecast accuracy because they had zero or very few reported cases. We thus chose to stratify counties by size to minimize any bias from aggregation. Performance results for most county forecasts are presented in S3 Appendix, and state-level spatial correlation is presented in S4 Appendix.

Defining epidemic phases

For every state and DC, we independently classified each forecast week based on the estimated time-varying reproduction number (R_t) for that given week. State-level R_t estimates were obtained from https://github.com/epiforecasts/covid-rt-estimates [49]. We extracted the R_t estimate for the Wednesday of each week from all available files. Because R_t estimates were updated on a rolling basis in near real-time, there were multiple estimates generated for the same date; we calculated the median estimated R_t per date for the upper and lower 90% credible interval and the median value (August 1, 2020 –January 15, 2022, or 2020 EW 31–2022 EW 2, reflecting 77 weeks in total). Each forecast week was then classified into one of the following categories based on the R_t estimates: increasing, peak, decreasing, and nadir.

Increasing and decreasing phases reflect weeks in which R_t had a 90% probability of being greater than or less than 1.0, respectively. There were several periods of rapid transmission in certain jurisdictions where R_t dipped above/below the 1.0 threshold but did not remain on an upward or downward trajectory. Thus, we classified weeks between two increasing phases as increasing and weeks between two decreasing periods as decreasing. Weeks between increasing and decreasing phases were classified as peaks, whereas nadirs were defined as periods between decreasing and increasing phases. Periods at the beginning or the end of an analysis period were classified as a continuation of whichever phase preceded or followed them. The proportion of weeks classified as each epidemic phase and graphical depictions of R_t are provided in S5 Appendix. The proportion of weeks suggests non-stable Rt trajectories in most locations and there is general concordance between R_t and reported cases.

Evaluation methodology

We evaluated probabilistic forecast accuracy using two different metrics, empirical prediction interval coverage rates and weighted interval scores (WIS) [43]. Coverage was calculated by determining the frequency with which the prediction interval contained the eventually observed outcome for the 50%, 80% and 95% intervals. WIS reflects a weighted estimate of sharpness (i.e., the range of the predicted interval) and calibration (i.e., precision or error) across the three prediction intervals and the median prediction, with higher WIS and indicating lower forecast skill. WIS also integrates measures of overprediction and underprediction, that is, the difference in the observed value and the lower or upper limit of the prediction interval. Importantly, WIS is highly correlated with the magnitude of observed and forecasted values. We used mean absolute WIS to assess forecast accuracy over time and mean relative WIS (rWIS) to assess forecast accuracy over space. Relative WIS was estimated by calculating the geometric mean of WIS across all sets of team forecasts and scaling that value to the WIS of a naïve model, the COVID-hub baseline. This approach eases interpretation, where values greater than 1.0 reflected worse accuracy than the baseline model and values below 1.0 reflected better model performance. Additionally, the pairwise relative comparison helps account for missing forecasts. Both coverage and WIS have been described in detail elsewhere [16,43]. Horizon-specific results for national, state/territory/DC, and large counties are presented in S6 Appendix.

To assess the association between WIS and epidemic phase for each team, we fitted separate Gaussian generalized estimating equation (GEE) models per team (Eq 1) with an independent working correlation structure at the state level. This structure assumes that observations are correlated within a state (denoted as l in the equations below), but not correlated over time in said state. Cases and weighted interval scores were log transformed and then standardized (subtracting the mean and dividing by the standard deviation) prior to fitting the model, as this transformation yielded more computationally and numerically stable estimates. We define those resulting variables as stdWIS and stdCases. The expected value for a standardized WIS for time (t) and location (l), with forecasts from a given team’s model, is as follows: (1)

Where p[t,l] is an index that reflects the phase of each time (t) and location (l), h is the horizon of the forecast in weeks, and ns(·) represents a natural spline with two degrees of freedom. The model intercept is represented by β₀ and error by ε_t.l.h. Using a regression model allows us to summarize patterns of overall average performance between teams while accounting for high correlation and variation in the scores. Comparisons of rWIS, in contrast, do not allow for formal inference with statistical hypothesis testing or interval-based inference. Prior to applying this regression model structure, our model building approach included exploratory analysis of several structures appropriate for longitudinal analysis. We examined model residuals, influential observations, goodness of fit metrics, and the impact of changing the functional form of the variables included in the model.

The inclusion of reported cases in models permitted flexible adjustment for the wide range in cases between and within jurisdictions, which led to a wide range of possible WIS values, as WIS values tend to be higher when counts are higher. Expected WIS values were computed by first obtaining a marginal mean from the GEE model and then undoing the transformations by exponentiating and un-standardizing the marginal mean. This was done separately for each team for all phases and for each team and each phase individually (see S7 Appendix for estimated team-specific marginal mean WIS relative to reported case counts). Additionally, we calculated whether the 80% confidence interval (based on Gaussian distributional assumptions) for each team’s expected WIS outcome (on the log-scale and normalized, as described above) was less than the average baseline model across all phases (i.e., the marginal mean WIS for the baseline model).

To determine the direction of the forecast predictions, each model’s 50% PI was compared to the last known incidence value. Forecasts were categorized as increasing predictions when the lower limit of their 50% PI was greater than the last known value, decreasing predictions when the upper limit of their 50% PI was less than the last known value, and stable/uncertain if their 50% PI contained the last known value.

All analyses were conducted using the R language for statistical computing (v 4.0.3) (50), and the following packages were used for the main analyses: scoringutils (44), covidhubUtils (51), geepack (52). Additionally, we included the EPIFORGE 2020 reporting guideline checklist in S8 Appendix to indicate each page in this evaluation that corresponds to each specific recommendation (13).

This activity was reviewed by the CDC and was conducted consistent with applicable federal law and CDC policy. See e.g., 45 C.F.R. part 46, 21 C.F.R. part 56; 42 U.S.C. §241(d); 5 U.S.C. §552a; 44 U.S.C. §3501 et seq.

CDC disclaimer: The findings and conclusions in this report are those of the authors and do not necessarily represent the official position of the Centers for Disease Control and Prevention.