Reporting delays: A widely neglected impact factor in COVID-19 forecasts

Abstract

Epidemic forecasts are only as good as the accuracy of epidemic measurements. Is epidemic data, particularly COVID-19 epidemic data, clean, and devoid of noise? The complexity and variability inherent in data collection and reporting suggest otherwise. While we cannot evaluate the integrity of the COVID-19 epidemic data in a holistic fashion, we can assess the data for the presence of reporting delays. In our work, through the analysis of the first COVID-19 wave, we find substantial reporting delays in the published epidemic data. Motivated by the desire to enhance epidemic forecasts, we develop a statistical framework to detect, uncover, and remove reporting delays in the infectious, recovered, and deceased epidemic time series. Using our framework, we expose and analyze reporting delays in eight regions significantly affected by the first COVID-19 wave. Further, we demonstrate that removing reporting delays from epidemic data by using our statistical framework may decrease the error in epidemic forecasts. While our statistical framework can be used in combination with any epidemic forecast method that intakes infectious, recovered, and deceased data, to make a basic assessment, we employed the classical SIRD epidemic model. Our results indicate that the removal of reporting delays from the epidemic data may decrease the forecast error by up to 50%. We anticipate that our framework will be indispensable in the analysis of novel COVID-19 strains and other existing or novel infectious diseases.

COVID-19 pandemic, reporting delays, epidemic forecasts, parametric optimization, SIRD compartmental model, time series analysis

Significance Statement

Accurate forecasts constitute the first step toward controlling an epidemic outbreak. However, the accuracy of such forecasts strongly depends on input data. Here, we develop a statistical framework to identify and remove reporting delays from the epidemic data. We use this framework to de-noise epidemic data of the first COVID-19 wave in eight regions worldwide. We demonstrate that the removal of reporting delays may increase the accuracy of epidemic forecasts by up to 50%. We anticipate that our framework will be invaluable in the analysis and forecasts of existing and future epidemics.

Introduction

The COVID-19 pandemic has been crippling the world’s health, economies, and quality of life for over 3 years. We are gradually understanding that COVID-19 is here to stay. Fast mutation rates of the virus and its overwhelming spreading capability make it extremely hard, if not impossible, to eradicate (1). COVID-19 is not the first and almost surely not the last pandemic to hit humanity. Therefore, in order to better prepare and withstand other contagious diseases in the future, we need to extract as many lessons as possible from the COVID-19 pandemic.

Public awareness is, arguably, the first line of defense against any infectious disease. Efficient collection of epidemic data and accurate epidemic forecasts allow for timely containment of the spread or flattening of the curve to win time for the development of pharmaceutical treatment methods. Due to the success of network epidemiology and the broad availability of data, significant advances in epidemic modeling and forecast methods (2–6) have been achieved. However, the accuracy of epidemic forecasts strongly depends on the accuracy and timeliness of the input data.

Are epidemic data—especially in the short-time period after the onset of the epidemic—devoid of inaccuracies? Multiple sources suggest the negative answer. Indeed, public health indicators are known to be biased, because they are sensitive to fluctuations in supply and demand for diagnostic testing (7–9). Likewise, inequalities in geographic accessibility have been documented to adversely affect the coverage and timelines of health indicators (10, 11). It is now well understood that ignoring epidemic data inaccuracies may lead to delays in deploying intervention policies resulting in dire consequences of an epidemic on the population’s health (12, 13).

One prominent challenge is the existence of temporal delays in data reporting (14, 15), defined as the time difference between the event when the person was affected by the virus and the time this event is accounted for. The reporting delays may occur due to many reasons, including patient hesitancy, medical testing delays, and reporting delays by public health authorities (16, 17). Consequently, reporting delays may vary not only across different diseases but also across different regions of interest. The median diagnosis delay for malaria is, for instance, approximately 4 days (18). Similarly, the research on the Middle East respiratory syndrome coronavirus (MERS-CoV) found a time difference between the symptom onset and confirmation of approximately four days (19). Hepatitis A, measles, and mumps data are usually reported after eight days (16). Reporting delays for some diseases are significantly longer: Hepatitis B, Shigellosis, and Salmonella can take on average 2–3 weeks (16). Recent COVID-19 studies reveal significant reporting delays of infections in China (20–25), Italy (26), Germany (27, 28), Singapore (29), the United States of America (30), and the United Kingdom (14).

While there is a plethora of scientific works aimed at the reconstruction of network data, including the reconstruction of network data for epidemic forecasts (31–33), we lack inference methods to assess and improve the accuracy of epidemic data itself. Instead, the prevalent epidemic forecast philosophy is focused either on data fusion (34–36)—or machine learning methods (37–39).

In our work, we develop a statistical framework to remove reporting delays. We demonstrate that the removal of reporting delays may significantly improve the accuracy of epidemic forecasts. Our statistical framework can be used in combination with any epidemic forecast method that takes infectious, recovered, and deceased epidemic data as input. Our work is organized as follows. We first present the evidence for the presence of reporting delays in the epidemic data extracted from the first COVID-19 wave. We then proceed to develop a statistical framework to remove reporting delays. After validating our framework on synthetic data, we move on to remove and analyze reporting delays in eight hotspots of the COVID-19 pandemic. We conclude our work with a discussion of the impact of delay removal on the accuracy of epidemic forecasts.

Results

Notation

Before presenting our findings, we introduce our notation for the epidemic data. Throughout the text, we operate with the infectious $I$⁠, recovered $R$⁠, and deceased $D$ data. Each dataset is a time series of values, each corresponding to a specific observation time. For brevity, we refer to the triplet of infectious, recovered, and deceased data as $Y=\{I,R,D\}$⁠. All values contained in the Y time series are fractions of individuals found in the corresponding state on a specific day. For instance, $I[k]$ corresponds to the fraction of individuals who are infectious on day k. Since COVID-19 reported data often is advertised in the form of changes in the number of epidemic cases, we find it convenient to introduce the daily changes in epidemic data as $\mathrm{\Delta }Y[k]$ for $Y=\{I,R,D\}$⁠. Further, in this work, we operate with reported epidemic data $\tilde{Y}$ and inferred data $\hat{Y}$⁠. Since reported and inferred data are expected to differ from the true data Y, we need to distinguish the three. We summarize our notation in Table 1.

