Incidence moments: a simple method to study the memory and short term forecast of the COVID-19 incidence time-series

Introduction

Since December 2019 the world has been going through the greatest global health crisis of the last decades. The COVID-19 pandemic had already produced more than 267million cases and about 5.3 million deaths at the end of December 2021. Each country that faces the threat of COVID-19 tries to inform its policy response using forecasting tools to determine the magnitude of the epidemic wave and the capacity of its health systems to plan ahead for the heavy demand for healthcare (WHO-Europe 2020). In a way, the quality of the response to the epidemic will depend on the ability of countries to adequately monitor and predict behavior of the epidemic, extracting useful information and implement timely effective public health interventions to reduce transmission (Ferguson et al. 2020).

However, the predictability of any epidemic is very uncertain and even more so in the case of an emerging new pathogen (Mangiarotti et al. 2020). Many models have been proposed, mainly of the SEIR type, considering social, biological, spatial and demographic characteristics in different countries, trying to predict the course of epidemics locally. An important limitation for the successful implementation of these models are data requirements and advanced mathematical modelling and epidemiological capacities, something that is not guaranteed in most settings within of low-and-middle income countries. However, even with complex models, the predictability show to be limited in the long term (Chen and Yu 2020; Jones and Strigul 2021). This is due, on the one hand, to the environmental, social, and demographic differences of the countries that affect the estimation of the different parameters, and, on the other hand, because the interaction dynamics between cases and susceptible may be non-linear reducing the predictability of the disease dynamic. The ability to predict COVID-19 dynamic has been very low, surprising us with sudden and unexpected changes in the number of cases. A recent study proposed that the dynamics of COVID-19 shows a chaotic behavior that would complicate the prediction of its dynamics (Chen and Yu 2020; Jones and Strigul 2021).

The epidemic is in a critical stage in Chile, with a high number of cases that already exceeds one million, with close to 30 thousand deaths and daily incidence rates that exceed 30 cases/10⁵ inhabitants, which has resulted in ICU occupation that exceeds 90% and repeatedly threatens to collapse the health system. Efforts have been made to monitor and model the course of the epidemic (Canals et al. 2020; Canals, Cuadrado, and Canals 2021; Córdova-Lepe, Gutiérrez-Aguilar, and Gutiérrez-Jara 2020; González et al. 2020; Guerrero-Nacuante and Maríquez 2020; Gutiérrez-Aguilar et al. 2020; Rojas-Vallejos 2020) proposing indicators that help decision-making. However, it is imperative to have predictions that allow the authorities to be alerted to the possible course of the epidemic and the overload of the health system. Maximum incidence-load models have been proposed (Canals, Cuadrado and Canals 2020) based on a previous short-term forecast models used for AH1N1 influenza (Canals 2010).

In order to have an indicator of the possible course of the epidemic, we study the temporal memory of the reported daily incidence time series and propose a simple model for short-term forecast of the incidence, which may be useful to estimate future ICU occupancy load. The simplicity of the proposed method is appealing since it allows a rapid implementation in different settings, even with limited technical epidemiological capacities and without requiring extensive computational power.

Methods

We conducted a study based on official daily public reports from the Ministry of Health of Chile, including daily new confirmed cases of COVID-19 at the national and sub-national levels (administrative regions) (MINSAL 2020). The series of daily cases in the world obtained from the daily information of the WHO (WHO 2021) was also analyzed, calculating the incidence under the assumption of a constant world population of 7 billion people. We analyzed these incidence time-series, proposing a model to study the memory of the time series, then we apply the method to incidence time series of 43 countries of the world from data gathered by Our World in Data (Our World in Data 2021) to assess its predictive capacity in different settings. These countries were selected based on the COVID data transparency index (CDTI index) (Total Analysis 2021). The countries with a CDTI index >40% were selected, as well as India, which although it has a lower index, is of great importance due to its high incidence and population size. Countries with less than 5 million people were not included and countries with registration problems or unreliable records, such as those with many zeros, were not analyzed. Besides Chile, of 58 countries with available data, 43 countries were finally analyzed. A period between March 24, 2020 and March 28, 2021 was analyzed for Chile, its regions, and for the 43 countries. For the global incidence time-series, a period between January 21, 2020 and May 15, 2021 was considered. We also included an incidence time series of Chile after reaching a 66% vaccination coverage of the population, from August 3 to October 28, 2021.

Model: To study the memory of the time series of reported incidence, we consider the following assumption: If on a day “t” there is I _t incidence and with an effective reproductive number R _t , then under the scenario of a population of constant size, it is expected that after a serial interval (τ), on t + τ, there is an incidence I _t R _t . If also R _t = R is a constant, in other words, the epidemic continues with its same growth rate (memory), at t + 2τ an incidence I _t R ² would be expected, and at t + 3τ an incidence I _t R ³ etc. The dependency of Rt overtime is therefore an element that is used to produce short-term forecast of the COVID-19 incidence.

