Applications of time series analysis in epidemiology: Literature review and our experience during COVID-19 pandemic

Abstract

Time series analysis is a valuable tool in epidemiology that complements the classical epidemiological models in two different ways: Prediction and forecast. Prediction is related to explaining past and current data based on various internal and external influences that may or may not have a causative role. Forecasting is an exploration of the possible future values based on the predictive ability of the model and hypothesized future values of the external and/or internal influences. The time series analysis approach has the advantage of being easier to use (in the cases of more straightforward and linear models such as Auto-Regressive Integrated Moving Average). Still, it is limited in forecasting time, unlike the classical models such as Susceptible-Exposed-Infectious-Removed. Its applicability in forecasting comes from its better accuracy for short-term prediction. In its basic form, it does not assume much theoretical knowledge of the mechanisms of spreading and mutating pathogens or the reaction of people and regulatory structures (governments, companies, etc.). Instead, it estimates from the data directly. Its predictive ability allows testing hypotheses for different factors that positively or negatively contribute to the pandemic spread; be it school closures, emerging variants, etc. It can be used in mortality or hospital risk estimation from new cases, seroprevalence studies, assessing properties of emerging variants, and estimating excess mortality and its relationship with a pandemic.

Key Words: Time series analysis, Epidemiology, COVID-19, Pandemic, Auto-regressive integrated moving average, Excess mortality, Seroprevalence

Core tip: Time-series analysis allows us to do easily and, in less time, precise short-term forecasting in novel pandemics by estimating directly from data. These models do not need extensive knowledge of pandemic mechanisms and interactions between peoples, societal structures, and pathogens. Its secondary but equally important role is distinguishing factors contributing to the spread or slowing it down. Of course, the time series analysis approach cannot give a forecast for an end of a pandemic, nor the precise moment of its peak, but it is invaluable for fast response based on sound statistical methodology.

INTRODUCTION

Time series analysis studies are consecutive collections of observations in time to predict or forecast behavior[1]. Prediction is related to studying the possible relationships between variables or factors that influence each other or are correlated. Forecasting predicts future values based on past values of the variable and possible future values of other variables or factors, towards which we regress[1]. Suppose we do not know the future values of the variables that influence our variable. In that case, we investigate different scenarios or the structure “IF–THEN”. If an external factor takes value A, we forecast the variable value B.

Since correlation cannot be estimated directly for nonstationary processes via the standard regression techniques, the so-called spurious correlations[2], the existing number of time series models such as Box-Jenkins оr Auto-Regressive Integrated Moving Average (ARIMA) models that can deal with nonstationarity. Why do we use them in epidemiology? First, these are linear models that are simpler and easier to use than the classical nonlinear epidemiological models such as Susceptible-Infectious-Removed (SIR), Susceptible-Exposed-Infectious-Removed (SEIR), etc., for which no closed-form exact analytical solutions exist and need numerical simulations or special techniques for approximation for the long-term behavior of the model[3].

However, the difficulty in predicting cases and deaths during the coronavirus disease 2019 (COVID-19) pandemic was also emphasized by Roda et al[4]. They raised concerns about utilizing the affirmed case information as nonidentifiability in model alignments[4].

Talking about the classical nonlinear epidemiological models for epidemics, SIR is able to propose a simple model based on the two-reaction mechanism. In this way, the conditions for epidemic development, the course of a simple closed epidemic, as well as the mitigation strategies could be explained[5]. SIR has been used successfully to estimate the number of cases and deaths in outbreaks such as influenza H1N1 (2009-2010) and Ebola (2014-2016) viruses, examining the early growth, including with modifications, such as SIR with reactive behavioral changes and SIR with inhomogeneous mixing.

The SEIR model was also used mostly for influenza epidemics. Zhan et al[6] additionally used the COVID-19 historical data of 367 Chinese cities to create the transmission mechanisms and contact topology utilizing a set of profile codes. Then the method was applied to South Korea, Italy, and Iran, to predict the infection peaks before the end of March 2020[7].

However, when comparing the SIR and the ARIMA models, it was shown that ARIMA outperformed the SIR in predicting the cases of COVID-19[8].

Abolmaali and Shirzaei[7] were among the first to compare multiple epidemiological methods that can be used during the COVID-19 pandemic to monitor and even prevent the spread of infection. The authors demonstrated that predictive models, such as SIR, SEIR, ARIMA, etc., have proven useful and effective in predicting the incidence of infections. However, it was shown that SIR could not provide helpful early prediction in cases of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, giving a significant error in estimating for countries such as Brazil, United States, India and Russia. Moreover, SIR was useful mainly for predictions in the short term. On the contrary, SEIR showed better results than SIR in the long term in forecasting COVID-19[9].