Table 1.

Naming convention for the epidemic data, $Y=\{I,R,D\}$⁠.

Cumulative quantities	Quantity increments
Y: fractions of cases	$\mathrm{\Delta }Y$⁠: fractions of new cases
$\tilde{Y}$⁠: fractions of reported cases	$\mathrm{\Delta }\tilde{Y}$⁠: fractions of reported new cases
$\hat{Y}$⁠: fractions of predicted cases	$\mathrm{\Delta }\hat{Y}$⁠: fractions of predicted new cases

Evidence for reporting delays in epidemiological data

We begin the exposition by considering daily reports on the fractions of infected $\mathrm{\Delta }\tilde{I}$⁠, recovered $\mathrm{\Delta }\tilde{R}$⁠, and deceased $\mathrm{\Delta }\tilde{D}$ individuals in Spain. We used the daily epidemic reports in Spain to construct the fraction of infected individuals as $\tilde{I}[k]=\sum _{\ell =0}^{k-1}(\mathrm{\Delta }\tilde{I}[\ell ]-\mathrm{\Delta }\tilde{R}[\ell ]-\mathrm{\Delta }\tilde{D}[\ell ])$⁠. While both the infected $\tilde{I}$ and the deceased data $\tilde{D}$ indicate that the first COVID-19 wave in Spain peaked in April 2020, the exact timings of the two peaks are, nevertheless, different. As observed in Fig. 1a, reported new deceased cases $\mathrm{\Delta }\tilde{D}$ reached their peak on 2020 April 1st, while the highest fraction of infectious individuals $\tilde{I}$ was observed 22 days later on 2020 April 23. This observation is not specific to Spain: the peaks in the number of reported new deceased cases $\mathrm{\Delta }\tilde{D}$ precede those of infectious cases by more than 1 week in most regions, Fig. S1. We make similar observations for COVID-19 daily recovery reports $\mathrm{\Delta }\tilde{R}$⁠, which exhibit their peaks after those of the deceased cases $\mathrm{\Delta }\tilde{D}$⁠, Figs. 1a and S1. Furthermore, when plotted as a function of daily deceased cases $\mathrm{\Delta }\tilde{D}$⁠, infectious cases $\tilde{I}$ and daily recovered cases $\mathrm{\Delta }\tilde{R}$ form “loop” patterns, see Figs. 1b, c and S2.

Fig. 1.

Evidence for reporting delays in COVID-19 epidemic data. a) Reported infectious $\tilde{I}[k]$⁠, recovered $\mathrm{\Delta }\tilde{R}[k+1]$ and deceased $\mathrm{\Delta }\tilde{D}[k+1]$ cases for the first COVID-19 wave in Spain. The $k=0$ day corresponds to 2020 February 25. Note the difference in peaks of the reported data. b) and c) display pairwise color-coded scatter plots of $\mathrm{\Delta }\tilde{D}[k+1]$ vs. $\tilde{I}[k]$⁠, and $\mathrm{\Delta }\tilde{R}[k]$ vs. $\mathrm{\Delta }\tilde{D}[k]$ for Spain. Colors, from light to dark green, reflect different days in the data ranging, respectively, from $k=0$ to $k=80$⁠. The scatter plots in b) and c) form, respectively, clockwise and counterclockwise loop patterns. d)–f) display the epidemic data in Spain after the removal of reporting delays. g) Synthetic epidemic data generated with the SIRD model by solving Eq. 1 with parameters $\beta =0.5$⁠, ${\gamma }_r=0.2$⁠, ${\gamma }_d=0.05$⁠, $R[0]=D[0]=0$⁠, $I[0]={10}^{-6}$⁠, and $S[0]=1-I[0]$⁠. Synthetic reporting delays were generated with the Pólya-Aeppli distributions using Eq. 3 with parameters ${\lambda }_I=3.68$⁠, ${\theta }_I=0.5$⁠, $E[T_I]=7.36$⁠; ${\lambda }_R=9.46$⁠, ${\theta }_R=0.4$⁠, $E[T_R]=23.65$⁠; ${\lambda }_D=0.3$⁠, and ${\theta }_D=0.3$⁠, $E[T_D]=1$⁠. h) and i) display pairwise color-coded scatter plots of $\mathrm{\Delta }\tilde{D}[k+1]$ vs. $\tilde{I}[k]$⁠, and $\mathrm{\Delta }\tilde{R}[k]$ vs. $\mathrm{\Delta }\tilde{D}[k]$ for the SIRD epidemic data. j)–l) display the synthetic epidemic data after the removal of reporting delays.

The patterns observed in Figs. 1a–c, S1, and S2 may indicate the presence of reporting delays. Indeed, epidemic models view recoveries and deaths of patients as stochastic processes with effective rates proportional to the number of infected individuals I, $\mathrm{\Delta }R[k]\propto {\gamma }_{r}I[k]$⁠, $\mathrm{\Delta }D[k]\propto {\gamma }_{d}I[k]$⁠. Since the SARS-CoV-2 virus has hardly changed during the first wave of the pandemic, we expect that recovery ${\gamma }_r$ and death ${\gamma }_d$ rates are approximately constant during this period (40–42). This observation suggests that changes in the fractions of recovered $\mathrm{\Delta }R$⁠, and deceased $\mathrm{\Delta }D$ data are proportional to the fraction of infectious individuals $I$ and, therefore, reach maximum values at the same time step, contradicting Fig. 1a–c.

We hypothesize that the disagreement between $\mathrm{\Delta }\tilde{R}$⁠, $\mathrm{\Delta }\tilde{D}$⁠, and $\tilde{I}$ is due to reporting delays. If each $\mathrm{\Delta }\tilde{R}$⁠, $\mathrm{\Delta }\tilde{D}$⁠, and $\tilde{I}$ time series are reported with different delays, their peaks are expected to differ. Likewise, the loop patterns of Fig. 1b, c may also be the result of an effective time shift between two nonmonotonous time series.

