COVID-19 contagion forecasting framework based on curve decomposition and evolutionary artificial neural networks: A case study in Andalusia, Spain.

Many types of research have been carried out with the aim of combating the COVID-19 pandemic since the first outbreak was detected in Wuhan, China. Anticipating the evolution of an outbreak helps to devise suitable economic, social and health care strategies to mitigate the effects of the virus. For this reason, predicting the SARS-CoV-2 transmission rate has become one of the most important and challenging problems of the past months. In this paper, we apply a two-stage mid and long-term forecasting framework to the epidemic situation in eight districts of Andalusia, Spain. First, an analytical procedure is performed iteratively to fit polynomial curves to the cumulative curve of contagions. Then, the extracted information is used for estimating the parameters and structure of an evolutionary artificial neural network with hybrid architectures (i.e., with different basis functions for the hidden nodes) while considering single and simultaneous time horizon estimations. The results obtained demonstrate that including polynomial information extracted during the training stage significantly improves the mid- and long-term estimations in seven of the eight considered districts. The increase in average accuracy (for the joint mid- and long-term horizon forecasts) is 37.61% and 35.53% when considering the single and simultaneous forecast approaches, respectively.

Introduction

Since the first outbreak detected in Wuhan, China, the new coronavirus has widely and rapidly spread around the world due to its powerful human-to-human transmission capacity (Sanche et al., 2020), leading to an exponential growth in the number of infected people in all countries. This new coronavirus, which produces the disease known as coronavirus disease 2019 (COVID-19), has put human health at risk by provoking fever, cough, and myalgia as common symptoms, and potentially leading to complications, such as acute respiratory distress syndrome in a significant percentage (Chen et al., 2020, Huang et al., 2020), among others. On the 11th of March 2020, when many countries were in an emergency situation, the World Health Organization (WHO) declared this virus a global pandemic thereby, forcing the countries to adopt prevention measures such as nationwide lockdowns, mandatory use of facial masks and human mobility controls to suppress virus transmission. Many studies have shown that reducing social interaction by applying control measures helps to mitigate the spread of contagions (Kharroubi and Saleh, 2020, Kraemer et al., 2020).

In each country, the number of outbreaks, their intensity and duration depend on a wide variety of parameters, such as the size of population, the control measures applied and even specific climatological features (Livadiotis, 2020, Malki et al., 2020). Nevertheless, the evolution of the outbreaks are similar regardless of the location where they take place. The dynamics of the pandemic are as follows: at the beginning, a small number of people became infected. Then, due to the high contagion capacity, if social interaction is not limited, the number of infections increases linearly in a short period of time. After that, when a sufficient number of people contract the virus, infections begin to rise exponentially, and the outbreak quickly spirals out of control. At this time, control measures are usually applied to bring the outbreak under control. Hence, it is of huge importance to predict the contagion rate to avoid the inflexion point when the spread velocity starts to increase exponentially. Forecasting will allow the anticipation of uncontrolled situations and thus can help in adopting health system preparedness measures to avoid hospital overcrowding and in devising social restrictions to minimize the number of infections and the economic impact (Maital & Barzani, 2020).

To devise suitable economic, social and health care strategies to mitigate the pandemic effects, anticipating the evolution of the outbreaks has thus become a crucial task. Multiple forecasting techniques have been proposed since data about COVID-19 impact (regarding infected, deceased and recovered people) started to be available. In this connection, epidemiological models (Kermack & McKendrick, 1927) are the most popular approach to estimating the evolution of infectious diseases. These models use mutually exclusive compartments or states and assign them to a population of individuals to describe the dynamics of the population. Individuals flow through the compartments according to the parameters of the model. Compartmental models have been mainly applied at the beginning of the pandemic period in different locations. For example, in Anastassopoulou, Russo, Tsakris, and Siettos (2020), a four-compartmental SIDR (Susceptible–Infectious–Recovered–Dead) model was used to estimate the most important epidemiological parameters, such as the basic reproduction number, , and the rates of infection, mortality and recovery, using the COVID-19 incidence data of Hubei (China) from January 11 to February 10, 2020. Moreover, in Hauser et al. (2020), also using data from Hubei, an age-stratified SEIR (Susceptible–Exposed– Infected–Removed) model was fitted to estimate the symptomatic case-fatality ratio (sCFR) and infection-fatality ratio (IFR) over six regions of Europe. An extension of this model, including predictions of the number of cases using different time series forecasting techniques, was presented in Katris (2021).

Alternatively, machine learning (ML) techniques have also been applied to model some threatening aspects of the COVID-19 pandemic. ML approaches have a long history in solving real-world problems in different fields, including health care, economy and natural language processing. Concerning the health care area, and particularly the new SARS-CoV-2 virus, several ML and artificial intelligence expert systems have been proposed for different purposes. On the one hand, with the aim of identifying screening and management of SARS-CoV-2 positive diagnoses, methods have been developed to augment traditional identification tools, using radio imaging technology to detect abnormalities associated to COVID-19 infections (Huang et al., 2020, Ng et al., 2020) as an alternative to conventional tests, or even as standalone methods when viral testing is not an option. In Song et al. (2021), a convolutional neural network (CNN) was applied to distinguish between COVID-19 infected patients, pneumonia infected patients and healthy patients using computed tomography (CT) images. A deep neural network called COVID-Net is presented in Wang, Lin and Wong (2020) for identifying positive SARS-CoV-2 diagnoses using chest X ray (CXR) images. In Abbas, Abdelsamea, and Gaber (2021) a class decomposition (Abbas, Abdelsamea, & Gaber, 2020) technique is applied in combination with a pretrained CNN, improving the performance of the classifiers when the class decomposition layer is included as a preprocessing step. A combination of deep learning and classical ML algorithms is proposed in Sethy and Behera (2020), where deep features of CXR images are extracted from the fully connected layer of a CNN and used to feed a support vector machine for distinguishing between healthy, COVID positive and pneumonia images. In Tamal et al. (2021), a set of radiomics features from CXR images were selected and used to train three classical ML classifiers, providing an accurate, fast and automatic method that can be integrated with standard X-ray reporting systems. Most recently, in Garg, Salehi, Rocca, Garner, and Duncan (2022), the performances of 20 different CNNs trained for classifying patients into three and two classes using chest CT images achieved an accurate and very efficient classification model.

Apart from the new imaging diagnostic assistance mechanisms, multiple algorithms and statistical techniques have been applied to obtain an accurate prognosis of the pandemic rates, the peak of outbreaks in different countries or even estimations of specific pandemic wave scopes and to create analytical models that could act as decision support systems. In this sense, in Benvenuto, Giovanetti, Vassallo, Angeletti, and Ciccozzi (2020), a classical autoregressive integrated moving average (ARIMA) is applied to analyze the trend of COVID-19 prevalence and incidence and to perform short-term prediction about the spread of the virus. In Wang, Zheng, Li and Zhu (2020), three critical points of infections and recovered cases are estimated in different countries by considering a hybrid forecast model based on a logistic curve fit and Prophet (Taylor & Letham, 2018) application. A comparison of six time series forecasting techniques applied to active COVID-19 cases is performed in Papastefanopoulos, Linardatos, and Kotsiantis (2020). Following a similar methodology, in Ribeiro, da Silva, Mariani, and dos Santos Coelho (2020), the predictive capacity of different machine learning regression models is measured using data from 10 Brazilian states. More complex models are used in Verma, Mandal, and Gupta (2022) for forecasting purposes, where long short-term memory (LSTM) recurrent neural networks are designed for predicting the contagion rate in 4 Indian states. Multiple regressors were applied in An et al. (2020) to predict the risk of mortality using different characteristics of an infected person, concluding that variables such as advanced age or taking metformin are important predictors that influence the output probability. An assisted fuzzy case-based reasoning (FCBR) algorithm for determining patient attention priorities based on eight factors is presented in Geetha, Narayanamoorthy, Manirathinam, and Kang (2021) with the aim of improving medical assistance and reducing the mortality of COVID-19 patients. In Fidan and Yuksel (2022), an unsupervised study is carried out to analyze the importance of city-related parameters, such as population and environmental variables, in addition to the number of cases to establish restrictions to contain the spread of the virus. In Desai (2021), the effects of news sentiments are included in the forecast task, concluding that negative news sentiments could help to reduce the contagion rate. Finally, in Khan et al. (2021), a very comprehensive review of research regarding the use of diverse ML algorithms to combat the COVID-19 pandemic is presented, where applications such as diagnosis, detection or forecasting are carried out using different types of data in several countries.