Regarding ARIMA, it was shown that this model requires further information for a more precise point-by-point forecast, but the spread of COVID-19 was forecasted ultimately. Additionally, ARIMA itself helps the data to remain stationary, resulting in modeling flexibility while capturing any changes at every stage. Last but not least, the ARIMA model managed to maintain minimum error in the forecasting. Based on all these outcomes, ARIMA outperforms all other compared models[7].

Many other research groups confirmed the effectiveness and usefulness of the ARIMA model for COVID-19 pandemic forecasting[10-12].

Additionally, the ARIMA model has three parameters: p, d and q, which define the order of its three components-autoregression, differencing and moving average. The autoregressive part captures the dependence on past values or the process dynamics-the variable is regressed on its previous or lagged values. The moving average part indicates that regression errors are not independent but are a linear combination of error terms from different earlier moments[13]. The differencing part is used to transform the process into a broad sense stationary process vs consecutive differencing of its values. The combinations of these three parts with proper orders depending on the process allow us to make linear regression (LR) that accounts for the dependence of error terms for nonstationarity and assess dependence with previous moments of time. The ARIMA model can be combined with external regressors or factors. Thus, it can be applied as regression with ARIMA errors for the process, made stationary by differencing its values. The use of LR for prediction is, therefore, straightforward.

One example is the model in our recent research[13]. We predicted some weekly deaths from COVID-19 in Bulgaria from new cases by age groups with different lags. We used the viral variants as external regressors.

Another example is the analysis of the incidence and death of AIDS and HIV in China[14]. Time series models such as ARIMA need historical data to be able to extrapolate future trends. They are beneficial for short-term forecast[15] and for prediction, even for complex dynamics, such as for COVID-19, with different policies to control the spread and different viral variants in various stages of the pandemic[16]. In addition, there is a practical advantage of time series analysis over theoretical epidemiological models, and no detailed knowledge of the causal mechanism is required; instead, the structure is inferred from real-life data. This allows fast development and deployment of models during an ongoing pandemic with limited knowledge of key parameters related to the influence of measures or the characteristics of newly emerging variants.

Based on the characteristics and effectiveness of ARIMA, we created the following objectives for our study: (1) To apply time series analysis for assessment of SARS-CoV-2 variants spread; (2) to assess the excess mortality from COVID-19 in Bulgaria; (3) to demonstrate the importance of time series analysis during the COVID-19 pandemic in radiology departments; and (4) to determine the SARS-CoV-2 seroprevalence with time series analysis.

TIME SERIES ANALYSIS FOR ASSESSMENT OF SPREAD OF SARS-CoV-2 VARIANTS

Time series regression studies have been extensively used in environmental epidemiology, especially to assess short-term associations between exposures[15]. The most commonly studied factors are air pollution, weather, and pollen, but health outcomes like mortality or disease-specific hospital admissions could also be investigated. Typically, data are available at regular time intervals for both exposure and result (e.g., daily pollution levels and daily death counts). The goal is to investigate short-term relationships between them.

In line with this, since COVID-19 has expanded globally, this results in a continuous pandemic, imposing limitations and expenditures on many governments. Therefore, anticipating the number of new cases and fatalities throughout this period can be crucial in predicting future expenses and facilities[17]. Furthermore, it is necessary to analyze the data during the pandemic to expect one or other intervention strategies, mitigations, etc. In such a way, the pandemic could be monitored, controlled, and properly managed. Many studies demonstrated that among the epidemiological models, the ARIMA model showed the desired precision in predicting the number of cases and fatalities with minimum error[10,18].

Although some researchers hesitantly avoided ARIMA for analyzing COVID-19 epidemiological data, Alabdulrazzaq et al[19] demonstrated that the ARIMA technique showed accurate and valid forecasting; especially the ARIMA best-fit model for predicting the confirmed and recovered cases of COVID-19. Despite the many dynamic aspects based on the novelty of the virus and the nature of the disease, the actual values for most of the periods were within the model prediction of a 95% confidence interval. Pearson’s correlation showed high correlations between the forecast points and the actual recorded data (r = 0.996)[19]. This confirms why ARIMA is one of the best-suited models with satisfactory results and minimum error.