To check if reporting delays may result in the observed patterns, we consider the compartmental Susceptible-Infectious-Recovered-Deceased (SIRD) epidemic model. Within the SIRD model, the population is split into four compartments: susceptible S, infectious I, recovered R, and deceased D. Compartment S denotes the fraction of susceptible individuals, who can be infected by infectious individuals. Compartment I denotes the fraction of individuals, who have been infected but have not recovered or are deceased. Compartments R and D are respectively the fractions of individuals, who have recovered or are deceased. The SIRD model assumes that recovered individuals become immune and cannot be infected by the virus in the future. Further, the SIRD model assumes the uniform mixing of the Infectious and Susceptible sub-populations. As a result, the discrete-time transitions between the compartments are governed by first-order difference time equations

$$\begin{array}{rl}& I[k+1]-I[k]=\beta I[k]S[k]-({\gamma }_r+{\gamma }_d)I[k],\\ & R[k+1]-R[k]={\gamma }_{r}I[k],\\ & D[k+1]-D[k]={\gamma }_{d}I[k],\\ & S[k]+I[k]+R[k]+D[k]=1,\end{array}$$

(1)

where β, ${\gamma }_r$⁠, and ${\gamma }_d$ are the infection, the recovery, and the deceased probabilities, respectively. We used Eq. 1 to generate infectious $I[k]$⁠, recovered $R[k]$⁠, and deceased $D[k]$ time series epidemic data, and then added synthetic delays to the time series data, see Fig. 1g–i and Supplementary Material B. Since reporting delays were manually added, we observe that $\mathrm{\Delta }R[k+1]$⁠, $\mathrm{\Delta }D[k+1]$⁠, and $I[k]$ reach maxima at different time steps, Fig. 1g.

Upholding our hypothesis, we observe the formation of loop patterns in the $\mathrm{\Delta }\tilde{D}[k+1]$ vs. $\tilde{I}[k]$ and $\mathrm{\Delta }\tilde{R}[k+1]$ vs. $\tilde{D}[k]$ in Fig. 1h, i, similar to those of Spain, Fig. 1b, c. The observed loop patterns are due to the effective horizontal shifts of the corresponding times series due to reporting delays. Indeed, let us split the observation time window into three windows formed by the maxima of the $\mathrm{\Delta }\tilde{D}[k+1]$ and $\tilde{I}[k]$ curves, as shown in Fig. 1g. In window i, both $\tilde{I}[k]$ and $\mathrm{\Delta }\tilde{D}[k+1]$ increase as a function of discrete time k. Since reporting delays in the synthetic $\mathrm{\Delta }\tilde{D}[k+1]$ data are smaller than those in $\tilde{I}[k]$⁠, this time window corresponds to the upper branch of the $\mathrm{\Delta }\tilde{D}[k+1]$ vs. $\tilde{I}[k]$ loop in Fig. 1h. In window ii, $\tilde{I}[k]$ increases while $\mathrm{\Delta }\tilde{D}[k+1]$ decreases. Thus, window ii corresponds to the top (decreasing) section of the $\mathrm{\Delta }\tilde{D}[k+1]$ vs. $\tilde{I}[k]$ loop, Fig. 1h. Finally, in window iii both $\mathrm{\Delta }\tilde{D}[k+1]$ and $\tilde{I}[k]$ decrease as a function of time step k resulting in the lowest section of the $\mathrm{\Delta }\tilde{D}[k+1]$ vs. $\tilde{I}[k]$ loop, Fig. 1h. Combined, all sections correspond to the loop pattern of Fig. 1h with points progressing in the clockwise direction. Similar considerations explain the counterclockwise loop pattern in the $\mathrm{\Delta }\tilde{R}[k+1]$ vs. $\tilde{D}[k]$ scatter plot. The counterclockwise progression of points in this loop pattern is due to $\mathrm{\Delta }\tilde{R}[k]$ lagging behind the $\mathrm{\Delta }\tilde{D}[k]$ time series, Fig. 1i.

A null model for reporting delays

To uncover reporting delays, we employ the following null model: each individual i is endowed with random event time ${\tau }_{Y_i}$ of either getting infected, ${\tau }_{I_i}$⁠, recovered, ${\tau }_{R_i}$⁠, or deceased, ${\tau }_{D_i}$⁠. Figure 2 depicts the random delay time for event $Y_i=\{I_i,R_i,D_i\}$ as $T_{Y_i}=\{T_{I_i},T_{R_i},T_{D_i}\}$ and the corresponding random reported time as $K_{Y_i}={\tau }_{Y_i}+T_{Y_i}$⁠. Within our discrete-time setting, these random times are discrete random variables. Assuming that the reporting delays $T_{Y_i}$ are independent of the events $Y_i$⁠, we obtain for the reported time $K_{Y_i}$

$$Pr[K_{Y_i}=k]=\sum _{y=0}^{\mathrm{\infty }}\sum _{m=0}^{\mathrm{\infty }}Pr[{\tau }_{Y_i}=y]Pr[T_{Y_i}=m]{\delta }_{y+m,k}=\sum _{m=0}^{k}Pr[{\tau }_{Y_i}=k-m]Pr[T_Y=m].$$

(2)

Fig. 2.

A schematic representation of epidemic events and corresponding delay times.

By taking the arithmetic mean of Eq. 2, we find that the expected reported fraction $\mathrm{\Delta }\tilde{Y}[k]=\frac{1}{N}\sum _{i}Pr[K_{Y_i}=k]$ of individuals in state Y is a discrete convolution

$$\mathrm{\Delta }\tilde{Y}[k]=\sum _{m=0}^{k}Pr[T_Y=m]\mathrm{\Delta }Y[k-m],$$

(3)

where $Y=\{I,R,D\}$⁠. Equation 3 serves as the foundation for uncovering reporting delays and improving epidemic forecasts.

Statistical framework to uncover reporting delays