Among all the state-of-the-art techniques in ML, there is one option that stands out: feed-forward neural networks (Bishop et al., 1995). This is a highly adaptative technique that can be used to model most nonlinear problems, mainly due to its universal approximation capability (Hornik, Stinchcombe, & White, 1989). These models are layer-divided structures, where each layer is composed of computational units that are connected to the nodes of the next layer by weighted connections. Artificial neural networks (ANNs) have been trained to solve multiple classification and regression problems (Abiodun et al., 2018). In this sense, backpropagation is the most popular algorithm for adjusting the network parameters in the training stage. This method computes the gradient of the loss function with respect to the weights of the network connections. However, even though this approach can efficiently adjust the weight magnitudes, a global optimal solution is not guaranteed. The reason behind this is the large extent of the multidimensional error surface and the presence of local optima where the backpropagation can get stuck. To perform a more exhaustive search over this surface, evolutionary algorithms have been presented as an alternative to discover and explore different regions (trying to escape from local optima) where more accurate solutions could be located. For forecasting purposes, evolutionary artificial neural networks (EANNs) have been applied in different fields such as renewable energy (Gómez-Orellana et al., 2022, Mason et al., 2018), economy (Au et al., 2008, Chiroma et al., 2015, Yi-Hui, 2007) or environment (Lu, Fan, & Lo, 2003). When modeling a future value as a target variable, it is important to consider the existing trade-off between accuracy and anticipation. Short-term events are usually easier to model but do not offer enough anticipation time, while long-term horizon forecasts are very useful for anticipating important events, but their accuracy decreases when the forecasting horizon increases. It is worth mentioning that the concepts of short and long horizons are relative and depend on the field of study where the prediction is carried out. Thus, in Gómez-Orellana et al. (2022), a 48-h forecast is considered a long-term horizon, while in Au et al. (2008) they refer to weekly events as short-term events. In the field of epidemiology, since the parameters associated with an outbreak evolve rapidly and the information is reported daily, 1, 3, and 5 days can be considered to be short-, mid- and long-term horizons, respectively.

In this study, we make the following contributions. First, we propose a novel methodology for the analysis of the contagion curve by extracting evolution characteristics. The extraction procedure is performed through an iterative process of polynomial model fitting. The coefficients of the polynomial model fitting the cumulative contagion curve describe the shape and evolution of the curve, where points represent days since the beginning of each wave. This procedure allows the collection of inherent features of a process that can be used as related information, which is especially useful when process exogenous related information is scarce. Second, we apply this methodology to real data involving the number of infected people in different locations of Andalusia (Spain) to build a transfer function per zone that can predict the evolution of the contagion rate on the different stages of the pandemic. By means of this methodology, we build a multivariate time series dataset per sanitary district by considering the real cumulative contagion curve data as the dependent series and the polynomial coefficient series as independent correlated series. Third, we estimate the parameters of these transfer functions using EANNs, with data from different periods of the pandemic waves that took place in Andalusia. Given that every district presents different contagion rates in the considered outbreaks, different architectures of EANNs are considered using distinct basis functions in the nodes belonging to the hidden layer with the goal of finding the best network scheme for each district. By using distinct periods for the estimation of the parameters, we build a forecast model per sanitary district using a different order of lags that serves to obtain an accurate prediction about the contagion rate in different stages of an outbreak. This results in models that can identify the beginning and ending of new outbreaks and predict the evolution of the contagion rate during the most critical phase. The parameter estimation is carried out by considering mid- and long-term time forecast horizons ( and days, respectively). The reason behind the choice of these time horizons is that they offer an acceptable trade-off between anticipation and accuracy. To serve as decision support for adopting prevention measures, 3- and 5-day predictions of the number of contagions are enough to devise suitable containment schemes. Larger horizons lead to less accurate results. It is not appropriate to base sensitive decisions such as the application of restriction measures on such results. Additionally, we consider a combined estimation where both time forecast horizons are involved in a single model using a multitask evolutionary artificial neural networks (MuEANNs) approach. The resulting models can act as a decision support system in new emergency situations caused by new viruses or strains with similar contagion rate behavior. The results obtained with the models that are trained with the extracted polynomial information significantly outperform those obtained with models trained with pure autoregressive data in almost all districts. Specially, the accuracy is increased by up to 73.55% in the case of the Almería district.

The remainder of this paper is organized as follows: in Section 2, we describe the data source and perform an analysis of the data used in this study. The proposed methodology is explained in Section 3 where the polynomial curve fit, the dataset building process, the autoregressive forecast models and the EANNs are presented. The experimental design and results obtained for each architecture, dataset and forecasting approach are presented in Section 4. Finally, we conclude the paper in Section 5.

Data description

The data used in this study and some considerations are specified in this section. This article focuses on eight areas of Andalusia, the country’s second largest and most populated autonomous region, located in southern Spain. The data were obtained from the official website of the Andalusia government,1 where the number of COVID-19 diagnosed, cured and deceased people were reported daily and segregated into sanitary districts. Sanitary districts are administrative divisions that have local health management competencies in specific zones of Andalusia. These divisions are represented in Fig. 1, where the eight districts that include the provincial capitals are highlighted in green and their populations are specified. For this study, the daily reported information in these eight sanitary districts from July , , to August , , is used, resulting in a total of observations, each representing one day of an outbreak. During this period of time, four different waves took place in Andalusia. Although the data reported by the Andalusia government about positive diagnosed contagions are available since March , , the data until July , correspond to the first wave of the pandemic. During this wave, the country wide lockdown resulted in a low intensity and short time wave. Hence we decided to exclude this period from the analysis.

Fig. 1

Geography of the sanitary districts of Andalusia. The districts including provincial capitals are highlighted in green.

The reported information presents a weekly pattern of fewer diagnosed contagions during weekends, making the positive diagnoses time series highly noisy. This effect is common to most countries (Ricon-Becker, Tarrasch, Blinder, & Ben-Eliyahu, 2020) and is due to the lower number of patients tested on Saturdays and Sundays, differences in testing timings and reporting delays. Missing data on weekends are usually included in the first days of the weekly reports thereby preserving the real number of positive diagnoses. Fig. 2, Fig. 2 show, respectively, the daily and cumulative positive COVID-19 diagnosis time series belonging to the district of Córdoba, one of the districts considered in this study, from July , to August , .

Fig. 2

Daily (a) and cumulative (b) reported positive COVID-19 diagnoses in the district of Córdoba from July , , to August , .

In Fig. 2(a), the noise produced by the weekly cycles is observable, resulting in a visible sawtooth effect in the contagion time series. The representation of the cumulative case time series is shown in Fig. 2(b) where the noise is considerably mitigated, leading to a soft and monotonically nondecreasing time series. The temporal trend of this district is repeated for the other districts, while varying in magnitude according to specific district parameters, such as the geographic location or the total population. The temporal positive diagnoses returned four different waves, each one defined by an increment and followed by a drop in the contagion rate.

As mentioned above, the first wave has not been considered due to its low intensity and short wave time. Therefore, for this analysis, the following waves have been used: the second wave spans from July , , to December , , and its long-term evolution is characterized by a high peak of positive contagions and a gradual increment and drop in the contagion rate. During this wave, the average age of the people needing hospital care decreased compared to the first wave and the mean hospitalization time was also reduced (Iftimie et al., 2021). The third wave spans from December , , to March , . The beginning of the contagion rate of this wave is highly accelerated because of the removal of control measures due to the Christmas holidays. During this period, mobility and social interaction increased, resulting in a quick increment in the number of infected people and the maximum number of COVID-19 cases was reached. At the end of the holiday period, when control measures were reset, the contagion rate slowed down. The fourth wave, lasting from March , , to June , , is defined by a soft initial diagnoses increment owing to the beginning of the vaccination campaign. This wave presents two peaks provoked by the different changes in the containment measures policies. Finally, the fifth wave took place in summer, from June , , to August , , a period where control policies were loose. However, the high vaccination rate produced a wave with weak magnitude and short time period. Considering these partitions, the second, third, fourth and fifth waves consist of , , and observations, respectively. These dates have been chosen based on expert knowledge and the behavior observed in the daily contagion curves of the considered districts.

The nature of these four waves is very different in terms of evolution because of the distinct social policies, control measures and vaccination campaigns that took place along their respective developments, resulting in a heterogeneous set of rates for each wave. For this reason, with the main objective of modeling the complete period of the pandemic, including the most diverse information possible about all stages is crucial. In this study, different periods for each wave are considered: beginning, growth and stabilization. Including information about these three stages in the parameter optimization phase results in models that can detect imminent increases in contagion rates that can potentially lead to local outbreaks, accurately predict the evolution of an initiated outbreak and recognize in advance the end of virus transmission.

Proposed methodology

In this section, the methodology proposed in this paper and the process to generate the datasets used in the experimentation are detailed. To provide an overview and facilitate understanding, Fig. 3 presents a flowchart containing all the steps detailed in this section.

Fig. 3

Overview flowchart describing the applied methodology.

Curve fitting

Many empirical data obtained from natural processes describe patterns that may be defined using algebraic expressions. Curve fitting techniques help to understand the inherent behavior of a process and may be useful in description, clustering or forecasting tasks (Motulsky & Ransnas, 1987). If the nature of the underlying process is known, an accurate parametric function working as an approximator results in a soft curve whose parameters constitute a characterization of the empirical data evolution. The information extracted from the fitted curve can be used in ML tasks as important features, both in supervised (Hamidi, Ghassemian, & Imani, 2018) and unsupervised (Abraham et al., 2003, Martin et al., 1998) contexts.

In this paper, we consider the approximation by a polynomial function of the cumulative curve of contagions of the districts involved in this study by the least squares method. The shape of the curve described in each wave is defined by a soft initial growth. Then, when a sufficiently high number of infected people is reached, an exponential increment occurs due to the high contagion rate. The waves finish with a linear evolution when the outbreak is under control, but residual positive cases are still present. According to this shape, a third-order polynomial is the most interpretable and lowest-variance model that can accurately fit the cumulative curve of each wave individually. The polynomial regression is a type of fit where the dependent variable linearly depends on the powers of a single independent variable, in this case, the number of days from the beginning of the wave. A -degree polynomial model, which is composed of parameters, is defined as: where is the cumulative number of contagions after days from the beginning of an outbreak, , are the parameters to be estimated and is the error term. For this study, the chosen -degree polynomial model is composed of parameters to be estimated during fit and has a single inflexion point along with its domain, where the function changes its concavity.