Definition: We define IR ⁿ as the nth moment of incidence, where I is the moving average of the incidence of the last 7 days (I=(I _t  + I _t−1 +  … I _t−6)/7). We propose that the nth moment of incidence will be a predictor of the incidence reported in “nτ” days. Thus, for example, the second moment of incidence (IR ² ) will be the expected incidence in two serial intervals, 2τ=10 days, under the assumption of incidence and an effective Rt will be constant (R).

Statistical Methods. We estimate the effective reproductive number (Rt) with the method of Cori et al. (2013) and with the RKI method (an der Heiden and Hamounda 2020; Cori et al. 2013), with a serial interval of 5 days (Alenne et al. 2021; Izadi et al. 2020; Zhang, Wang, and Xie 2020) and we calculate the first five moments of the incidence (IR, IR ² , IR ³ , IR ⁴ and IR ⁵ ) representing predictions for 5, 10, 15, 20 y 25 days respectively. We compared the observed values with those expected for each of the first five moments using the mean arctangent absolute percentage error (MAAPE). MAAPE = 100 1 N ∑ t = 1 N arctan I t − F t A t , where N is the number of data points, It and Ft are the actual and forecasted incidence rates. This metrics has been proposed as an alternative for the mean absolute percentage error (MAPE) avoiding the disadvantage that MAPE produces infinite or undefined values for zero or close to zero actual values (Kim and Kim 2016). We consider MAAPE values ≤30 good, and values between 30 and 50 reasonable forecasting. This analysis was carried out for the country and for each of the administrative regions of Chile, and also for the incidence series of the world. The moment with largest long-term predictive capacity, with MAAPE≤30 for Chile and the world, was selected (α ^th moment (IR ^α )). The predicted incidence series for Chile was constructed using the chosen incidence moment.

We studied the correlation between population size and CDTI index with MAAPE value with Spearman correlation coefficient (Rs).

Results

We did not find substantial differences between the estimates of the effective reproductive number with the Cori et al. (2013) method and with the RKI method, so for simplicity we show the results using RKI (Suppl Inf.1).

For the Global and Chile series, we found good forecasting between the observed incidence values and the expected values for the first three moments, decreasing significantly for the fourth and fifth moments in Chile and for the fifth moment in the World. It can be seen that the third moment (IR ³ ) it is the one with a long-term predictive capacity (15 days). We do not find differences in predictive value between pre and post vaccine periods in Chile. The MAAPE for IR ³ were 30.5 and 28.7 respectively (Figure 1).

Figure 1:

Mean arctangent absolute percentage error (MAAPE) between the reported and the expected incidence of COVID-19 in the World and Chile, by means of the incidence moments method for 5 (IR), 10 (IR ² ), 15 (IR ³ ), 20 (IR ⁴ ) and 25 days (IR ⁵ ).

At the level of administrative regions of Chile, most of them with population sizes lower than one million people, we found only reasonable forecasting up to 15 days in the larger regions (Table 1). We choose the third moment of incidence (IR ³ ) as the longest-term predictor for Chile, without compromising its ability to predict future cases (Figures 2 and 3).

Figure 2:

Evolution of the reported incidence (blue) and the expected values of Chile, by the third incidence moment method (IR ³ ) (red), from March 28, 2020 to March 28, 2021.

Figure 3:

Evolution of the reported incidence (blue) and the expected values of eight selected regions of Chile, from North to South, by the third incidence moment method (IR ³ ) (red) from March 28, 2020 to March 28, 2021.

In 24 of 44 countries we found reasonable forecasting with IR ³ (Table 2, Figure 4). We found a lower predictive capacity (i.e. larger MAAPE) for countries with low population size (Rs = −0.60, p < 0.05). We did not find a strong correlation between MAAPE and CDTI index (Rs = 0.21, p > 0.05).

Figure 4:

Evolution of the reported incidence (blue) and the expected values of the World and the seven countries with better correlations (excepting Chile), by the third incidence moment method (IR ³ ) (red) from March 28, 2020 to March 28, 2021.

Discussion

We proposed the concept of incidence moments, derived from the definition of the effective reproductive number as the average number of new cases that a case produces in a serial interval. Thus, the expected cases with a time interval equal to a serial interval elapses will be the product of the number of cases and the effective reproduction number. As the process occurs very quickly, the population can be considered invariant and this concept can be extended to incidence rates (new cases per population at risk and time). As the serial interval of COVID-19 is short, estimated in this study as 5 days, IR estimates the future incidence in τ days. We used the moving average of the last 7 days as an estimator of the current incidence, which dampens the periodic oscillations in the weekly report due to weekend testing and reporting biases (Bukhar et al. 2020) and occasional stochastic fluctuations due to holidays or unexpected events. The same occurs with the calculation of effective reproduction number Rt with the RKI method and with the Cori et al. method, which in our case used a lag of 7 days (an der Heiden and Hamouda 2020; Cori et al. 2013). Thus, we propose that IR ⁿ is a proxy for the incidence in “nτ” days. This assumes that I and R will remain constant in this period, which is not necessarily true. We note that IR ⁿ depends linearly on I, but it depends on the nth power of the reproductive number, which makes it very sensitive to the estimation and variations of this parameter. The higher the value of the reproductive number, the greater the sensitivity of the predicted incidence. The time for which the prediction is valid (nτ) depends on the serial interval. Although there are estimates of the serial interval between 1.9 and 7.5 days (Alenne et al. 2021; Izadi et al. 2020; Zhang, Wang, and Xie 2020), most of the estimates are in a value close to 5 days. For example, in the meta-analysis of Izadi et al. (2020) an average value of 4.45 days was estimated, in that of Zhang, Wang, and Xie (2020) a value of 5.35 days and in that of Alenne et al. (2021) it is reported 5.2 days, so the value of 5 days used in this study is reasonable. However, the prediction time is sensitive to this parameter and variations in the serial interval can affect the validity of the prediction. For example, it has been reported that the serial interval can be shortened over time by non-pharmaceutical interventions (Alí et al. 2020). Additionally, if emerging variants modify the infectious dynamic of the agent (e.g., shorter incubation periods or increased infectivity in pre-symptomatic period) this can have an impact on serial interval, therefore affecting the predictive capacity of the method if this parameter is specified incorrectly.