The delay distribution $Pr[T_Y=m]$ in (3) can be determined, in principle, as the solution of a large set of quadratic equations, deduced from the governing epidemic equations, as shown in Supplementary Material D, provided that we assume a same distribution for each delay $T_{Y_i}=\{T_{I_i},T_{R_i},T_{D_i}\}$⁠. Since that solution is rather demanding, we propose a simplification. We assume that the delay distribution $Pr[T_Y=m]$ is generated by a two-parameter probability distribution with parameters $\kappa =\{\lambda ,\theta \}$⁠, which is sufficiently general to incorporate the real probability distribution of the delay.

We consider three families of two-parameter discrete probability distributions: the negative binomial distribution, the Neyman type A distribution, and the Pólya-Aeppli distribution, see Materials and methods A. These distributions are quite general and contain some well-known 1-parameter distributions as special cases. The logarithmic, geometric, and Poisson distributions are, for instance, all special cases of the negative binomial distribution (43). The three distributions are sufficiently general to generate data with variable mean and variance and of varying skewness (43, 44)—properties expected of the epidemic delay data, see Materials and methods A. For brevity, we focus on the Pólya-Aeppli distribution in the main text and report the results for the other two distributions in Fig. S4.

In our statistical framework, we assume that the reporting delays correspond to three datasets, $Y=\{I,R,D\}$⁠, which are all characterized by the Pólya-Aeppli distribution, albeit with different parameters ${\mathit{\kappa }}_Y\equiv \{{\lambda }_Y,{\theta }_Y\}$⁠. Hence, there are in total six parameters for three Pólya-Aeppli distributions, $\mathit{\kappa }=({\lambda }_I,{\theta }_I,{\lambda }_R,{\theta }_R,{\lambda }_D,{\theta }_D)$⁠. With the choice of the Pólya-Aeppli distribution, the reporting delays $\mathrm{\Delta }\tilde{Y}$ in our basic equation Eq. 3 can be determined in the parameterized form $\mathrm{\Delta }{\tilde{Y}}_{\mathit{\kappa }}$⁠.

Given the incremental time series $\mathrm{\Delta }{\tilde{Y}}_{\mathit{\kappa }}$⁠, the cumulative time series ${\tilde{Y}}_{\mathit{\kappa }}$ are given by

$$\begin{array}{rl}& {\tilde{I}}_{\mathit{\kappa }}[k+1]={\tilde{I}}_{\mathit{\kappa }}[k]+\mathrm{\Delta }{\tilde{I}}_{\mathit{\kappa }}[k]-\mathrm{\Delta }{\tilde{R}}_{\mathit{\kappa }}[k]-\mathrm{\Delta }{\tilde{D}}_{\mathit{\kappa }}[k],\\ & {\tilde{R}}_{\mathit{\kappa }}[k+1]={\tilde{R}}_{\mathit{\kappa }}[k]+\mathrm{\Delta }{\tilde{R}}_{\mathit{\kappa }}[k],\\ & {\tilde{D}}_{\mathit{\kappa }}[k+1]={\tilde{D}}_{\mathit{\kappa }}[k]+\mathrm{\Delta }{\tilde{D}}_{\mathit{\kappa }}[k],\end{array}$$

(4)

for $k\ge 0$⁠, and $\widetilde{I_{\mathit{\kappa }}}[0]=\widetilde{R_{\mathit{\kappa }}}[0]=\widetilde{D_{\mathit{\kappa }}}[0]=0$⁠.

The main assumption of our reporting delay removal framework is that the increments in new recovered $\mathrm{\Delta }R$ and deceased $\mathrm{\Delta }D$ individuals are proportional to the cumulative fraction of infectious individuals I.

Therefore, we determine the “best” parameters $\overline{\mathit{\kappa }}$ that maximize the product of pairwise correlations among the three epidemic time series in $\widetilde{Y_{\mathit{\kappa }}}$⁠:

$$O_b(\widetilde{Y_{\mathit{\kappa }}})\equiv O_b(\widetilde{I_{\mathit{\kappa }}},\mathrm{\Delta }\widetilde{R_{\mathit{\kappa }}},\mathrm{\Delta }\widetilde{D_{\mathit{\kappa }}})=\rho (\mathrm{\Delta }\widetilde{R_{\mathit{\kappa }}},\mathrm{\Delta }\widetilde{D_{\mathit{\kappa }}})\rho (\widetilde{I_{\mathit{\kappa }}},\mathrm{\Delta }\widetilde{R_{\mathit{\kappa }}})\rho (\widetilde{I_{\mathit{\kappa }}},\mathrm{\Delta }\widetilde{D_{\mathit{\kappa }}}),$$

(5)

where $\rho (X,Y)$ is the Pearson correlation coefficient between time series X and Y. The objective function $O_b(\widetilde{Y_{\mathit{\kappa }}})$ in Eq. 5 reaches its maximum value of 1 when all three pairwise correlations among the $\widetilde{I_{\mathit{\kappa }}}$⁠, $\mathrm{\Delta }\widetilde{R_{\mathit{\kappa }}}$ and $\mathrm{\Delta }\widetilde{D_{\mathit{\kappa }}}$ time series are 1, which we expect when recovery ${\gamma }_r$ and deceased ${\gamma }_d$ epidemic probabilities are constant. Due to the nature of the Pearson correlation coefficient, the objective function $O_b(\widetilde{Y_{\mathit{\kappa }}}[k])=O_b(\widetilde{Y_{\mathit{\kappa }}}[k-T])$ is invariant under the constant time shift T of the epidemic data. As a result, we can only infer the reporting delays up to a constant time shift T.