Different methods in time series analyses can be used, such as the hybrid machine learning approach (using multiple simple algorithms to complement and facilitate each other) to anticipate the number of infected people and mortality rate[20]; LR (based on using regression models enabling subject-matter interpretation of the data); Least Absolute Shrinkage and Selection Operator (a model that uses over regression methods for more accurate predictions); support vector machine (using optimal hyperplane in an N-dimensional space, separating the data points in different classes); exponential smoothing (forecasting univariate time series data) to determine the affected by the virus people and the deceased cases[21]; numerical modeling to assess the effect of the population age on the mortality rate[22]; numerical modeling methods such as polynomial regression (fitting of a nonlinear relationship between the value of something and the condition mean of other); Bayesian Edge (estimating probability influenced by the belief of the likelihood of a certain outcome) and long short-term memory (having the ability to learn long term sequences of observations) to estimate the prevalence of SARS-CoV-2 infection and to predict the scale of the pandemic along with the mortality rate[23]; a deep learning system for the prediction of the COVID-19 time series[24]; mathematical model about the spread of COVID-19[25]; a stochastic model considering comorbidities and age[26]; an SIQR model made stochastic, considering the uncertainty of infection progress[27]; a fractional-order dynamical system[28]; fractional calculus and natural decomposition[29]; Caputo-Fabrizio fractional derivative[30,31]; and a nonpharmaceutical intervention approach to reduce the outbreak of COVID-19[32].

EXCESS MORTALITY FROM COVID-19 IN BULGARIA: OUR EXPERIENCE

We describe an example of application of time series analysis in epidemiology: Analysis of overall mortality and its dependence on the number of registered cases of COVID-19. In other words, we used ARIMA and tried to establish which part of the COVID-19 deaths were unregistered and which part of the overall mortality was influenced by the pandemic. In our modeling, we did not try to define and determine which part of the mortality was excess-this needs serious theoretical modeling that includes seasonal climate variations, dependence on population growth models, identification of other causes of mortality and their distribution to determine expected mortality, etc. With time series analysis, we could identify the excess part directly from data, by testing how much of the mortality variance for 2020-2021 could be explained by new COVID-19 cases. We used weekly data for deaths from the National Statistical Institute[33] and from John Hopkins University for total cases[34].

Material and methods

We added to the model categorical variable to account for the changes of dominant variants in Bulgaria during the first 2 years of the pandemic wild type, alpha and delta. The coevolutionary arms race between the acquired immunity for the survival of infections and reinfections and the mutating virus did not allow capturing of a clean and powerful connection between mortality and the disease severity caused by the virus variants because mortality was the product of interaction between the two coevolving agents.

Since the immune system adapts to decrease mortality, the mortality from previous waves also suppressed further mortality via natural selection. Thus, any positive correlation between variants and mortality in time-series analysis models carries more information than it appears purely from the coefficients and their standard estimation errors. Even the slightest edge of the variants over the immune system and the process of natural selection should be treated as significant here. This was the reason to include not only weekly COVID-19 cases with lags of 0-7 and 7-14 d (L1 and L2) but also variants as factor variables with the same lags in our optimal model.

The model, described in Tables 1 and 2 and Figures 1 and 2, resulted from selecting different lags of the two chosen variables that produced the best fit (Figure 2). Also, the model was developed with the language R. It contained differencing of order I, as suggested by the ndiffs function via different unit root tests to achieve stationarity and meaningful correlations[35].

Figure 1 Residuals of our model of mortality. ACF: Autocorrelation function.

Figure 2 Model prediction of mortality based on new coronavirus disease 2019 weekly cases. COVID-19: Coronavirus disease 2019.

Table 1 Correlations between the variables for Model I-factors contributing to the spread among children.

Variables	Cases L1	Cases L2	Variants L1	Variants L2
Cases L1	1	0.69	0.02	0.03
Cases L2	0.69	1	0.02	0.02
Variants L1	0.02	0.02	1	-0.02
Variants L2	0.03	0.02	-0.02	1

Table 2 Regression models with Auto-Regressive Integrated Moving Average (0, 1, 1) errors-Model I-overall mortality and coronavirus disease 2019.

Coefficient	Estimate	Standard error
Cases L1	0.0348	0.0103
Cases L2	0.0394	0.0105
Variants L1	106.4794	134.3759
Variants L2	96.3354	134.5631
MA1	0.2594	0.1100
R2	0.953
RMSE	188.83
Bias	-2.632
MAPE	0.152

RMSE: Root mean square error; MAPE: Mean absolute percentage error.

Results

Factor variables such as variants L1 and L2 were not different. There was a significant influence on cases with lags L1 and L2 with a coefficient ratio to standard error over 3:1 (2:1 is required as a rule) and a sufficiently small correlation coefficient between their first differences of 0.69. Variants L1 and L2 had positive values but high standard errors of estimation.