Our results show that the prediction capacity, represented by the MAAPE metrics between the observed incidence and the expected values for 5, 10, 15, 20 and 25 days according to the first five moments, decreases notably for the fourth and fifth moments, particularly when the population size of the analyzed region is small. This occurs for two reasons, because of the dependence on a high power of the reproductive number and because of the low probability that both the incidence and the reproductive number remain constant for so long. Thus this system has a very short temporary memory.

The incidence remains stable over time in an endemic state with Rt=1. If Rt is different from 1 but constant, the incidence will increase or decrease but in a predictable way. However, with a variable reproductive number, the ability to predict the future is quickly lost. In our case, the predictive capacity using the third moment of incidence remained relatively high until day 15 (three serial intervals). However this was dependent on the size of the population. When we consider the data of the world and those of Chile, with more than 19 million people, this effect is very clear, but when we analyze each administrative division the high predictive capacity is maintained only in the Metropolitan Region that includes the capital of Chile, which contains more than 40% of the population of Chile. However, although the predictive capability are low in the regions with a smaller population, a clear decrease in predictability is found for the fourth and fifth order moments, with correlations that do not even reach reasonable forecasting.

The graphs showed that over estimations occur mostly in the high incidence periods. The graphs also show that the over-estimates occur after abrupt accelerations in the infection rate, that is, after significant increases in the effective reproductive number. This is explained by the aforementioned dependence on the third power of R. In several cases, such rapid changes in Rt could be explained by problems with data reporting, a limitation that is inherent to surveillance of an ongoing epidemics for and emerging disease.

For 24 analyzed countries (54.5%), we found 15-day reasonable forecasting capacity with our method. The performance is particularly good for India, Indonesia, Russia and The United States, countries with high population sizes. Also when we used this method for the 15-day prediction of the global incidence, under the assumption of a constant population of 7 billion individuals, we obtained good forecasting. This good value was not an expected result, given that there is great heterogeneity and asynchrony in infections between the different countries that make the assumption of population homogeneity for calculating the effective reproductive number impossible. It is probable that in this case the law of large numbers operates, and local variabilities are dampened, converging to an average value. Changes in incidence in one direction in a given country are offset by changes in the opposite direction in another country, making the predicted incidence equivalent to that which would occur in a homogeneous population.

We observed that occasionally the value of the incidence was overestimated in the COVID-19 global incidence and also in Chile and in the Chilean regions. This occurred when there were sharp fluctuations in incidence. For example if there is a drop in the incidence, or actual reporting problems, and then a sharp rise of this, the value of Rt increases abruptly and the prediction dependent on R ³ , overestimates the observed incidence.

The rapid loss of predictive capacity over time can be explained because the dynamics between susceptible and infected populations is non-linear and sensitive to initial conditions (Anderson and May 1992; Jones and Strigul 2021; Mangiarotti et al. 2020; Schaffer and Kott 1985). This loss of predictive capacity appears to be intrinsic to the COVID-19 dynamics without any detectable changes in predictability without and with high vaccination coverage (≥66%) in Chile. Low predictability associated with chaotic dynamics has been detected in the dynamics of many infectious diseases (Canals 1996; Canals and Labra 1997; Canals and Labra 1999; Grenfell et al. 1995; Olsen and Schaffer 1990; Olsen, Truty, and Schaffer 1988). For COVID-19 after a transient period, the time-series trajectory in China and Italy reached a steady regime characterized by a chaotic attractor (Mangiarotti et al. 2020). Even Jones and Strigul in 2020 propose that the chaotic properties of the system, rather than the data and the parameterization, are what impede an adequate prediction of epidemic behavior.

This study shows that the epidemic dynamics of COVID-19 has a very short prediction window at national and subnational levels, probably associated with an intrinsic chaotic behavior of its dynamics. The present analyzes also show that the incidence moments modeling approach, in conjunction with other approaches, could be useful as a simple tool for decision makers to short-term monitor the course of an epidemic in large populations that can be easily implemented across different settings. This approach could also be used to follow the course of the epidemic to prevent the possible load on the hospital system, the possible increases in ICU occupancy and the mortality that follows later, probably increasing the time window by 10–15 days.