To maximize $O_b(\widetilde{Y_{\mathit{\kappa }}})$⁠, the random search (45) method is applied: we conduct a large set of independent random iterations, $\ell =1,\dots ,L$⁠. At each iteration ℓ, we select the elements of the parameter vector $\overline{\kappa }$ uniformly at random from the prescribed domain of values. Further, at each iteration ℓ, we use the selected parameter set ${\overline{\kappa }}_{\ell }$ to reconstruct the original epidemic data $\mathrm{\Delta }{\tilde{Y}}_{{\overline{\kappa }}_{\ell }}$ by solving Eq. 3 and convert incremental data $\mathrm{\Delta }{\tilde{Y}}_{{\overline{\kappa }}_{\ell }}$ to cumulative data ${\tilde{Y}}_{{\overline{\kappa }}_{\ell }}$ by solving Eq. 4. We then compute the pairwise correlations in ${\tilde{Y}}_{{\overline{\kappa }}_{\ell }}$ to obtain the objective function $O_b({\tilde{Y}}_{{\overline{\kappa }}_{\ell }})$⁠. After completing all random search iterations, the resulting delay parameter $\hat{\kappa }$ is the one maximizing the objective function, $\hat{\kappa }={\mathrm{argmax}}_{{\overline{\kappa }}_{\ell }}{O}_b({\tilde{Y}}_{{\overline{\kappa }}_{\ell }})$⁠, see Materials and methods B.

Tests on synthetic data

Before uncovering reporting delays in real epidemic data, we test our framework on synthetic datasets, which we generate with the SIRD compartmental model, Eq. 1.

To test our framework, we generated 50 SIRD epidemic model datasets with different epidemic parameters and added synthetic reported delays to the obtained times series, as prescribed by Eq. 3. The resulting delay-perturbed synthetic data is fed to our statistical framework that infers the reported synthetically added delays. Figure 1d–f and j–l display, respectively, the reconstructed real and synthetic epidemic data. We find strong correlations between the reconstructed $\mathrm{\Delta }\tilde{R}$⁠, $\mathrm{\Delta }\tilde{D}$⁠, and $\tilde{I}$ data. Further, Fig. 3a, b indicate that the inferred delay parameters are in good agreement with the true parameters that generated the datasets. Figure 3 illustrates that the inference errors decrease fast as a function of the number of iteration steps saturating at the value of 2 days after ${10}^7$ iterations. To assess the robustness of the inference procedure, we have conducted cross-inference experiments by generating synthetic data with one delay distribution and inferring delays using another distribution, arriving at similar results in Fig. S4.

Fig. 3.

Uncovering reporting delays in synthetic data. We generate epidemic data using the SIRD model, Eq. 1 with parameters $\beta =0.23$⁠, ${\gamma }_r=0.079$⁠, and ${\gamma }_d=0.11$ and variable delay parameters. The schematics on top illustrates the process of generation and inference on synthetic data. a) The mean absolute inference errors as functions of the number of iteration steps in the random search. Due to the nature of the Pearson correlation coefficient, it is only possible to infer relative delays, $E[T_R]-E[T_D]$⁠, $E[T_I]-E[T_D]$⁠, and $E[T_R]-E[T_I]$⁠. b) The relationship between the true and the inferred values of delay parameters. c, d) An example of an SIRD synthetic data c) before and d) after the removal of synthetic reporting delays.

Uncovering reporting delays in real data

After testing our inference framework on synthetic data, we moved on to uncover reporting delays in real epidemic data that we have collected from eight regions worldwide, Supplementary Material A. Figure 1d–f and j–l display the reconstructed epidemic data for Spain and for synthetic SIRD model after the identification and removal of reporting delays, see Table S1 for the inferred parameters. Consistent with our expectations, the reconstructed time series of infected $I[k]$⁠, new recovered $\mathrm{\Delta }R[k+1]$⁠, and new deceased $\mathrm{\Delta }D[k+1]$ individuals are strongly correlated, Fig. 1e, f and k, l, and their peaks co-occur, as seen in Fig. 1d, j.

We summarize the properties of uncovered reporting delays in the eight regions in Table 2. We find that the infectious and recovered data delays are longer than those for deceased cases and vary from several days to several weeks. This observation is hardly surprising. Indeed, there are several factors contributing to the delays of infectious cases. One factor is the delay between an individual becoming infectious, and the symptom onset (46). Another factor, particularly significant in the first COVID-19 wave, is the delay between the symptom onset and the test result (47). In their turn, delays in the recovered events are likely caused by hospital discharge policies. Indeed, the recovered data are usually derived from hospital discharge events, which occur after the patient’s recovery.

Table 2.

Inferred reporting delays.

Regions	$E[T_R]-E[T_D]$	$E[T_I]-E[T_D]$	${\sigma }_I$	${\sigma }_R$	${\sigma }_D$	$\rho (\mathrm{\Delta }\overline{R},\mathrm{\Delta }\overline{D})$	$\rho (\overline{I},\mathrm{\Delta }\overline{R})$	$\rho (\overline{I},\mathrm{\Delta }\overline{D})$	$O_b(\overline{Y})$
Italy	28.52	6.59	5.08	27.79	0.72	0.91	0.94	0.99	0.84
Spain	22.66	6.36	8.63	28.39	0.73	0.99	0.99	1.00	0.98
Wuhan	14.80	20.64	51.67	16.13	23.87	0.78	0.89	0.91	0.63
Turkey	14.13	24.85	43.89	12.47	11.21	0.91	0.95	0.98	0.85
Hubei	9.96	21.07	78.62	10.89	54.14	0.85	0.90	0.92	0.71
Romania	3.30	19.05	43.08	90.12	0.29	0.81	0.89	0.97	0.70
Germany	2.72	16.55	39.25	106.89	4.31	0.87	0.88	0.99	0.76
Denmark	0.17	25.07	36.90	2.71	28.25	0.83	0.93	0.94	0.72

The table displays the inferred differences between the expected delay times $E[T_R]-E[T_D]$⁠, $E[T_I]-E[T_D]$⁠, standard deviations ${\sigma }_I$⁠, ${\sigma }_R$ and ${\sigma }_D$⁠, and optimized values for the objective function $O_b(\overline{Y})$⁠. See Tables S1 and S2 for the corresponding parameters of the Pólya-Aeppli distributions.

Based on the expected delay values, one can naturally split the eight ROIs into three categories: (i) large infectious and small recovered delays: Romania, Germany, and Denmark, (ii) small infectious and large recovered delays: Italy and Spain, and (iii) large infectious and recovered delays: Wuhan, Hubei, and Turkey.