Following an iterative procedure, for each time period between the start of the outbreak (Day ) and a day , a -degree polynomial model is fitted. Once the polynomial model has been fitted, the estimated parameters are used to describe a point of the curve, which is represented as: , where includes the four parameters of the third order polynomial. Note that each point of the curve represents a single day belonging to the th outbreak, where . Moreover, it is worth mentioning that the sanitary districts are considered independently. A more detailed description of the variables representing each point is given in Table 1.

Table 1

Descriptors of a point of the cumulative curve of contagions.

Variable	Description
α_1,d,i	Estimated α_1 after d days of the ith wave
α_2,d,i	Estimated α_2 after d days of the ith wave
α_3,d,i	Estimated α_3 after d days of the ith wave
α_4,d,i	Estimated α_4 after d days of the ith wave
y_d,i	Cum. contagions after d days of the ith wave

Given that for building accurate models a minimum number of days is needed, can only take values in the range . In this way, is the minimum number of observations needed to build accurate models, and is the total duration, in days, of the th outbreak. Therefore, the number of polynomial models to be fitted for outbreak is .

For the estimation of , the evolution of the quality of the fitted polynomial models has been measured for the sanitary districts under study, using the coefficient of determination as an evaluation metric, defined in Eq. (2) for the th wave. statistic represents the proportion of the variance of the dependent variable explained by the independent variables: where , a real number, is the cumulative contagions after days, is the number of cumulative contagions predicted by the model, and is the mean of observed cumulative contagions. This metric ranges from to , reaching its maximum value when the model is able to entirely explain the variance present in the real data, e.g., fitting perfectly to the empirical nature of the observed contagion rate.

Analyzing the evolution of the mean magnitude of the fitted polynomial models of the eight districts, the minimum feasible choice for each outbreak is based on the instant from which it begins to stabilize. At the beginning of the waves, the cumulative number of contagions is reset to , considering the waves individually, e.g., ignoring the cumulative number of infected people of past waves. Fig. 4 shows the average, maximum and minimum evolution of the coefficient of determination for the eight sanitary districts, and the chosen is indicated with a vertical line. The minimum number of days where the determination coefficient variance begins to stabilize corresponds to , , and for the second, third, fourth and fifth waves, respectively. Considering these initial thresholds for fitting polynomial curves, from the observations of each district, a total of days are used as the minimum number of observations to build accurate regressions, resulting in a total of polynomial fits on the entire pandemic period for each district. In Table 2, the , and standard deviation, , of all approximated models for each district are presented. As shown in Table 2, the selected polynomial models fit analytically well for all districts and waves.

Fig. 4

Mean, maximum and minimum evolution of the eight sanitary districts for each wave.

Table 2

, and standard deviation, (), of the statistics for the fitted polynomial models.

District	R^2
	Mean	SD
Almería	0.9954	0.0043
Cádiz Bay	0.9961	0.0023
Córdoba	0.9944	0.0075
Granada	0.9947	0.0047
Huelva Coast	0.9932	0.0044
Jaén	0.9929	0.0166
Málaga	0.9959	0.0029
Sevilla	0.9967	0.0028

Autoregressive models

Formally, a time series of length is an ordered sequence of values . The timestamps conform to a sequence of positive and ascending discrete values, which, in most cases, are equally spaced. Time series are produced by any kind of sequenced phenomenon, whether natural, such as the climate in a particular location, or artificial, such as stock prices of a given enterprise over time. Due to its intrinsic temporal nature, time series have been the object of investigation in multiple fields of data mining, such as forecasting, classification, clustering, and outlier analysis (Fu, 2011). A time series may be composed of more than one time-dependent variable measured at the same time instant. In this case, time series are called multivariate, and its analysis involves modeling the inherent relationship between the dimensions that compose it. Mathematically, a multivariate time series of length , is composed of dimensions; thus, the value is a vector of length . In this sense, the point descriptors, obtained from the polynomial fit iterative process described in Section 3.1 and detailed in Table 1, are used to build a multivariate time series dataset per district.

Within the scope of time series forecasting, autoregressive (AR) models are a useful technique for modeling future values by using the information of the past instants. A -order AR model, denoted as , can be applied to the target variable as follows: where is the target value obtained from past lagged values ( for ), is a constant that serves as the intercept of the model, is the th coefficient of the model, and is the error term. AR models represent any value of the series as a combination of its past values. Specifically, in this study, AR models will be used for the prediction of the number of cumulative contagions of the th outbreak of ; hence, (note that is the number of days since the beginning of the th outbreak). Furthermore, when the forecasting of the number of cumulative contagions is carried out by using multivariate time series, a vector autoregressive model (VAR) (Zivot & Wang, 2003) is used. These VAR models generalize the AR model by including information of the different dimensions of the time series (independent terms) to model the target time series (dependent term). models include lags of all the dimensions of a given time series, in our case, lags of the point descriptors presented in Table 1. Consequently, the model composed of dimensions is defined as follows: where is the vector of model coefficients and is the vector containing the th lagged representation of the curve point belonging to the th outbreak, i.e, the th lagged values of the variables described in Table 1. Note that the vector has length . In the same way that happens with AR models, VAR models will be used for the prediction of the number of cumulative contagions of the th outbreak of ; hence, .

Artificial neural networks

Artificial neural network (ANN) models (Bishop et al., 1995) have been applied to a large number of regression and classification tasks, demonstrating great performance due to their capability of being universal approximators (Hornik et al., 1989). In this sense, multiple ANN models have been proposed in recent years: one of the first approaches was the multilayer perceptron (MLP) (Bishop et al., 1995) with sigmoid units (SUs) as basis functions for the nodes in the hidden layers. Other alternatives include the radial basis function (RBF) (Broomhead and Lowe, 1988, Poggio and Girosi, 1989) neural networks, which usually make use of Gaussian transfer functions in the nodes of the hidden layers, and multiplicative networks, which use product units (PUs) (Durbin & Rumelhart, 1989), computing a weighted product instead of a weighted sum for a node input.

In this paper, given that mid- and long-term forecasts are considered, the use of single and simultaneous forecasts is proposed. In this sense, two different models of ANNs are developed. On the one hand, ANNs aim to perform each forecast task separately (monotask), and on the other hand, ANNs achieve both forecast tasks simultaneously (multitask). Accordingly, monotask models only have one output, whereas multitask models have two outputs. Considering that the aim of this work is to forecast the cumulative number of contagions, two different horizons have been considered: and . These two forecast horizons are focused on predicting the cumulative number of contagions at mid- and long-term, i.e.,  and days ahead, respectively.

Therefore, the ANNs models applied to perform the forecast tasks can be defined as follows: where is the model output for the forecast task considering the output , is a vector including the basis function parameters corresponding to the synaptic weights: including the weights from the hidden layer to the th output node, and , containing the weights of the connections from the input nodes to the th hidden node, being the number of neurons in the hidden layer. represents the basis function of the th hidden neuron, being the set of inputs of the ANN model (Table 3 shows a complete description of the different variables used).

Table 3

Inputs () included in the EANNs for day belonging to the th outbreak. is the total number of inputs for the generated datasets. Note that the different datasets are named according to the AR models used (either AR or VAR).

Lags	Dataset name	x_d,i	I
p=1	AR(1)_y	{y_d−1,i,d_i}	2
	VAR(1)_y,α_4	{y_d−1,i,α_4,d−1,i,d_i}	3
	VAR(1)_y,α_3,α_4	{y_d−1,i,α_3,d−1,i,α_4,d−1,i,d_i}	4
p=2	AR(2)_y	{y_d−1,i,y_d−2,i,d_i}	3
	VAR(2)_y,α_4	{y_d−1,i,α_4,d−1,i,y_d−2,i,α_4,d−2,i,d_i}	5
	VAR(2)_y,α_3,α_4	{y_d−1,i,α_4,d−1,i,α_3,d−1,i,y_d−2,i,α_4,d−2,i,α_3,d−2,i,d_i}	7

With respect to the type of neurons of the hidden layer, in this study, two different basis functions (BFs) are considered:

Radial Basis Function (RBF), as a kernel BF, with a Gaussian Transfer Function (GTF), defined as follows: where and being the centroid and radius of the GTF, respectively.

Product Unit (PU), as a projection BF, defined as follows: where are the weights from the input layer to the hidden layer.

Furthermore, combining BFs in the hidden layer has some advantages, such as providing flexibility to the decision rules. Any continuous function can be decomposed into two different types of functions: one belonging to the projection group and the other belonging to the kernel group (Donoho & Johnstone, 1989). In this sense, hybrid ANNs (both monotask and multitask) have also been considered. Such hybrid models combine PU (projection) and RBF (kernel) BFs in the hidden layer, i.e., they perform a linear combination of both types of BFs. These hybrid ANN models are defined as follows: where and are the numbers of hidden neurons of the first and second types, respectively. contains the coefficients of the ANN model, being the coefficients between the hidden layer and the th output node, and and being the weights connecting the input layer to the th hidden neuron of the first and the second type, respectively. and represent the basis functions of each type, defined in Eqs. (6), (7).