Nonetheless, we considered their appearance significant after three different variants and the enormous increase of mortality over the 2 years of nearly 25% (on average) over the mean for 2015–2019 and with a high number of officially registered cases which was 11.5% of the overall population, with polymerase chain reaction (PCR) positivity on the average 12.14% (maximum 33.75%) indicating many more unregistered cases. Even this slight edge that we detected here indicated increasing severity with variants up to and including delta. Our model explains > 95% of the variation in deaths even though there was considerable variation in the mean age of new cases during these 2 years, and the exponentially increasing hospitalization (and therefore, mortality) risk with age[36]. This is evidence that most of the mortality increase in 2020 and 2021 was due to COVID-19. The model had shallow bias and mean absolute percentage error.

Conclusions from this model. Time-series analysis can serve as a first step in studying causal connections between an epidemic and excess mortality. Models such as ARIMA can show whether two or more nonstationary processes are moving together so that we can predict one behavior from another. Our models allowed us to catch when several different processes contributed with additional time lags to a resulting process, helping us uncover the link between two or more processes invisible to the naked eye. In this case, new weekly cases from the previous 2 wk, together with the changes in viral variants (factor variables), could explain 95% of the excess mortality. This is relevant as an answer to the often-appearing question of which part of the mortality during an epidemic is hidden from the official figures.

Moreover, during a pandemic with a high hospital burden, other patients have delayed treatment and are collateral victims of the pandemic. Time-series analysis can help quantify the excess mortality caused by a pandemic. It can help answer the question: Is the excess mortality due to closures or other mitigation measures, or due to the pandemic itself (although we did not try to answer this question in our study)? Closures could be added as a factor variable via the Oxford Stringency Index[37], and in a similar fashion, viral variants were added by us. Our research tried to check if variant evolution contributed to increased deaths. Still, it cannot be conclusively shown-the standard error of the relevant parameters, variants L1 and variants L2, is not small enough for that purpose (it should be at most half of the absolute value of the coefficient). This is possibly due to the only delta variant being significantly deadlier for that period[38].

However, it is essential to mention that virus variants change over time along with the clinical and epidemiological picture, as it was discussed recently by Miteva et al[39].

RADIOLOGY DEPARTMENTS DURING SARS-CoV-2 PANDEMIC AND THE IMPORTANCE OF TIMES SERIES ANALYSIS

Medical imaging is crucial for initial diagnosis, staging, and follow-up. Therefore, the organization in departments of radiology is a significant topic because there has to be a separation of COVID-19 and possible COVID-19 patients from other patients in the hospital. Usually, it is done by arranging so-called COVID corridors in the hospital, when only COVID-19 patients are being scanned. This causes an interruption of regular hospital activity.

Usually, there is a delay in diagnosing patients with other diseases, which is a problem, especially with emergency patients. In the COVID corridor, the personnel in the radiology department are fully equipped to diagnose patients. However, the equipment was not always available, especially at the beginning of the pandemic. Moreover, it was costly which limited its use. Also, deep cleaning of the department is done after the end of the corridor, which causes even more delays in the other patients’ diagnoses and more expenses.

Forecasting the COVID-19 waves is crucial because the management of the radiology departments can be done according to it. A durable prediction model allows the departmental heads to organize the necessary COVID-19 corridors according to the expected wave. The duration and the exact hours of the corridors can be correctly adjusted, thus providing the required diagnosis of COVID-19 and non-COVID-19 patients with minimization of the delay of diagnosis for each group.

Also, the personnel shifts can be arranged according to the predicted model, providing enough X-ray technicians and radiologists. Furthermore, the necessary equipment will be provided in advance, reducing the needed time for changing clothes, and lowering expenses. The heads of the departments in which patients require diagnostic imaging can organize their work according to the COVID-19 corridor active hours, thus providing a safe and calm environment. Ambulatory patients can also arrange their examination according to the available hours.

In 2020, a genetic programming prediction model was introduced and further developed into a gene expression programming model. This model predicts the cases according to two parameters-confirmed cases and number of deaths[40]. Another prediction model used in India is ARIMA, which has value in predicting cases and shows the effect of unlocking after lockdown. The ARIMA model relies on the number of positive cases, the number of performed tests per day and the average positive percentage[41]. In the United Kingdom, weighted interval scoring was used for the prediction model, which used the data from the linear progression of 7-day cases[42]. In Chile, ARIMA (henceforth), exponential smoothing techniques, and Poisson models for time-dependent count data are used[43].