Small infectious delays imply that COVID-19 tests are timely and accurate. This seems to be the case for Italy, which executed more tests per capita in April 2020 than other countries (48). On the contrary, the testing ability of the Hubei province was significantly insufficient during the first wave of the COVID-19 outbreak.

Large delays in the recovered data may be attributed to strict hospital discharge policies. Discharge policy in Italy, for instance, was based on the negative test (49), and it has been shown that COVID-19 tests may stay positive for an extended time after COVID-19 symptoms disappear. In contrast, discharge policies in Denmark were based not on the negative test but on patient symptoms (49), likely leading to shorter reporting delays in the recovered data.

Large standard deviations in the reporting delays may indicate irregularities in reporting mechanisms. As an example, Hubei province expanded its daily testing capacity from 200 to 2,000 individuals from the beginning of the pandemic until January 27th (50). As a result, fewer individuals were tested late at the end of the first COVID-19 wave compared to its beginning, likely resulting in the large standard deviation of infectious data delays.

The small standard deviation in the delays of recovered data observed in Hubei and its capital Wuhan is likely the consequence of the strict discharge criteria (51). Although strict discharge policies cause significant delays, these delays are similar, resulting in relatively smaller ${\sigma }_R$ values. In contrast to Chinese regions, the recovery data for Germany are not reported directly but instead are estimated by a not explained algorithm (52), resulting in large errors and, consequently, larger ${\sigma }_R$ values. Based on the optimized values $O_b(\overline{Y})$ as shown in Table S2, the optimization performances of our algorithm for Italy, Spain, and Turkey are better than the other countries.

Improving epidemic forecasts

Is accounting for reporting delays likely to improve epidemic forecasts? To answer this question, we designed two experiments. Experiment 1 aims to forecast the epidemics ignoring reporting delays and serves as a baseline for Experiment 2, which uncovers reporting delays prior to forecasting epidemic data.

In both experiments, we split the reported epidemic data into two parts, which we call the training and the testing sets, respectively, see Fig. 4a. In Experiment 1, we fit the training set with the SIRD epidemic model, obtaining model parameters β, ${\gamma }_r+{\gamma }_d$⁠, and the fraction of initial infected cases $I[0]$⁠. We then use the SIRD model with the obtained parameters to forecast the epidemic data, Supplementary Material C.

Fig. 4.

Accounting for reporting delays improves epidemic forecasts. a) The schematic diagram for the two forecast experiments. b), c) display the results of epidemic forecasts in Spain with different forecast start dates. The forecasts of experiment 2 (yellow) are closer to the reported data (red) than the forecasts of experiment 1 (blue). The forecast results for all eight countries or regions are shown in Figs. S5 and S6. d), e) All forecast results for Spain and Germany. Shown are the RMSE forecast errors and their ratios. h) Average forecast errors $ϵ$ for the 8 regions. The benefits of removing reporting delays vary by region. g) The average RMSE forecast error of experiment 2 as a function of SIRD model fitting error for all experiments. h) The mean relative forecast error as a function of SIRD model fitting error for all experiments.

In Experiment 2, we first use the reported data $\tilde{Y}$ in the training set to infer the parameters of the delay distributions $\overline{\mathit{\kappa }}$⁠. We rely on these parameters to remove reporting delays from the testing set and obtain reconstructed data $\mathrm{\Delta }\overline{I}$⁠. In the next step, we fit the reconstructed data $\mathrm{\Delta }\overline{I}$ with the SIRD model, obtaining β, ${\gamma }_r+{\gamma }_d$ and $I[0]$ spreading parameters. We use these spreading parameters to forecast the epidemic data $\mathrm{\Delta }\hat{I}$⁠, which we compare to the reported data $\mathrm{\Delta }\tilde{I}$ in the testing set. Since $\mathrm{\Delta }\tilde{I}$ in the testing set contains the reporting delays, while the forecast $\mathrm{\Delta }\hat{I}$ does not, we added the inferred reporting delays back to the $\mathrm{\Delta }\hat{I}$ using Eq. 3, obtaining $\mathrm{\Delta }{\hat{I}}^{*}$⁠, See Fig. 4a and Supplementary Material C.

Figure 4b, c presents the results of the forecast experiments for Spain, indicating that correcting for reporting delays does improve the forecast accuracy. To evaluate the forecast errors in a systematic way, we evaluated the root mean square errors between the forecasts $I_F$ and the testing $I_T$ sets:

$$\mathrm{RMSE}(I_F,I_T)={[\frac{1}{n}\sum _{k=0}^{n-1}(I_F[k]-I_T[k]{)}^2]}^{\frac{1}{2}},$$

(6)

where n is the size of the testing set. Further, to quantify the benefits of accounting for reporting delays, we used the ratio of the root mean square errors measured in forecasts with and without reporting delays, $ϵ\equiv \frac{\mathrm{RMSE}({\hat{I}}_2^{*},I)}{\mathrm{RMSE}({\hat{I}}_1,I)}$⁠, where ${\hat{I}}_1$ and ${\hat{I}}_2^{*}$ are the epidemic forecasts obtained in experiments 1 and 2, respectively. The smaller the ratio, the smaller the relative forecast error.

We measured the accuracy of epidemic forecasts with different start dates for Spain in Germany. While we did observe a substantial improvement in forecast accuracy for Spain, Fig. 4d, this was not the case for Germany, where the removal of reporting delays only improved the forecasts by a small margin, Fig. 4e. On a broader scale, we observed that the benefits of correction for reporting delays vary across all regions of interest, Figs. 4f, S5, and S6. While there is nearly a two-fold improvement in the forecast accuracy for Denmark, there is little improvement for Wuhan and Hubei.

Apart from only measuring and testing a small part of the population—which might be too small or not sampled adequately—there are at least two factors that may hinder epidemic forecast accuracy. The first factor is the insufficient accuracy of reporting delay removals. The second one is the inability of the SIRD model to reproduce the COVID-19 dynamics accurately. We can quantify the former by the maximum attained value of the objective function $O_b(\overline{Y}{)}_{\max }$ that we aim to maximize when removing reporting delays. While there is no direct way to quantify the goodness of the SIRD model in epidemic forecasts, as an indirect measure, we fit the training set with the SIRD model and compute the fitting error.