The general forecast task proposed in this study considers the real number of cumulative contagions as the continuous variable to be forecasted. Consequently, for the evaluation of the ANNs, the mean square error (MSE) of the network output layer with respect to the real values is used as the loss function, defined as follows for the monotask models: where and are the real and forecasted cumulative numbers of contagions, respectively. Note that specifies the general output of the monotask model () forecasting the cumulative number of contagions for the day of the th outbreak.

Moreover, as mentioned above, in this study, two forecast horizons are considered: or . Accordingly, the loss function to evaluate the multitask models is defined as: where, in this case, and .

Evolutionary artificial neural networks

Even though ANNs are considered to be universal approximators, for the different BFs, the required training time and final complexities of the networks may vary significantly. The speed of convergence may be a problem for those cases involving a sufficiently high number of patterns in the training set or when the complexity of the network increases, becoming the main challenge to be addressed (Livni, Shalev-Shwartz, & Shamir, 2014). In most studies, gradient-based methods have been applied to optimize ANN parameters. However, due to the nonlinear nature of ANNs, the optimization process may converge to one of the many local optima existing on the error surface of these models (Sutton, 1986).

To perform a more exhaustive exploration process over the error surface, Yao (1993) presented a thorough review of evolutionary algorithms (EAs) applied to the optimization of ANN parameters. EAs use a set of different candidate solutions (population) to: (1) simulate the evolutionary process and natural selection and (2) provide the algorithm with a global search capability with the objective of discovering regions on the error surface where good performance solutions are located. Therefore, EAs represent an efficient method to optimize both the structure and weight connections of ANNs.

In this work, the EA proposed in Martínez-Estudillo, Martínez-Estudillo, Hervás-Martínez, and García-Pedrajas (2006) is applied to optimize the ANN models described in Section 3.3. In Algorithm 1 the pseudocode of the EA is shown. As the goal of the EA is to minimize the MSE of the ANNs, the fitness function guiding the evolutionary process is defined as a strictly decreasing transformation of the MSE:

Consequently, the fitness for the monotask or multitask ANNs is calculated using the MSE defined in Eq. (9) or Eq. (10), respectively.

The EA begins creating an initial population of ANNs. Each ANN (individual) of the population is randomly generated, i.e., the number of neurons in the hidden layer, the number of connections of each hidden neuron (to link both the input layer to the hidden layer and the hidden layer to the output layer), and the weights of those connections are randomly generated considering the parameters described in Table 4 of Section 4.2. Next, the ANNs are evaluated and ranked according to their fitness.

Table 4

Parameter values that have been used in the EA for all models (PU, RBF and RBFPU for both MoEANNs and MuEANNs).

Parameter	Value
Independent runs	40
Stopping criteria:
(1) maximum number of generations	1500
(2) consecutive generations without improving individuals	10
Population size	1000
Number of hidden layers of each individual	1
Minimum number of hidden neurons (initialisation)	2
Maximum number of hidden neurons (initialisation)	3
Maximum number of hidden neurons (whole process)	4
Range of hidden neurons to be added or deleted	[1,3]
Range of links to be added or deleted	[1,5]
Range for weights between input and hidden layer	[0.1,0.9]
Range for weights between hidden and output layer	[−10,10]

After that, ANN optimization (evolutionary process) is performed as long as the stopping criteria are not met. Specifically, two stopping criteria are considered, as shown in Table 4 of the following Section 4.2. Hence, if either of them is reached, then the EA stops. In each generation of the evolutionary process, the worst 10% of ANNs are replaced with a copy of the best 10% of ANNs, which represents an elitist pressure since the best 10% of individuals (after being cloned, they constitute the 20% of the population) are optimized in a different way, as described below. Next, all individuals are evolved by simultaneously applying two types of mutations: parametric and structural.

Parametric mutation is applied to the best 10% of ANNs and updates their weights by adding Gaussian noise, whose variance decreases throughout the evolutionary process, i.e., the strength of the changes decreases as the individuals become better. In that way, this adaptive variance (particularly to each individual) dynamically modifies the intensity of the exploration of the error surface, favoring exploitation as the performance of the individual increases.

On the other hand, structural mutation (applied to the remaining 90% of ANNs) alters the structure of individuals. More specifically, this type of mutation modifies the number of hidden neurons and their number of connections between input and output layers, i.e., it explores a wider area of the search space while trying to maintain the diversity of individuals. In particular, the EA applies five types of structural mutations: Add Neuron, Delete Neuron, Add Link, Delete Link and Neuron Fusion (Gutiérrez, Hervás, Carbonero, & Fernández, 2009), which are applied sequentially to each individual.

Thus, the performance of individuals is improved throughout the evolutionary process, generation after generation, maximizing their fitness (and hence minimizing their MSE). Finally, the ANN with the best fitness is selected as the final solution when the EA stops. Since the population is randomly initialized, the EA has a stochastic component. Therefore, different runs will lead to distinct populations and, consequently, to different final solutions.

In this way, and regardless of the BFs used in the hidden layer, the developed evolutionary artificial neural networks (EANNs) will be denominated as monotask EANNs (MoEANNs) and multitask EANNs (MuEANNs) depending on whether they perform the forecast tasks separately or simultaneously, respectively.

Experiment settings and results

The settings used in the experimentation and the results obtained using EANN models are presented in this section. In addition, a statistical analysis is performed to evaluate significant differences between the models and datasets used, with respect to the distinct sanitary districts considered.

Datasets

To analyze the effect of the polynomial information in the forecasting task, for each district, different multivariate time series datasets are constructed by combining different sets of input variables and autoregressive orders. Experimental results have empirically demonstrated that using the four polynomial coefficients simultaneously generates overfitted models, performing poorly on the generalization set. For this reason, only the coefficients defining the critical point of the -degree polynomial model are considered in the dataset building process. The critical point is produced when the function changes its concavity which is related to the instant where the outbreak gets out of control, and, therefore, the number of positive cases begins to rise exponentially. This moment is reached when the second derivative of the polynomial is equal to zero, i.e., . Consequently, only these two coefficients will be considered as polynomial information in the datasets because they characterize the inflexion moment with having more importance.

Table 3 shows the inputs of the EANN models. Regarding the coefficients extracted from the polynomial approximation, except for the coefficient , which is always included in the AR/VAR models, the other two coefficients, and , can be included or not. Moreover, apart from these coefficients, the number of days from the beginning of the th wave is always included as an input variable with the aim of temporally locating the fitted evolution. Notably, the and models are applied with . Therefore, the simplest dataset is built with only two input variables ( lagged time and ), being identified as . On the other hand, the most complex dataset is built with seven input variables (being , and lagged two times and ), which is identified as . Note that the total number of input variables is specified in column of Table 3.

Regarding the target or output variables, as mentioned in the previous Section 3.3, the cumulative number of contagions at two different time horizons is employed, and . For the MoEANN approaches, two different datasets are built as , depending on the forecasting horizon, or . On the other hand, for the simultaneous forecasting with and carried out by the MuEANN approach, the dataset is built as . Note that is one of the six combinations presented in Table 3.

The train and test partitions have been generated while taking into account the three stages previously mentioned in Section 2. From each wave, a specific period has been included in the generalization test, aiming to include as many different cases as possible to build robust models. The rest of the stages are included in the training set. Specifically, the testing partition is composed of the beginning stage belonging to the second wave, the growth stage of the fourth wave and the stabilization stage from the third wave and fifth wave. This splitting methodology is shown in Fig. 5, where periods for training and testing sets are represented over the cumulative contagion curve and the evolution of polynomial coefficients ( and ) for the district of Sevilla. This hold-out splitting methodology responds to the need to involve as much information as possible about the different stages of all the waves in the training set. With these nonrandom partitions, we aim to create a heterogeneous and representative set of observations for training and testing sets, given that the application of cross-validation would have increased unaffordably the needed runs of the algorithm. Therefore, from the patterns that were involved in the datasets for each district, the final datasets used in the experimentation contain and observations for the two considered autoregressive orders and , respectively. Note that for , the first observation from each of the waves is needed to build the final datasets, whereas in the case of , the initial observations needed are . From these final datasets, and for , (and and for ) patterns are used for building the training and testing sets, respectively (representing 80% and 20% of the total data).

Fig. 5

Train and test partitions of the district of Sevilla.

Experimental settings

The contagion curve is specific for every district, varying in magnitude and transmission rate according to geographical and population parameters. To find the best model and dataset for each specific district, the experimentation carried out in this paper makes use of all generated options detailed in Table 3 for all districts. In this way, the considered BFs for the ANNs detailed in Section 3.3 have been applied to estimate the weights from the input layer to the hidden layer of the models, resulting in different combinations per district: ANNs (using each one a different BF combination: PU, RBF or the hybrid one RBFPU) applied to the different datasets generated from the AR/VAR models. In addition, depending on whether the two forecast tasks are approached simultaneously or not, the two different approaches proposed in this work are compared: (1) MoEANN, which considers the two tasks separately ( or days ahead as different tasks), and (2) MuEANN, which groups both tasks carrying out them simultaneously. Consequently, different models ( MoEANNs and MuEANN) were considered for the combinations of each sanitary district, resulting in a total of experiments being executed for each district.