TIME SERIES ANALYSIS FOR SARS-CoV-2 SEROPREVALENCE

SARS-CoV-2 serology is used to identify previous infections, both in individuals and in populations. For this purpose, changes in antibody levels against SARS-CoV-2 are monitored to assess how they change over time and how long the protective immunity is preserved[44-46].

Many seroepidemiological studies are aimed at specific populations, such as health staff, police officers, and hospitalized patients with chronic diseases and COVID-19. Sometimes they use poorly and not well-validated laboratory methods[47] and mainly aim to study only the immunoglobulin G (IgG) response[48]. Therefore, various antibody responses began to be tested to improve the evaluation and diagnosis.

Many studies demonstrated that people with confirmed COVID-19 infection developed IgA, IgM and IgG against the S1 domain of the spike protein and nucleocapsid protein within 2 wk of symptoms[49,50]. Specific IgM antibodies are detected after 5-7 d of symptoms. After approximately 14 d, IgG begins to appear. IgA responses are detected almost simultaneously with IgM or earlier[51-53].

It was also found that the levels of antibodies correlated with the severity of the disease[54-57]. Previous research has found that with time, immunity to SARS-CoV-2 natural infection is short-lived and leads to a risk of reinfection[58-60].

А cohort study was conducted among health workers at the first SARS-CoV-2 epidemic peak in London[61]. They were tested weekly for symptoms, with RT-PCR and blood samples for 16-21 wk. Serological analysis was for IgG to the S1-domain of the spike protein and nucleocapsid protein. Asymptomatic or mild SARS-CoV-2 infection has been shown to elicit faster heterogeneous responses, and antibodies are cleared more quickly, which may affect the longevity of humoral immunity to SARS-CoV-2.

Another population-based study in Catalonia was conducted from February to November 2020. A multiplex serological test was used on 5000 participants from blood samples. Responses to 15 isotype-antigen combinations were monitored, and seroprevalence of 18.1% was found in adults and 15.3% in a simulation of the total population of Catalonia[62]. Based on the severity of the disease, immune profiles reveal that with increasing severity of infection, serum responses are more stable. The age and sex of the participants, overweight/obesity, compared with normal weight, and if they were or were not smokers when included in the study. There were no significant differences in seroprevalence between the two sexes. The seroprevalence among children was lower than in adults. Children were at lower risk of seropositivity than their parents in one family. Overweight/obese participants had higher antibody levels than those with normal weight. This is confusing because of suggesting that higher levels were adjusted for the severity of the disease[62].

A study on the seroprevalence of SARS-CoV-2 is currently being conducted in Barcelona from February 2021 to March 2022. The SeroCAP sentinel monitoring system is being used. IgG detection against SARS-CoV-2 spike protein will be performed monthly from blood samples collected from three hospitals. About 3000 samples are taken per month, and the prevalence will be assessed by age, sex and time in the three health zones in Barcelona. A complete analysis of the prevalence of SARS-CoV-2 infections will be performed, considering the demographic, social and economic factors. The correlation between seroprevalence confirmed cases of COVID-19. All measures applied so far will be studied[63].

A systematic review analyzed 47 studies from 23 countries[64]. In addition, other representative population-based studies at the national or regional level have been published[65-68]. The data showed that the SARS-CoV-2 seroprevalence in the general population varied from 0.37% to 22.1%. Biological, behavioral and social factors, including vaccine coverage, influence these percentages. We must acknowledge that the titers of protection against SARS-CoV-2 are currently unknown. However, virus-neutralizing antibodies are needed to protect and control the infection[69].

CONCLUSION

Mathematical modeling of pandemics is of vital importance for several reasons. One is to study the mechanisms of interactions between people, societal structures, and pathogens for past epidemics to develop knowledge that would help predict and control future ones. These are the classical epidemiological models, such as SIR and SEIR. Another reason is to enable us to manage a current epidemic by distinguishing productive from unproductive measurements and delivering precise short-term forecasts for the number of new cases, the number of new hospital admissions, expected deaths, etc.

This is the application of time series analysis, which relies on past values to predict future ones without extensive use of theoretical knowledge with all its uncertainties during an ongoing epidemic. It is used in mortality risk estimation, seroprevalence studies, reliable short-term forecasting for the healthcare system burden, and excess mortality estimation and analysis. It has a broad spectrum of linear and nonlinear, and single and multidimensional models. This allows one to choose ease of use vs capabilities according to different contexts where they are applied. Time series analysis supplements the classical epidemiological models in predictive and forecasting capabilities. It reinforces our decisions on how to act vs an epidemic with additional analytical approaches and results on which to step on.