Figure 4g shows that the forecast error is strongly correlated with the SIRD fitting error, Pearson $r=0.89$⁠, $p<{10}^{-27}$⁠, indicating that the highest epidemic forecast accuracy is attained when the model fit is accurate. The forecast error ratio depends both on the $O_b(\overline{Y}{)}_{\max }$ and SIRD model fitting, Fig. 4h. We observe the largest error ratio $ϵ$ in Wuhan, Hubei, and Romania, Fig. 4h. Wuhan and Hubei correspond to the largest SIRD fitting errors, while Romania corresponds to the lowest SIRD fitting errors. At the same time, all three regions are quantified by the lowest $O_b(\overline{Y}{)}_{\max }$ values. The other five regions, corresponding to smaller error ratios, are characterized by significantly larger $O_b(\overline{Y}{)}_{\max }$ values, Fig. 4g and Table S2.

These observations indicate that $O_b(\overline{Y}{)}_{\max }$ may be used as an early indicator of the benefit of delay removals. Indeed, lower $O_b(\overline{Y}{)}_{\max }$ values may indicate that reporting delays were not removed successfully or the initial assumption of the proportionality among I, $\mathrm{\Delta }R$⁠, and $\mathrm{\Delta }D$ times series does not hold. We note that high $O_b(\overline{Y}{)}_{\max }$ values correspond to high benefits of delay removals regardless of the SIRD fitting error. In the case of Romania, for instance, SIRD fitting errors are among the smallest, resulting in low forecast errors. Yet, the removal of reporting delays does not result in even lower forecast errors.

There might be multiple reasons for suboptimal reporting delay removals in the case of Hubei, Wuhan, and Romania. One possibility is the nonstationary nature of delays. The main assumption of our framework is that delay times are independent and drawn from the same distribution. While COVID-19 did not significantly mutate over the short time span of the first wave, our ability to handle the infections improved significantly. The PCR test capacity in China has grown remarkably during the first COVID wave, possibly explaining the limited effect of reporting delays on improving epidemic forecasts in Hubei and Wuhan.

Discussion

Data delays are ubiquitous in data sciences and adversely affect data analysis (53). The unique property of the epidemic data, enabling us to identify and filter out reporting delays, is the proportionality between the number of infectious individuals I and the rates of change in the deceased $\mathrm{\Delta }D$ and the recovered $\mathrm{\Delta }R$ individuals.

We relied on this proportionality property to develop a parametric statistical framework to uncover reporting delays in the first COVID-19 wave and applied it to eight regions significantly affected during the first COVID-19 wave. The character of uncovered delays varies across studied regions and can be explained by region-specific medical capacities and epidemic policies.

Concerning the epidemic forecasts, we found that the benefits of using curated epidemic data, as opposed to the raw data, are maximized in situations when reporting delays are removed efficiently. One of the main factors hindering the efficiency of delay removal—the nonstationarity of reporting delays—can be caused by either varying spreading properties of the virus or by rapid changes in medical capacities or epidemic restrictions across the regions. While the spreading properties of COVID-19 did not change significantly during the first wave of the pandemic, both medical capacities and epidemic policies were the subjects of constant updates as the communities learned how to handle the pandemic.

We expect that our framework will prove most useful in the early stages of an epidemic when accurate epidemic forecasts are essential to plan intervention strategies and raise public awareness. Common sense dictates that data collection and reporting challenges are the most significant during that time.

Our statistical framework is based on the assumption that transition probabilities ${\gamma }_r$ and ${\gamma }_d$ quantifying transitions from the infectious to the recovered, $I\to R$ and the infections to the deceased, $I\to D$ states, respectively are constant across the population and do not change over the observational period. While these assumptions are approximations, we expect them to hold better in early epidemic stages than in later stages. In the later stages, the population becomes more heterogeneous, e.g. due to natural or vaccine-induced immunity (54). In these situations, our framework can be generalized by splitting the infectious compartment I into sub-compartments and considering different rates, e.g.

$$\begin{array}{rl}& \mathrm{\Delta }R[k]={\gamma }_{1,r}{I}_1[k]+{\gamma }_{2,r}{I}_2[k],\\ & \mathrm{\Delta }D[k]={\gamma }_{1,d}{I}_1[k]+{\gamma }_{2,d}{I}_2[k],\\ & I[k]=I_1[k]+I_2[k],\end{array}$$

(7)

where $I_1[k]$ and $I_2[k]$ are two infectious sub-compartments characterized by different recovery and death probabilities.

At the same time, we stress that our statistical framework does not make strong assumptions about the spreading mechanisms of a pathogen and can be used in combination with any forecasting method that requires infectious, recovered, and deceased data as input.

In conclusion, our framework can be used as a preliminary filter for any epidemic forecast tool that takes infected, recovered, and deceased data as input. Our framework holds for any compartmental epidemic model, provided that each compartment can be measured, because a time series of each compartment is indispensable. Accurate and timely epidemic forecasts are of immense value for society and policymakers to minimize the adverse effects of the virus.

Materials and methods

Types of distributions

To determine the family of reporting delay distributions that best suit our data, we consider three different two-parameter discrete distributions (55) below:

(I) Negative binomial distribution.

$$Pr[T=m]=(\genfrac{}{}{0ex}{}{m+r-1}{m})(1-p{)}^m{p}^r.$$

(8)

The negative binomial distribution with parameters $r>0$ and $p\in [0,1]$ has mean value $E[T]=r(1-p)/p$ and variance $\mathrm{Var}[T]=r(1-p)/p^2$⁠.

(II) Pólya-Aeppli distribution is also known as the geometric Poisson distribution.