Regardless of the district, dataset and computing error perspective, the evolutionary algorithm used to optimize the EANN models is applied with the same configuration. More specifically, in Table 4, the parameters of the EA are shown, together with their values.

The weight values of the connections from the hidden layer to the output layer are randomly generated regardless of the type of BF of the model. However, the weight values of the connections from the input layer to the hidden layer are randomly generated for PU units, whereas for RBF, units the weight values are initialized using K-means (Lloyd, 1982) as a clustering method to determine the Gaussian function centroids. For hybrid EANN models combining RBFPU units, both types of BFs have the same probability of being selected in the initialization process of the individuals, and the final number of both types of units will depend on the evolutionary process.

Input features for all datasets are scaled in the interval . The output layer is composed of and linear units for the two MoEANN and the MuEANN approaches, respectively, and their output values are scaled in the interval .

The described configuration for performing the optimization of the models follows the guidelines published in Gutiérrez et al. (2009) and Martínez-Estudillo et al. (2006).

The EANN models containing SU basis functions in the hidden layer, either pure or combined with PU or RBF, were affected by a high overfitting and, consequently, performed poorly in the test set. Thus, this basis function has been excluded from this study.

Results

The results presented in this section aim to compare the model performance differences produced when polynomial coefficients are included as training features and the extent to which multitask models are able to take advantage of using the same structure and parameters for tackling several tasks simultaneously. The results are expressed in terms of root mean square error (RMSE, applying the square root to Eqs. (9), (10)). As explained in Section 3.4 the EANN models are stochastic. Thus, runs are performed and the results are expressed as the mean and standard deviation of the RMSE obtained for the generalization sets, i.e., patterns unseen during the training stage. Given the interest in knowing which BF is better and the best value for in the AR/VAR models, the results are segregated considering the different options, as well as for the individual districts.

In addition, the performances of the models for the two forecasting horizons are also compared. Note that for single horizon forecasts, only the corresponding output neuron error of the MuEANN model is considered. Regarding individual forecasting horizons, Table 5, Table 6 show the results obtained for and , respectively, in the eight districts considered. On the other hand, worth mentioning that when simultaneous time horizon forecasting tasks are approached, the results achieved by the MoEANN approaches (one tackling and the other tackling ) are the average of the errors of the two single MoEANN models. Table 7 shows the results obtained for all zones and datasets considering these two horizons simultaneously. In Table 5, Table 6, Table 7, best results per forecasting approach and district are highlighted in bold, whereas the second best result is highlighted in italics.

Table 5

Performances of the MoEANN and MuEANN models trained with the different combinations of input features and autoregressive orders, , evaluating the errors of the forecast horizon. The results are expressed as RMSE of the generalization set, and SD stands out for Standard Deviation.

District		BF		p		MoEANN				MuEANN
						AR(p)_y	VAR(p)_y,α_4	VAR(p)_y,α_4,α_3		AR(p)_y	VAR(p)_y,α_4	VAR(p)_y,α_4,α_3
Almería		PU		1		312.61_13.72	95.62_13.24	85.42_9.52		275.97_18.92	119.6723.75	108.78_16.77
				2		300.24_12.58	103.39_15.85	90.7212.90		275.81_22.80	136.72_30.20	126.75_26.41
		RBF		1		312.11_13.72	176.52_29.38	181.87_49.91		332.09_18.92	257.04_67.98	304.01_80.76
				2		293.53_12.58	184.22_46.05	203.34_61.35		353.01_22.80	255.39_84.67	301.07_72.12
		RBF+PU		1		131.90_34.75	107.12_26.05	98.69_28.41		147.53_24.31	142.82_34.69	146.51_42.91
				2		136.14_54.19	97.71_25.18	99.61_46.13		152.75_24.91	164.12_45.26	171.78_45.95
Cádiz Bay		PU		1		358.68_4.41	223.44_7.76	222.20_10.90		397.65_20.54	237.84_32.91	247.07_24.69
				2		345.91_12.52	234.57_15.14	235.53_22.00		354.00_33.69	228.06_26.10	229.83_41.18
		RBF		1		292.56_4.41	168.0735.62	173.03_33.94		295.57_20.54	268.97_91.29	281.15_68.42
				2		299.97_12.52	153.89_37.15	195.45_63.40		303.96_33.69	333.86_382.05	271.46_78.10
		RBF+PU		1		306.97_42.29	201.65_27.82	194.85_39.42		226.17_38.12	185.51_57.84	174.02_59.10
				2		295.82_58.44	210.14_38.15	191.68_40.64		217.75_49.77	184.24_67.41	175.1265.75
Córdoba		PU		1		195.73_4.50	82.58_5.87	76.00_10.42		194.04_8.90	95.40_9.57	93.70_8.34
				2		182.71_13.97	80.43_7.57	74.41_12.31		186.74_9.00	89.40_9.87	86.68_11.13
		RBF		1		189.78_4.50	155.67_32.39	166.90_33.23		226.67_8.90	196.20_33.28	209.91_40.78
				2		186.41_13.97	161.03_29.05	176.94_30.49		235.62_9.00	197.97_36.07	220.06_30.83
		RBF+PU		1		116.65_20.27	82.28_28.21	71.9616.16		120.36_16.08	74.58_10.96	77.1516.61
				2		117.15_15.00	77.14_13.39	70.52_13.79		120.75_15.61	78.10_12.98	78.60_15.63
Granada		PU		1		198.09_5.62	112.25_9.39	118.49_13.14		213.20_14.01	121.39_11.73	127.45_12.34
				2		157.57_40.79	114.70_19.86	113.14_12.95		210.26_17.62	113.15_12.98	128.26_17.07
		RBF		1		203.49_5.62	182.22_45.85	169.50_42.52		318.22_14.01	227.12_43.82	229.42_46.94
				2		186.01_40.79	168.01_31.91	181.62_41.68		316.98_17.62	246.96_46.30	270.13_51.13
		RBF+PU		1		148.56_45.33	130.61_104.41	108.39_69.97		173.50_36.95	105.6921.87	109.03_26.16
				2		152.28_59.65	112.0852.14	115.38_103.97		165.85_35.86	104.71_22.91	123.64_64.24
Huelva Cost		PU		1		147.36_18.14	61.30_10.94	68.92_10.70		139.29_9.39	62.534.05	67.97_8.01
				2		147.93_18.95	69.72_9.51	76.85_14.70		133.15_7.50	60.94_4.69	70.00_11.91
		RBF		1		135.65_18.14	87.27_18.56	109.96_46.24		158.09_9.39	114.14_38.05	141.25_41.63
				2		133.62_18.95	83.45_19.20	121.65_45.67		154.99_7.50	108.66_35.66	142.12_57.78
		RBF+PU		1		89.77_33.70	63.5417.01	71.23_29.77		77.04_15.55	69.55_20.88	79.07_20.38
				2		75.00_23.14	67.65_21.77	74.14_23.85		75.46_17.85	68.14_19.17	69.24_16.99
Jaén		PU		1		100.27_12.54	78.89_6.54	81.39_6.44		117.33_4.59	98.80_12.75	105.99_8.80
				2		86.57_16.14	78.30_6.98	77.30_8.38		113.97_12.73	96.97_10.78	103.66_11.09
		RBF		1		109.06_12.54	96.61_12.49	88.38_13.97		112.81_4.59	99.44_13.78	112.94_32.96
				2		108.08_16.14	106.57_12.88	96.79_19.70		114.19_12.73	115.94_25.38	119.28_29.18
		RBF+PU		1		74.0335.94	70.64_61.72	94.16_49.77		70.8318.46	84.69_43.50	78.41_21.37
				2		80.66_47.86	74.21_34.16	74.11_20.98		70.64_23.42	85.71_54.24	77.53_19.70
Málaga		PU		1		428.79_14.34	312.39_20.39	305.80_23.28		391.10_28.60	296.15_27.29	293.25_33.05
				2		423.88_43.40	300.56_29.19	293.06_39.57		370.87_41.41	281.45_40.26	276.09_43.11
		RBF		1		365.13_14.34	467.07_85.02	499.11_84.02		451.47_28.60	559.18_150.48	559.78_105.88
				2		395.36_43.40	499.41_68.21	528.72_68.64		448.60_41.41	564.19_107.76	531.75_107.13
		RBF+PU		1		358.95_33.49	316.25_38.65	332.67_148.19		288.31_60.17	255.6856.73	259.13_46.60
				2		352.63_26.08	312.10_48.37	298.1572.09		274.45_65.26	243.80_65.70	267.21_127.55
Sevilla		PU		1		268.46_18.32	155.93_19.72	156.78_27.77		292.42_14.57	147.5919.40	162.60_16.60
				2		226.27_22.92	174.82_23.72	187.60_37.87		270.86_16.85	141.32_22.04	158.98_19.26
		RBF		1		388.28_18.32	320.38_86.43	348.93_123.80		473.52_14.57	446.58_102.65	459.12_96.61
				2		365.89_22.92	331.61_91.40	371.61_82.34		483.75_16.85	442.58_112.32	518.32_71.35
		RBF+PU		1		263.98_47.13	171.73_39.08	146.11_34.11		273.38_58.92	150.53_23.62	152.44_21.21
				2		256.68_35.30	163.81_34.53	150.7120.47		274.94_52.97	159.38_30.53	156.38_31.06

Table 6

