Real-time forecasting of COVID-19-related hospital strain in France using a non-Markovian mechanistic model

Abstract

Projects such as the European Covid-19 Forecast Hub publish forecasts on the national level for new deaths, new cases, and hospital admissions, but not direct measurements of hospital strain like critical care bed occupancy at the sub-national level, which is of particular interest to health professionals for planning purposes. We present a sub-national French framework for forecasting hospital strain based on a non-Markovian compartmental model, its associated online visualisation tool and a retrospective evaluation of the real-time forecasts it provided from January to December 2021 by comparing to three baselines derived from standard statistical forecasting methods (a naive model, auto-regression, and an ensemble of exponential smoothing and ARIMA). In terms of median absolute error for forecasting critical care unit occupancy at the two-week horizon, our model only outperformed the naive baseline for 4 out of 14 geographical units and underperformed compared to the ensemble baseline for 5 of them at the 90% confidence level (n = 38). However, for the same level at the 4 week horizon, our model was never statistically outperformed for any unit despite outperforming the baselines 10 times spanning 7 out of 14 geographical units. This implies modest forecasting utility for longer horizons which may justify the application of non-Markovian compartmental models in the context of hospital-strain surveillance for future pandemics.

Fig 1. Structure of the underlying COVID-19 epidemic discrete time model with vaccination.

The three shaded areas represent the general population or community (left), the critically-ill hopitalised patients (center) and individuals removed by either recovery or death (right). The first subscript (i) of compartment densities (capital letters) indicates the age class, while the second denotes the time (in days) elapsed since the entry into the compartment—g, h, r, and u are thus the maximum number of days possible to remain in a compartment represented by contiguous boxes (see [21] for their computation). Individuals in S are susceptible, J are non-critically infectious, Y are infected and will eventually require hospitalisation, H are hospitalised in a critical care bed for a long stay, W are in non-ICU beds, R are recovered and D are deceased. S^V, J^V and Y^V are the vaccinated counterparts of S, J and Y. Λ_i is the daily force of infection (a probability). δ_i is the daily vaccination rate. μ_i, ψ_i, and θ_i are transition rates between compartments, the latter being reduced by the factor ν_C for vaccinated individuals. Arrows between boxes show the daily flow of individuals between compartments where dotted arrows occur with probability 1. For the sake of simplicity, only one age group is depicted here and only one of the two complementary probabilities is shown for each bifurcating transition.

Fig 2. Qualitative inspection of COVIDici forecasts.

A) Overlay of national forecasts of ICU occupancy produced by $\mathsf{\text{COVIDici}}$ with list of governmental interventions issued (national level only). B) All forecasters at the 4-week horizon plotted with the observed ICU occupancy for metropolitan France.

Fig 3. Artistic representation of the survivor bias that may occur after a restrictive governmental intervention.

The dashed curve is the counterfactual ICU occupancy that would have occurred if the intervention did not happen. The solid curve is the ICU occupancy that we observed because the intervention did occur.

Fig 4. Standard evaluation metrics for ICU occupancy.

A) Empirical coverage rate across all geographic units and forecast horizons (dashed line is optimal). Auvergne-Rhône-Alpes = ARA, Bourgogne-Franche-Comté = BFC, Bretagne = BRE, Centre-Val de Loire = CVL, Corse = COR, Grand Est = GES, Hauts-de-France = HDF, Île-de-France = IDF, Normandie = NOR, Nouvelle-Aquitaine = NAQ, Occitanie = OCC, Pays de la Loire = PDL, Provence-Alpes-Côte d’Azur = PAC. B) Mean weighted interval score (WIS) of all geographic units over time for the four-week forecast horizon scaled by the Naive model. Shaded areas represent 95% bias-corrected and accelerated (BCa) bootstrap confidence intervals with 10000 replicates per target date.

Fig 5. Median absolute error (AE) for COVIDici relative to other models.

Ratio = median AE for COVIDici / median AE for baseline model. 95% confidence interval (CI) is the bias-corrected and accelerated (BCa) bootstrap confidence interval generated by a nonparametric bootstrap with 10000 replicates. The p value is the smallest alpha such that 1 is not contained in the corresponding 1 − α CI. See Fig 4 for region code definitions.

Fig 6. Median weighted interval score of COVIDici relative to other models.

Ratio is the median weighted interval score (WIS) for COVIDici divided by median WIS using all other models as baselines. Each 95% confidence interval (CI) is based on the bias-corrected and accelerated (BCa) bootstrap confidence interval generated by nonparametric bootstrap with 10000 replicates. The p value is the smallest alpha such that 1 is not contained in the corresponding 1 − α CI. See Fig 4 for region code definitions.

Fig 7. Binary metrics for ICU occupancy overload by capacity threshold.

Metrics are conditioned separately for dates where ICU overload was observed or not observed respectively relative to arbitrary capacity thresholds. Thresholds are defined as a proportion of the peak ICU occupancy observed in a geographic unit. Shaded areas correspond to 95% bias-corrected and accelerated (BCa) confidence intervals based on 10000 replicates of a non-parametric bootstrap that assumes both spatial and temporal independence. Shaded areas completely below the blue line indicate that COVIDici statistically outperformed that model for that metric, threshold and forecast horizon.