$$Pr[T=m]=\{\begin{array}{cc}{\displaystyle \sum _{j=1}^m{e}^{-\lambda }\frac{{\lambda }^j}{j!}(1-\theta {)}^{m-j}{\theta }^j(\genfrac{}{}{0ex}{}{m-1}{j-1}),}& m>0\\ e^{-\lambda },& m=0\end{array}.$$

(9)

The Pólya-Aeppli distribution with parameters $\lambda >0$ and $\theta \in [0,1]$ has mean value $E[T]=\lambda /\theta$ and variance $\mathrm{Var}[T]=\lambda (2-\theta )/{\theta }^2$⁠.

(III) Neyman type A distribution.

$$Pr[T=m]=\frac{{\mu }^m{e}^{-\xi }}{m!}\sum _{j=0}^{\mathrm{\infty }}\frac{{(\xi e^{-\mu })}^j}{j!}{j}^m.$$

(10)

The Neyman type A distribution with parameters $\xi >0$ and $\mu >0$ has mean value $E[T]=\xi \mu$ and variance $Var[T]=\xi \mu (1+\mu )$⁠.

Inferring reporting delays

To uncover reporting delays in epidemic data, we determine parameters of delays distributions ${\mathit{\kappa }}_Y$ maximizing the pairwise correlations between the epidemic time series, Eq. 5. This optimization problem can be compactly written as:

$$\begin{array}{rl}\underset{\mathit{\kappa }}{\text{arg max}}& O_b(Y)\equiv \rho (\mathrm{\Delta }R,\mathrm{\Delta }D)\rho (I,\mathrm{\Delta }R)\rho (I,\mathrm{\Delta }D)\\ \text{s.t.}& \mathrm{\Delta }I={\mathrm{\Psi }}_I^{-1}\mathrm{\Delta }\tilde{I},\mathrm{\Delta }R={\mathrm{\Psi }}_R^{-1}\mathrm{\Delta }\tilde{R},\mathrm{\Delta }D={\mathrm{\Psi }}_D^{-1}\mathrm{\Delta }\tilde{D},\\ & min(\mathrm{\Delta }I[k],\mathrm{\Delta }R[k],\mathrm{\Delta }D[k],I[k])\ge 0,\text{for}k=1,\dots ,T.\end{array}$$

(11)

Here T is the size of epidemic time series, $\tilde{Y}\equiv \{\tilde{I},\tilde{R},\tilde{D}\}$ and $Y\equiv \{I,R,D\}$ are reported and reconstructed epidemic data, respectively, while $\mathrm{\Delta }Y={\mathrm{\Psi }}_Y^{-1}\mathrm{\Delta }\tilde{Y}$ are the matrix solutions of Eq. 3. Indeed, Eq. 3 can be written in the matrix form as $\mathrm{\Delta }\tilde{Y}={\mathrm{\Psi }}_Y\mathrm{\Delta }Y$ for $Y=\{I,R,D\}$⁠, where

$$\begin{array}{rl}& {{\mathrm{\Psi }}_Y}_{i,j}\triangleq \{\begin{array}{cc}Pr[T_Y=i-j]& \mathrm{if}i\ge j.\\ 0& \mathrm{otherwise}.\end{array},\end{array}$$

Then, ${\mathrm{\Psi }}^{-1}$ is the inverse of ${\mathrm{\Psi }}_Y$ and $\mathrm{\Delta }Y={\mathrm{\Psi }}_Y^{-1}\mathrm{\Delta }\tilde{Y}$⁠.

In the main text, we assume that reporting delays are iid random variables drawn from three Pólya-Aeppli distributions, Eq. 9, with distinct parameters $\{{\lambda }_I,{\theta }_I\}$⁠, $\{{\lambda }_R,{\theta }_R\}$⁠, and $\{{\lambda }_D,{\theta }_D\}$⁠. In the case of Pólya-Aeppli distributions, the parameter vector takes the form of $\mathit{\kappa }\equiv \{{\lambda }_I,{\theta }_I,{\lambda }_R,{\theta }_R,{\lambda }_D,{\theta }_D\}$⁠. We solve the optimization problem given by Eq. 11 and the random search (45). For each iteration $\ell =1,\dots ,L$⁠, we treat the expected values of the Pólya-Aeppli distributions, $E[T_Y]={\lambda }_Y/{\theta }_Y$⁠, as iid random variables and draw from the uniform pdfs $U[0,30]$ for $Y=\{I,R,D\}$⁠. Similarly, we draw ${\theta }_Y$ parameters independently at random from uniform pdfs $U[0,1]$ and then determine ${\lambda }_Y$ parameters as ${\lambda }_Y={\theta }_{Y}E[T_Y]$⁠, obtaining ${\mathit{\kappa }}_{\ell }$⁠. In our experiments, we set the maximum number of search iterations to $L={10}^7$⁠.

For each ${\mathit{\kappa }}_{\ell }$⁠, we use Eq. 3 to reconstruct original data $Y_{\ell }$⁠, which we then use to compute the objective function $O_b(Y_{\ell })$⁠. After completing all iterations, the thought parameter vector $\hat{\kappa }$ describing reporting delays is the one corresponding to the maximum $O_b(Y_{\ell })$ value, $\hat{\kappa }={\text{arg max}}_{\mathit{\kappa }}{O}_b(Y)$⁠.

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

This research has been funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 101019718). L.M. is thankful for the support from the China Scholarship Council. M.K. has been supported by the Dutch Research Council (NWO) grant OCENW.M20.244.

Author Contributions

L.M. conceived the research idea. L.M., P.V.M., and M.K. designed research. L.M. and Z.Q. performed research. All authors discussed results and wrote the manuscript.

This manuscript has been posted as a preprint in the ArXiv repository https://arxiv.org/abs/2304.11863.

Data Availability

All epidemic data used in this work are publicly available from the original sources, see Supplementary Information A. The extracted epidemic time series for the eight regions of interest are deposited in FigShare (https://doi.org/10.6084/m9.figshare.22639519.v1). Code to generate SIRD epidemic time series and to uncover reporting delays is available on GitHub (https://github.com/qzhszl/Reporting-delays-a-widely-neglected-impact-factor-in-COVID-19-forecasts.git).

Author notes

Competing Interest: The authors declare no competing interest.