Performances of MoEANN and MuEANN models trained with the different combinations of input features and autoregressive orders, , evaluating the errors of the forecast horizon. The results are expressed as RMSE of the generalization set, SD stands out for Standard Deviation.

District		BF		p		MoEANN				MuEANN
						AR(p)_y	VAR(p)_y,α_4	VAR(p)_y,α_4,α_3		AR(p)_y	VAR(p)_y,α_4	VAR(p)_y,α_4,α_3
Almería		PU		1		426.60_20.90	137.2415.40	124.94_12.66		322.39_15.91	140.4022.01	131.47_15.30
				2		415.68_27.30	147.07_21.20	138.45_18.57		326.93_33.21	169.14_35.54	146.16_25.96
		RBF		1		421.46_20.90	254.38_57.33	261.43_76.61		463.63_15.91	272.98_60.76	310.54_84.80
				2		426.25_27.30	248.08_70.96	247.53_71.40		466.01_33.21	288.35_75.45	307.70_77.09
		RBF+PU		1		182.97_59.55	147.01_31.56	152.29_43.97		215.80_37.53	156.12_32.20	154.41_27.30
				2		180.70_55.17	148.54_34.73	137.52_42.93		228.23_36.45	174.99_42.88	176.23_39.83
Cádiz Bay		PU		1		501.23_25.16	325.03_11.60	322.42_24.88		409.40_18.17	276.86_29.40	269.10_22.70
				2		498.95_19.77	312.92_13.03	312.07_26.27		429.54_34.05	278.18_25.54	257.23_43.86
		RBF		1		460.89_25.16	219.5240.75	242.40_67.64		396.46_18.17	324.42_93.69	353.20_86.84
				2		454.45_19.77	217.76_48.41	258.11_57.20		410.00_34.05	381.32_401.25	326.98_77.12
		RBF+PU		1		502.89_147.86	285.99_47.07	275.16_50.38		293.54_37.96	230.03_55.89	216.9153.80
				2		442.87_107.47	286.29_52.97	267.10_67.05		311.51_112.80	220.24_58.14	212.68_60.66
Córdoba		PU		1		264.52_6.69	132.21_5.69	110.73_18.68		252.26_4.47	135.56_11.81	124.95_15.13
				2		251.62_11.44	128.85_8.86	115.51_18.41		255.80_6.53	133.14_10.32	128.61_12.48
		RBF		1		263.82_6.69	207.74_37.88	216.63_46.81		328.34_4.47	249.11_38.39	260.34_43.14
				2		258.87_11.44	223.00_26.99	228.65_28.77		338.04_6.53	258.78_41.17	265.93_30.29
		RBF+PU		1		174.55_27.13	111.63_17.62	111.06_18.58		184.47_15.81	125.0313.78	128.51_20.28
				2		172.91_30.25	121.88_19.65	110.9516.27		190.21_18.62	135.81_20.29	131.58_20.06
Granada		PU		1		279.01_21.82	152.27_9.72	164.82_18.47		301.14_12.14	139.00_7.85	158.66_12.52
				2		192.38_45.54	152.2414.59	154.67_23.55		311.78_21.74	143.1712.28	165.91_15.85
		RBF		1		294.82_21.82	220.14_37.12	212.94_52.10		412.32_12.14	284.42_41.23	300.08_64.77
				2		281.19_45.54	242.33_38.38	237.80_46.39		407.60_21.74	314.26_34.37	324.77_54.32
		RBF+PU		1		216.08_44.87	176.55_118.66	126.36_31.26		240.34_61.51	168.62_24.40	170.27_36.28
				2		207.33_44.49	196.73_123.74	158.56_75.74		241.67_53.77	169.84_27.97	196.56_69.05
Huelva Cost		PU		1		202.60_19.02	86.37_6.78	93.89_10.47		198.09_10.09	90.979.85	92.13_9.49
				2		203.57_15.97	79.06_7.22	89.88_8.48		196.18_9.77	89.03_9.64	96.41_10.44
		RBF		1		204.09_19.02	114.35_23.67	135.86_40.22		226.78_10.09	140.19_33.66	164.75_36.37
				2		189.28_15.97	105.75_22.14	137.26_39.83		227.02_9.77	134.19_30.46	169.98_67.43
		RBF+PU		1		127.05_66.98	84.60_25.29	95.66_39.06		116.44_22.26	97.03_21.84	101.24_21.02
				2		108.34_39.04	79.9220.69	96.17_32.21		114.94_25.31	96.38_22.83	94.87_28.81
Jaén		PU		1		155.56_8.00	117.36_10.04	122.96_7.68		149.90_5.03	129.34_14.86	117.03_10.98
				2		128.23_22.05	111.72_8.35	115.08_12.26		148.61_15.77	127.90_13.89	123.06_15.93
		RBF		1		140.76_8.00	121.14_13.64	121.52_14.38		181.86_5.03	160.27_20.03	166.23_27.05
				2		138.02_22.05	128.71_16.94	124.17_24.20		188.39_15.77	169.52_24.90	173.19_33.89
		RBF+PU		1		106.77_54.71	102.8334.25	109.76_36.07		108.3941.12	124.17_91.61	113.36_36.17
				2		98.82_31.06	115.55_103.98	116.29_33.37		102.82_30.77	115.75_45.87	108.93_24.21
Málaga		PU		1		511.47_33.48	410.57_22.16	395.97_34.33		308.50_26.07	326.36_27.25	325.70_30.92
				2		528.86_76.06	385.95_36.31	365.23_46.20		351.93_96.41	306.17_37.17	308.01_44.44
		RBF		1		542.03_33.48	608.60_67.13	639.00_70.48		528.34_26.07	572.49_69.65	581.85_87.50
				2		570.64_76.06	628.99_46.10	687.12_61.91		552.15_96.41	594.26_80.60	586.21_92.53
		RBF+PU		1		492.83_55.81	401.03_37.86	407.21_57.16		342.61_82.43	288.6358.38	289.45_45.87
				2		472.38_41.33	419.93_73.59	385.0262.96		352.97_78.50	268.53_62.05	303.64_111.44
Sevilla		PU		1		378.64_22.45	229.26_25.04	209.9838.14		374.89_11.51	206.52_23.67	218.54_18.82
				2		318.66_27.24	260.88_34.39	233.25_28.07		385.27_28.82	216.9826.45	226.48_25.73
		RBF		1		553.19_22.45	436.69_104.15	427.74_91.88		654.71_11.51	556.51_109.41	563.96_121.52
				2		365.89_27.24	470.94_89.03	494.17_73.37		665.62_28.82	550.15_126.87	613.40_78.47
		RBF+PU		1		394.14_61.60	227.19_66.71	196.44_56.24		449.79_92.50	252.38_41.93	251.01_30.25
				2		394.38_79.74	251.60_74.07	226.51_44.38		435.62_88.04	271.02_39.57	258.87_42.70

Table 7

Performances of MoEANN and MuEANN models trained with the different combinations of input features and autoregressive orders, , simultaneously evaluating the errors of the and forecast horizons. The results are expressed as RMSE of the generalization set, SD stands out for Standard Deviation.

District		BF		p		MoEANN				MuEANN
						AR(p)_y	VAR(p)_y,α_4	VAR(p)_y,α_4,α_3		AR(p)_y	VAR(p)_y,α_4	VAR(p)_y,α_4,α_3
Almería		PU		1		369.61_57.15	116.71_21.15	104.92_11.46		299.18_29.05	130.0425.13	120.13_11.46
				2		357.96_57.89	125.23_22.26	114.5824.19		301.37_38.27	152.93_36.74	136.45_27.93
		RBF		1		366.09_57.15	214.94_39.49	221.15_40.57		397.86_29.05	265.01_64.96	307.28_82.87
				2		359.89_57.89	216.15_32.84	225.44_23.55		409.51_38.27	271.87_81.87	304.38_74.72
		RBF+PU		1		157.43_26.44	127.65_20.66	125.89_27.47		181.67_46.53	149.47_34.12	150.46_36.18
				2		158.42_23.48	123.13_25.99	118.56_20.09		190.49_48.97	169.55_44.42	174.01_43.06
Cádiz Bay		PU		1		441.32_72.28	274.24_50.89	271.64_50.29		403.52_20.26	257.35_36.81	258.09_26.15
				2		422.43_76.62	273.75_39.36	273.80_38.58		391.77_50.73	253.12_35.98	243.53_44.69
		RBF		1		375.60_72.28	193.7926.46	207.71_35.41		346.01_20.26	296.70_96.56	317.17_86.08
				2		378.18_76.62	185.42_11.46	226.38_32.28		356.98_50.73	357.59_392.49	299.22_82.43
		RBF+PU		1		412.70_98.75	245.51_42.65	237.12_40.77		259.86_50.81	207.77_61.08	195.4760.45
				2		369.35_74.09	248.22_38.67	229.87_38.42		264.63_98.98	202.24_65.47	193.90_11.46
Córdoba		PU		1		230.13_34.48	108.05_24.94	93.58_17.78		223.15_29.95	115.48_22.78	109.32_19.84
				2		217.16_34.64	104.64_24.38	95.22_20.92		221.27_35.41	111.27_24.09	107.65_24.07
		RBF		1		226.80_34.48	182.39_26.71	192.40_25.66		277.50_29.95	222.66_44.62	235.13_48.97
				2		222.64_34.64	192.41_31.43	202.79_26.42		286.83_35.41	228.37_49.22	243.00_38.21
		RBF+PU		1		145.60_29.36	96.96_15.44	91.5119.99		152.42_35.80	99.80_11.46	102.8331.67
				2		145.03_28.29	99.79_22.74	90.74_11.46		155.48_38.75	106.96_33.51	105.09_32.01
Granada		PU		1		238.55_40.63	132.5220.25	141.66_23.50		257.17_45.88	130.1913.31	143.06_19.95
				2		174.97_18.60	133.47_19.22	133.91_21.20		261.02_54.48	128.16_11.46	147.09_25.01
		RBF		1		249.73_40.63	201.18_20.02	191.22_22.79		365.27_45.88	255.77_51.29	264.75_66.69
				2		233.60_18.60	205.17_37.63	209.71_28.87		362.29_54.48	280.61_52.86	297.45_59.41
		RBF+PU		1		182.75_34.42	153.58_25.28	117.49_11.46		206.92_60.76	137.16_39.07	139.65_44.02
				2		179.81_28.46	155.49_43.36	136.97_23.58		203.76_59.38	137.27_41.40	160.10_76.00
Huelva Cost		PU		1		174.98_27.96	73.8312.88	81.40_12.90		168.69_30.97	76.7516.09	80.05_14.94
				2		175.75_28.13	74.39_5.49	83.36_7.35		164.67_32.70	74.99_11.46	83.21_17.32
		RBF		1		169.87_27.96	100.81_14.29	122.91_14.53		192.44_30.97	127.17_38.21	153.00_40.81
				2		161.45_28.13	94.74_12.04	129.46_10.18		191.00_32.70	121.43_35.53	156.05_64.32
		RBF+PU		1		108.41_19.94	74.07_11.49	83.44_13.55		96.74_27.51	83.29_25.40	90.15_23.49
				2		91.88_17.58	73.79_11.46	85.15_12.22		95.20_29.48	82.26_25.37	82.05_26.90
Jaén		PU		1		127.92_27.83	98.13_19.45	102.44_20.96		133.62_16.98	114.07_20.61	111.51_11.38
				2		107.40_21.28	95.01_16.94	96.43_19.16		131.29_22.48	112.43_19.84	113.36_16.81
		RBF		1		124.91_27.83	109.04_12.78	104.95_16.99		147.34_16.98	129.86_34.94	139.58_40.24
				2		123.05_21.28	117.36_11.73	110.48_14.47		151.29_22.48	142.73_36.74	146.24_41.56
		RBF+PU		1		90.61_17.70	86.73_11.46	101.96_10.18		89.6136.99	104.43_74.37	95.89_34.46
				2		89.7411.04	94.88_22.28	95.20_21.72		86.73_11.46	100.73_52.43	93.23_27.09
Málaga		PU		1		472.89_41.72	361.48_49.30	351.46_45.40		349.80_49.54	311.25_31.17	309.48_35.88
				2		476.37_53.05	343.80_43.08	329.15_11.46		361.40_74.79	293.81_40.67	292.05_46.60
		RBF		1		453.58_41.72	537.83_71.30	569.05_70.50		489.90_49.54	565.83_117.44	570.81_97.75
				2		483.00_53.05	564.20_65.23	607.92_79.61		500.37_74.79	579.22_96.33	558.98_103.73
		RBF+PU		1		431.63_67.52	359.18_42.84	370.41_38.63		315.46_77.10	272.1659.87	274.29_48.66
				2		412.51_60.16	366.02_54.48	342.1444.21		313.71_82.17	256.17_11.46	285.42_121.15
Sevilla		PU		1		323.55_55.28	193.53_36.98	183.0427.22		333.65_43.27	177.05_11.46	190.57_33.13
				2		273.05_46.47	217.85_43.37	210.43_23.54		328.06_61.88	179.1544.99	192.73_40.69
		RBF		1		470.73_55.28	378.54_58.97	387.83_40.75		564.12_43.27	501.54_119.48	511.54_121.64
				2		365.89_46.47	401.27_70.31	432.89_61.91		574.68_61.88	496.37_131.34	565.86_88.79
		RBF+PU		1		332.48_65.59	199.81_28.67	171.59_11.46		361.59_117.44	201.45_61.25	201.72_55.78
				2		326.40_69.27	208.26_44.51	188.61_38.33		355.28_108.32	215.20_66.06	207.62_63.40

Apart from the model performances, the complexities of the two MoEANN and the MuEANN models are an important aspect to be compared. Since almost all experiments resulted in the maximum limited number of nodes in the hidden layer ( neurons), complexities are measured in terms of the total number of connections. Table 8 shows the average number of connections for each model, forecasting approach and dataset. For the sake of a fair comparison, model complexities are measured from a simultaneous forecast horizon perspective. For MoEANN models, the average of connections is computed for each input feature combination and autoregressive order. Then, the resulting averages for the two horizons are summed. Bold and italics represent the least complex and the second least complex architectures by forecasting approach, basis function and dataset, respectively. As expected, models become more complex when the number of autoregressive input features increases. Considering the average number of connections of the models, it can be said that the use of MuEANN models is justified. Even though the results obtained with this approach are slightly worse, the models are much simpler in all cases while their performances remain competitive.

Table 8

Comparison of MoEANN and MuEANN model complexities considering double and simultaneous forecasting, . The results are expressed in terms of the average number of links involved in the EANNs of the eight considered districts.

	BF	p	Input dataset
			AR(p)_y	VAR(p)_y,α_4	VAR(p)_y,α_4,α_3
MoEANN	PU	1	18.70	23.26	24.80
		2	20.76	26.74	29.42
	RBF	1	19.46	20.45	21.84
		2	21.14	25.50	27.74
	RBF+PU	1	20.81	24.49	26.80
		2	22.13	29.02	33.09
MuEANN	PU	1	13.04	16.20	17.03
		2	13.73	18.59	19.53
	RBF	1	14.33	14.91	15.48
		2	15.26	16.92	18.16
	RBF+PU	1	15.41	16.85	17.67
		2	15.82	19.17	20.59

In most cases, the results produced by the models involving one single output neuron (MoEANN) are slightly better than those obtained using the MuEANN approach, as seen in the results shown in Table 5, Table 6. This is to be expected since MoEANN model training is guided by considering only the error made for a single output variable, focusing the model training on a single objective. Nevertheless, it is important to point out the special case of the district of Málaga, where the results of the single and simultaneous forecast improve if the models are trained considering the errors of the two horizons in the training set. This phenomenon is due to the capability of the MuEANN model to infer relationships between the related tasks used to optimize itself.

Considering the results achieved, it can be concluded that there is no single combination of BF and dataset that performs accurately for all districts. However, some highlights can be extracted from the experimental results. Pure RBF models always produce poor performances when modeling from the MuEANN perspective, and except for Cádiz Bay, they also underperform using the MoEANN approach. In the particular case of simultaneous forecasting (see MuEANN columns of Table 7), PU architecture is the best choice when modeling the cumulative curve of contagion of the districts of Almería, Granada, Huelva and Sevilla, while an RBFPU architecture is needed for the rest of the locations. In this respect, a different number of input variable lags are needed to achieve the best results depending on the district to be modeled. Overall, considering Table 5, Table 6, Table 7, it can be seen that hybridization outperforms pure model architectures in the districts of Málaga, Jaén and Cádiz when modeling from a MuEANN approach. On the other hand, MoEANN approaches achieve better results in Sevilla, Jaén and Granada when modeling the considered horizons, either individually or jointly.

Moreover, concerning the complexities presented in Table 8, the average number of links for the eight considered districts is higher for the hybrid architectures, regardless of the autoregressive order and the number of model outputs. Pure PU models trained with are less complex than RBF models trained with the same dataset. Nevertheless, it is noteworthy that the opposite effect occurs when pure models are trained including polynomial information, i.e, when they are trained with the and datasets.

The results show that the performances obtained for both approaches (the MoEANN and MuEANN models) improve when some sort of autoregressive polynomial information derived from curve fitting is included in the model training set, except for the district of Jaén. This fact demonstrates that adding lags of polynomial coefficients has a positive effect on the forecasting task, not only in terms of accuracy but also reducing, in some cases, the variance of the results. The results obtained in Jaén may be because it is the least populated district, and the fitness of the polynomial models along its outbreaks has the higher standard deviation of the eight districts, as shown in Table 2. Fig. 6 shows a boxplot diagram where the RMSE achieved with the MuEANN model and the best dataset is compared with the results of the analogous model trained with no polynomial data considering the dual and simultaneous horizon forecast. The improvement in the results when autoregressive polynomial coefficients are included in the training data is graphically noticeable in the districts except Jaén, especially in Córdoba, Huelva Coast, Almería, Sevilla and Granada, where the median obtained with the best combination is distinctly lower than its analogous model. There is no distinguished difference in the dispersion of the results obtained for the Cádiz and Málaga districts. However, in these districts, the worst result obtained with the model is always better than the worst result obtained with the model. Additionally, to quantitatively summarize the aforementioned improvement, Table 9 presents the percentage accuracy gain of the double forecast results obtained with models trained with best with respect to training with best datasets. In terms of percentage, the results for all the districts except Jaén are notably improved when modeled with polynomial information, with Almería being the district that benefits the most from this additional information, increasing the accuracy by 73.55% and 59.46% for the MoEANN and MuEANN approaches, respectively.

Fig. 6

Comparison boxplots of simultaneous forecast results using MuEANN models including () and not including () polynomial information in the training set. The results are expressed as the RMSE obtained for the test set of the executions for each district.

Table 9

Percentage accuracy gain obtained in simultaneous double-horizon forecasting, and , training the models with the best datasets with respect to the best in all districts.

District	y_d+3,d+5
	MoEANN	MuEANN
Córdoba	37.43%	34.52%
Huelva Cost	19.68%	54.46%
Almería	73.55%	59.84%
Cádiz Bay	50.97%	26.72%
Sevilla	48.39%	46.93%
Granada	35.71%	50.90%
Jaén	4.28%	−7.49%
Málaga	30.90%	18.34%
Mean	37.61%	35.53%

To graphically illustrate the improvement produced by the inclusion of the polynomial curve point descriptors in the training phase of the EANNs, in Fig. 7, Fig. 8 the best performances obtained for the district of Sevilla with MuEANN models trained with and datasets are compared. More specifically, in Fig. 7, Fig. 8, the model test predictions are represented over time with the real curve of cumulative contagions of the four waves. For both forecast horizons, the model trained with the best dataset fits better than the one trained with the best . This improvement is especially noteworthy in the second and fourth waves for both horizons and in the third wave for the mid-term forecast. In Fig. 7, Fig. 8, where the real test values for mid- and long-term horizons are scattered with the values predicted by both models, the improvements are visibly easier to appreciate. Points located in the gray equality line represent a perfect prediction. Predictions resulting from models trained with models of the two considered forecast horizons are closer to this line, representing a more accurate prognosis of the cumulative number of contagions. Note that in Fig. 7, Fig. 8 the plotted data correspond, from left to right, to the second, fourth, fifth and third waves.

Fig. 7

test predictions of the cumulative number of Sevilla over time (a) and scattered with real values (b) using the best and datasets.

Fig. 8

test predictions of the cumulative number of Sevilla over time (a) and scattered with real values (b) using the best MuEANN model trained with and datasets.

Statistical analysis

Once all the results for the models and datasets combinations for each district have been descriptively analyzed, a statistical analysis has been carried out to obtain robust conclusions from several aspects of the experimentation. For simplicity purposes, this analysis is mainly focused on the results obtained in simultaneous horizon forecasting with the MoEANN and MuEANN approaches (see Table 7). As mentioned in Section 4.3, these models are less complex than using both MoEANNs (with one forecasting horizon each), while performances remain competitive. These tests assume normality in the data being compared. In this respect, a previous Kolmogorov–Smirnov test (Massey, 1951) was applied using the results from the different seeds before performing the comparison tests. In all cases, the -values obtained from the Kolmogorov–Smirnov test show that the null hypothesis of normality is accepted for the different sets of results used.

For reasons of clarity, several comparisons of mean tests have been performed to conclude significant differences within several aspects of the results obtained. First, the results of the best MuEANN model and training data combination for all districts are contrasted with the best results obtained from the MoEANN model, which performs better in most cases. Table 10 shows the -values obtained in a paired t-test for the mean results of the best models in both approaches. Considering a level of significance of , significant differences can only be found in the results of the districts of Córdoba, Almería and Málaga (i.e.,  out of sanitary districts). In this latter case, the MuEANN method is more accurate than the monotask method. For this specific sanitary district, in addition to being less complex in terms of total links, models that simultaneously predict mid- and long-term time forecasting horizons perform significantly better. Notwithstanding the additional difficulty of modeling two different output time horizons, this approach takes advantage of the inherent relationship between the two tasks addressed achieving similar results to the models with a single time horizon target. For the rest of the districts, the mean results obtained in both approaches are not significantly different, although the MuEANN approach is simpler in terms of the number of connections, and, thus, easier to implement.

Table 10

Statistical differences between the average mean RMSE test results of the MoEANN and MuEANN best models.

District	t	p-value
Córdoba	−3.26	2.31E−03a
Huelva Cost	−0.42	6.78E−01
Almería	−6.21	2.63E−07a
Cádiz Bay	−0.85	3.98E−01
Sevilla	−0.92	3.61E−01
Granada	−1.69	9.82E−02
Jaén	−1.18	2.46E−01
Málaga	6.75	4.69E−08b

Statistically significant differences favoring MoEANN method.

Statistically significant differences favoring MuEANN method.

To analyze the effect of the polynomial information on the performance of the models, a paired t-test was applied to the results obtained with the best and datasets for every horizon and model approach. The results, expressed as the -value obtained in the test, are presented in Table 11 for all districts, forecast horizon approaches and methodologies (MoEANN or MuEANN). With the exception of Jaén, where no significant differences can be found between the best results obtained using the or datasets, the rest of the models benefit significantly from the additional data derived from the polynomial curve fitting.

Table 11

Statistical differences between the best results obtained with the different models trained with the and datasets in the eight districts. The results are expressed as the resulting -value of a paired t-test.

District	y_d+3,i		y_d+5,i		y_d+3,d+5,i
	MoEANN	MuEANN	MoEANN	MuEANN	MoEANN	MuEANN
Córdoba	1.63E−15a	2.11E−18a	2.27E−12a	2.58E−21a	4.17E−15a	3.63E−21a
Huelva Cost	1.79E−03a	3.84E−06a	5.95E−05a	5.86E−07a	1.45E−05a	4.36E−07a
Almería	5.84E−10a	1.52E−09a	2.57E−07a	4.75E−16a	1.64E−09a	5.89E−14a
Cádiz Bay	1.91E−21a	6.66E−04a	3.59E−15a	2.54E−08a	1.56E−14a	3.67E−07a
Sevilla	6.57E−15a	3.64E−28a	5.26E−16a	2.01E−33a	6.29E−24a	1.39E−32a
Granada	4.90E−03a	8.82E−12a	4.35E−08a	7.73E−15a	1.84E−07a	1.26E−13a
Jaén	7.51E−01	1.56E−01	5.89E−01	3.22E−01	9.45E−01	2.23E−01
Málaga	2.00E−09a	3.93E−02a	8.02E−12a	3.06E−04a	3.76E−14a	1.98E−04a

Statistically significant differences favoring models trained with datasets.

Conclusions

In this paper, a novel approach for extracting inherent information of the cumulative curve of contagions produced by the COVID-19 transmission rate has been proposed. By iteratively fitting a -degree polynomial, a coefficient representation of the curve is generated for every day of a specific outbreak. The experimentation carried out, using the information of eight specific locations in Andalusia, Spain, demonstrates that using the coefficients that define the polynomial curve inflexion instant improves model performances when forecasting the mid- and long-term values of the original curve. In the experimentation carried out in this paper, evolutionary artificial neural networks (EANNs) are used to model the cumulative contagion curve of different districts using two different basis functions and considering the hybridization of both. For each location, different datasets are used to estimate the transfer function parameters using combinations of autoregressive information of the original curve and extracted polynomial features. Due to the unique and specific characteristics of the sanitary districts considered, there is no single combination of network architecture and data that perform well for all cases. However, on the whole, the methodology applied achieves excellent results.

However, from the results obtained, it can be concluded that including lags of the polynomial features to the forecaster significantly improves the results, except for the district of Jaén. In most cases, the best performances obtained in Jaén are achieved using only autoregressive data of the number of cumulative contagion time series. In this district, the methods that slightly benefit from the additional information of the polynomial coefficients (MoEANNs approaches modeling and ) return results that are not significantly different from the same approaches using only autoregressive information, as seen in the -values presented in Table 11.

Moreover, two different EANN approaches have been applied. On the one hand, monotask EANNs (MoEANN models) have been considered for modeling one target value at a time, i.e., two different models are needed to forecast on two different forecasting horizons. On the other hand, a multitask EANN (MuEANN model) has been used for forecasting both horizons simultaneously, resulting in a much simpler model. It is expected that models performing only one task are more accurate, as their optimization procedure is less complex. Nevertheless, there are no significant differences between the performances achieved by the MoEANN approaches and the MuEANN one for the sanitary districts of Huelva, Cádiz, Sevilla, Granada and Jaén, although the complexity of MuEANN models is lower in terms of average number links. Furthermore, in the district of Málaga, the MuEANN models perform significantly better than the average results of the two best MoEANN models. In this case, the MuEANN approach benefits from the intrinsic relationship between the two output targets.

Finally, it has been demonstrated that considering the features extracted from the inherent behavior of a process in the modeling stage can significantly improve the performance of the tasks carried out. Applying a similar methodology for any kind of problem, if the intrinsic problem nature is known, adapting a low-variance model as a description step may provide high quality descriptors that are able to support several supervised and unsupervised tasks.

CRediT authorship contribution statement

Miguel Díaz-Lozano: Formal analysis, Visualization, Methodology, Writing – original draft. David Guijo-Rubio: Validation, Writing – original draft, Writing – review & editing. Pedro Antonio Gutiérrez: Software, Methodology, Writing – review & editing. Antonio Manuel Gómez-Orellana: Software. Isaac Túñez: Project administration, Resources. Luis Ortigosa-Moreno: Project administration, Resources. Armando Romanos-Rodríguez: Project administration, Resources. Javier Padillo-Ruiz: Project administration, Resources. César Hervás-Martínez: Funding acquisition, Conceptualization, Methodology.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.