An Artificial Neural Network-Based Approach for Predicting the COVID-19 Daily Effective Reproduction Number R_t in Italy

Abstract

Since December 2019, the novel coronavirus disease (COVID-19) has had a considerable impact on the health and socio-economic fabric of Italy. The effective reproduction number R_t is one of the most representative indicators of the contagion status as it reports the number of new infections caused by an infected subject in a partially immunized population. The task of predicting R_t values forward in time is challenging and, historically, it has been addressed by exploiting compartmental models or statistical frameworks. The present study proposes an Artificial Neural Networks-based approach to predict the R_t temporal trend at a daily resolution. For each Italian region and autonomous province, 21 daily COVID-19 indicators were exploited for the 7-day ahead prediction of the R_t trend by means of different neural network architectures, i.e., Feed Forward, Mono-Dimensional Convolutional, and Long Short-Term Memory. Focusing on Lombardy, which is one of the most affected regions, the predictions proved to be very accurate, with a minimum Root Mean Squared Error (RMSE) ranging from 0.035 at day t + 1 to 0.106 at day t + 7. Overall, the results show that it is possible to obtain accurate forecasts in Italy at a daily temporal resolution instead of the weekly resolution characterizing the official R_t data.

1. Introduction

Since December 2019, the entire world has been fighting the COrona VIrus Disease-2019 (COVID-19) pandemic, whose first case was identified in Wuhan, China [1]. The virus responsible for this new infection is named “SARS-CoV-2” [2], and the main symptoms experienced by the infected population are fever, fatigue, and dry cough that, in severe cases, can lead to pneumonia, severe acute respiratory syndrome, and death [3]. World governments have applied several containment measures to limit the virus spread, which mainly occurs through close contact. In Italy, after a strict lockdown imposed from March to May 2020 [4], different containment measures have been applied, starting from 3 November 2020 [5], by identifying three risk scenarios related to increasing critical levels of the pandemic (yellow, orange, and red) in the various Italian regions. Specifically, the risk scenario for each region and the corresponding color are attributed by analyzing a set of epidemiological parameters, such as the effective reproduction number R_t and the number of regional hospitalizations. Monitoring the R_t parameter is both strategic and important, not only for establishing the specific risk level associated with each region, but also to assess the magnitude of the effort required to control the pandemic spread. Accurate estimates of the R_t value are crucial for decision makers to control the spread of the disease and to plan containment measures accordingly [6]. As a consequence, the temporal resolution of the parameters to be estimated has a fundamental role in develo** predictive models in order to facilitate the timely application of appropriate containment measures. For this reason, the present study focused on all 21 Italian regions, aiming to provide an approach based on Artificial Neural Networks (ANNs) for predicting the effective reproduction number (R_t). Specifically, six ANNs were designed, trained, and then compared to predict the R_t indicator 7 days ahead in each Italian region, starting from 21 epidemiological input variables at a given day. Indeed, the novelty of the present study consists in providing a methodology able not only to estimate the R_t indicator (as several approaches already do), but also to predict it forward in time at a daily resolution. Although the official R_t data are provided by the Italian National Institute of Health with a weekly temporal resolution, this indicator exhibits an inherent daily variability being directly related to several COVID-19 indicators (e.g., the number of positive cases, the number of hospitalized patients with symptoms, the intensive care units occupation, the number of people home confined, etc.), which are monitored with a daily temporal resolution. Therefore, obtaining daily predictions of the R_t represents a crucial aspect because a finer temporal resolution constitutes a benefit for tracking the disease spread and supporting the corresponding policy makers decisions.

Related Works

Several studies have exploited Artificial Intelligence (AI) and mathematical models to track and forecast the COVID-19 spread. In particular, a review of several techniques proposed by the scientific community to face the pandemic can be found in [7,8]. Specifically, Clement et al. [7] reviewed the mathematical models, such as compartmental and statistical ones, as well as machine learning models that are currently exploited for different tasks related to the COVID-19 context (e.g., disease diagnosis in infected people, forecast of several key epidemiological indicators, etc.). Therefore, the related works can be grouped into three main categories in respect of the used methodology: (i) mathematical-statistical modeling, (ii) machine learning approaches to support traditional modelling solutions, and (iii) approaches fully based on machine learning.

Regarding the use of mathematical models, Zou et al. [9] proposed an extension of the Susceptible-Exposed-Infected-Recovered (SEIR) epidemic model named SuEIR for predicting COVID-19 spread in the USA, including the computation of the basic reproduction number R₀. Moreover, Zhou et al. [10] exploited a compartmental SEIR model to compute the basic reproduction number R₀ related to COVID-19 pandemic in China. In detail, the SEIR model hyper-parameters were included in a closed formulation that allows the computation of the R₀ value. Each hyper-parameter was estimated by analyzing observed epidemiological data or, in the absence of recent scientific studies, by reporting the analogue value estimated for the SARS epidemic outbreak that occurred in 2002. Furthermore, Peirlinck et al. [11] applied a mathematical model and a Bayesian framework to infer the COVID-19 epidemiological characteristics by exploiting symptomatic case data combined with antibody seroprevalence studies. Besides computing the time-varying contact rate of the outbreak in real-time, the obtained model projects the temporal evolution as well as the credible intervals of the effective reproduction number, along with figures regarding the symptomatic, asymptomatic, and recovered population. In addition, Giuliani et al. [12] proposed an endemic–epidemic time-series mixed-effects generalized linear model for predicting spatio-temporal diffusion of the COVID-19 phenomenon in Italy.

Recently, machine learning approaches have also been widely exploited to support traditional modelling solutions. Indeed, Shorten et al. [8] demonstrated that AI models can be combined with traditional approaches such as statistical modelling, the Reverse Transcriptase-Polymerase Chain Reaction (RT-PCR) COVID-19 testing procedure, and Susceptible Infected Recovered (SIR) compartmental modelling. Regarding the task of forecasting key epidemiological indicators, a hybrid approach based on machine learning was developed in [13] for overcoming the inherent limitations of SIR/SEIR modelling. Moreover, Zheng et al. [14] also proposed a hybrid approach in which a Natural Language Processing (NLP) module and a Long Short-Term Memory (LSTM) network are embedded with an Improved-Susceptible-Infected (ISI) compartmental model to produce COVID-19 predictions. Deep learning techniques were also exploited by Deng [15] to estimate transmission parameters of a customized compartmental model, with the aim of simulating the dynamics of the COVID-19 pandemic in the United States and providing 35- and 42-day predictions. Furthermore, Abbasimer and Paki [16] combined multi-head attention models, Convolutional Neural Networks (CNNs), and LSTM neural networks for forecasting COVID-19 daily confirmed cases. A Bayesian optimization procedure was exploited to define the parameters of the machine learning models and enhance forecasting performance. Solanki and Singh [17] compared ARIMA, SARIMAX, Polynomial regression, and Long Short-Term Memory models when handling COVID-19 time-series data.

Regarding the approaches fully based on machine learning, Punn et al. [18] exploited a set of machine learning models (polynomial regression and recurrent neural networks) to forecast COVID-19 confirmed and recovered cases and deaths in several countries, dealing with John Hopkins University epidemiological data. Moreover, Zeroual et al. [19] compared different deep learning methods, such as LSTM NNs, Recurrent NNs, and Variational Auto-Encoders (VAEs), for predicting COVID-19 daily confirmed and recovered cases in Italy, Spain, France, China, the USA, and Australia. In addition, Davahli et al. [20] exploited Graph Neural Networks (GNNs) based on the graph theory to forecast the COVID-19 effective reproduction number R_t trend in the USA.

The present study provides an approach fully based on ANNs, thus not requiring a-priori knowledge for building and fitting the models, as opposed to traditional compartmental frameworks, which need a preliminary specification of disease-related hyperparameters. To the best of the authors’ knowledge, this is the first study providing a framework able to predict the effective reproduction number R_t up to seven days ahead at a daily resolution for each Italian region, contrary to other studies and official reports, which present only its estimation at a coarser temporal resolution.

2. Materials and Methods

2.1. Data Sources and Preprocessing

Data concerning the COVID-19 epidemiological indicators in the period from 24 February 2020 to 6 May 2021 were collected from the publicly available Italian Civil Protection Department (ICPD) repository [21]. In particular, the considered indicators are: (i) hospitalized patients with symptoms, (ii) intensive care units occupation, (iii) total hospitalized patients, (iv) home confinement, (v) total amount of daily positive cases, and since the beginning of the pandemic, (vi) recovered, (vii) deaths, (viii) total amount of positive cases, and (ix) number of swabs performed. These data were collected for all the Italian regions and autonomous provinces: Abruzzo, Aosta Valley, Apulia, Basilicata, Calabria, Campania, Emilia-Romagna, Friuli Venezia Giulia, Lazio, Liguria, Lombardy, Marche, Molise, Piedmont, Sardinia, Sicily, Autonomous Province (A.P.) of Bolzano, A.P. of Trento, Tuscany, Umbria, and Veneto. The ICPD provides daily updates to data on the repository by adding new records and eventually correcting past misreported records. Because the collected data are cumulative, a pre-processing step was performed to obtain incremental daily quantities from the cumulative ones (x–xviii). In addition to the nine cumulative indicators and their incremental daily counterpart, three further variables were considered, relating to the daily reproduction number estimates, including the upper and lower bounds of 99% Confidence Interval (CI) (xix–xxi), for a total of 21 variables. Furthermore, data concerning the R_t estimates are weekly updated by the Italian National Institute of Health and made available to the public through official reports [22]. These data, along with their CI, were gathered and then used in the present work for weekly reference against forecasted values of the daily effective reproduction number.

2.1.1. Daily Effective Reproduction Number Estimation

The effective reproduction number, R_t, is one of the most important epidemiological indicators. Its evaluation allows assessing both the actual epidemic risk and the effort required to control the spread of the disease. The R_t indicator is defined as the number of secondary infections that arise from a typical primary case in a partially immunized population. As a consequence, the magnitude of the effective reproduction number evaluates the risk related to the spread of an infectious disease. Specifically, when R_t > 1, a primary infection will cause more than one secondary infection, meaning that the disease is in a spreading phase; on the other hand, when R_t < 1, the infection is in a containment stage, because a primary infection causes less than one secondary infection. Hence, accurate estimates of this indicator are crucial to plan suitable containment measures to tackle the infection. Although the Italian National Institute of Health provides weekly estimates of the R_t index through a statistical approach based on the Monte-Carlo Markov Chain (MCMC) algorithm [23], effective reproduction number estimates at a finer daily temporal scale are needed for the purposes of the present study. Therefore, the daily estimates of the effective reproduction number, R_t, were obtained using the method proposed by Wallinga and Lipsitch [6] based on a rephrased formulation of the Lotka–Euler equation in the context of infectious disease epidemiology. Specifically, the Lotka–Euler equation was exploited to infer a relationship between the effective reproduction number (denoted as R in Wallinga–Lipsitch equations), the growth rate, r, and the generation interval distribution g(t):

$$\frac{1}{R}={{\displaystyle \int }}_{t=0}^{t=\infty }{e}^{-rt}g\left(t\right)dt$$

(1)

where r is defined as the per capita change in number of new cases per unit of time, g(t) is the probability density function that models the time interval between the infection of an individual and the secondary infection case caused by that individual, and t indicates the time unit (i.e., days in the present study). Contrary to g(t), which is a-priori defined with a particular probability distribution (for Italy, a gamma function with rate = 1.84 and shape = 0.28 [24,25]), the growth rate needs to be empirically estimated for each day, k, in the reference period, exploiting observed data concerning new daily positive individuals. The present study exploited an iterative procedure based on a fixed size sliding window ($w_{size}$ = 17), as reported in Algorithm 1. Notice that a greater window size allows us to obtain a more accurate estimate of the growth rate and, consequently, of the daily effective reproduction number. However, if the number of observational data is a concern (especially at early stages of the disease), a larger window size would not be appropriate because the higher the length of the window, the higher the time instants (in this work, days) for which the growth rate estimation cannot be performed.

Algorithm1 Daily growth rate estimation

Input:p (time series of new positive cases)

Output:rlist (growth rate for each day in the list)

Initializations:

1:       wsize = Integer odd number 2:       hwindow = (wsize−1)/2 3:       wdays = arange(0, wsize) 4:       rlist = [] 5:       L = length(p)

LOOP Process: 6:       for k in range(hwindow, L−hwindow, step = 1) do 7:                 left = k 8:                 right = k + hwindow 9:                 y = p[left: right] 10:               model = Poisson(y ~ 1 + wdays).fit() 11:               r = model.coeff(wdays) 12:               rCI = model.coeff(wdays).CI(99%) 13:               rlist.append([r,rCI]) 14:      end for 15:      return rlist

Being able to estimate a daily trend for the growth rate allows computing of the daily effective reproduction number in a specific day, k, as follows:

$$R\left(k\right)=\frac{1}{{\sum }_{t=t_1}^{t_2}g\left(t\right)e^{-r\left(k\right)t}} \forall k \mathrm{in} \left[t_1;t_2\right]$$

(2)

with $t_1=t_0+\frac{w_{size}}{2}; t_2=t_f-\frac{w_{size}}{2}$, where $t_0$ is the first day on which data regarding new positive cases are available, $t_f$ is the last day on which these data are provided, and w_size is the sliding window size in Algorithm 1. As an example, Figure 1 reports the Wallinga–Lipsitch (W–L) daily estimates of the COVID-19 effective reproduction number R_t in Lombardy, which turned out to be one of the most affected Italian regions. Data of new positive cases from 24 February 2020 to 19 April 2021 were exploited by Algorithm 1 for computing the daily growth rate, which was then used in (2) to estimate the daily effective reproduction number, along with its 99% CI in the period from 3 March 2020 to 11 April 2021. R_t estimates in the period from 24 February to 2 March 2020 and from 12 to 19 April 2021 are not available because the W–L algorithm requires a sliding window-based procedure that does not provide estimates for a number of days equal to half of the window length at both the beginning and the end of the considered time horizon. However, this sliding window approach allows one to obtain more accurate results with a narrow confidence interval compared to methods based on daily observations of epidemiological indicators, such as the Wallinga–Teunis approach [26].

2.2. Artificial Neural Network (ANN) Architectures

Because there was no a priori knowledge concerning the mathematical relationship between the 21 COVID-19 indicators and the target (i.e., the R_t temporal trend), different ANN architectures were adopted to exploit their characteristics and strengths. Specifically, three complementary architectures of ANN (each with two configurations of hidden layers) were selected, designed, and implemented to forecast the daily R_t index trend in all the 21 Italian regions and autonomous provinces. Each architecture was designed for the 7-day ahead prediction of the daily reproduction number by taking as input a set of 21 epidemiological indicators at a specific day, as defined in Section 2.1. In particular, the prediction task in this work made use of the following architectures, each with 1 to 2 hidden layers, according to their characteristics:

• a Fully Connected Neural Network (FCNN), which represents a baseline NN;
• a Mono-dimensional Convolutional Neural Network (1D CNN), used for extracting the inherent information (i.e., the internal representation of features) of a one-dimensional sequence of observations, such as time series data;
• a Long Short-Term Neural Network (LSTM), typically used for selectively retaining the long- to short-term temporal relationships included in sequence data;

Their design details are reported in

Table 1. Design details of the models.

Model	1-Layer FCNN	2-Layers FCNN	1-Layer 1D-CNN	2-Layers 1D-CNN	1-Layer LSTM	2-Layers LSTM
#Units (per layer)	42	21; 21	-	-	42	21; 21
#Filters (per layer)	-	-	21	10; 10	-	-
Kernel size (per layer)	-	-	1	1; 1	-	-
Stride (per layer)	-	-	1	1; 1	-	-
Activation function (per layer)	ReLU	ReLU; ReLU	Sigmoid	ReLU; ReLU	Tanh	Tanh; Tanh
Weights Regularizer (per layer)	L1 (0.001)	L1 (0.01); L1 (0.01)	L1 (0.001)	L1 (0.001); L1 (0.001)	L1(0.001)/L2(0.01)	L1 (0.001)/L2(0.01); L1 (0.001)/L2(0.01)
Bias Regularizer (per layer)	L1 (0.001)	L1 (0.01); L1 (0.01)	L1 (0.001)	L1 (0.001); L1 (0.001)	L1(0.001)/L2(0.01)	L1(0.001)/L2(0.01); L1 (0.001)/L2(0.01)
Weights Initializer (per layer)	Uniform	Uniform; Uniform	Uniform	Uniform; Uniform	Glorot	Uniform; Uniform
Bias Initializer (per layer)	Uniform	Uniform; Uniform	Uniform	Uniform; Uniform	Glorot	Uniform; Uniform

ReLU = Rectified Linear Unit, L1 = Lasso Regularization, Tanh = Hyperbolic Tangent, Sigmoid = Sigmoid Function, Glorot = Glorot Uniform initializer, FCNN = Fully Connected Neural Network, 1D-CNN = Mono-Dimensional Convolutional Neural Network, LSTM = Long Short–Term Memory.

.

2.2.1. Fully Connected Neural Networks

A fully connected ANN consists of a series of fully connected layers so that each neuron in a layer is connected to the neurons in the subsequent layer. Neuron output is typically regulated by non-linear activation functions that allow the ANN to learn non-linear patterns in the input data. Additionally, the number of layers, along with the number of neurons per layer, which are both design parameters, allow for more complex tasks [27]. The main advantage of these architectures is that they are “structure agnostic”, which means that there are no specific assumptions to be made on the input. Hence, the input may be represented by images, videos, or mono-dimensional data. As a consequence, this characteristic makes the fully connected ANN very broadly applicable to several case studies. On the other hand, such networks tend to have weaker performance than special-purpose networks, which are designed and tuned to the specific structure of a problem. In the present study, two fully connected ANNs were adopted, which differ for the number of hidden layers and the number of neurons in each layer. Specifically, the 1-layer FCNN is composed of an input layer with 42 neurons, whereas the 2-layers FCNN is designed with 21 neurons in each hidden layer. For both models, the output layer contains 7 neurons, one for each of the subsequent 7 days of prediction.

2.2.2. 1-D Convolutional Neural Networks (1D CNNs)

Convolutional Neural Networks (CNNs) were first proposed for dealing with multi-dimensional input, such as images, because these models succeed in processing large inputs by fusing the feature extraction and feature classification processes into a single learning schema through the application of convolutional operations. Mono-dimensional CNNs have recently been introduced for dealing with mono-dimensional (1D) data. The main difference between 1D and 2D CNNs consists of 1D arrays that replace 2D matrices for both kernels and feature maps. Indeed, each convolutional operation in a 1D CNN requires lower computational complexity with respect to 2D CNNs, because it is applied to vectors rather than matrices. Kiranyaz et al. [28,29] underlined how well 1D CNNs performed when dealing with limited labeled data with high variation, and the epidemiological data exploited in the present study have the same characteristics. A 1D-CNN architecture is composed of a series of mono-dimensional convolution layers featuring a kernel size, a stride, and an activation function, and it also has an output layer made up of fully connected neurons. The present work implemented two 1D CNNs with a single multi-layer perceptron output layer made up of seven fully connected neurons, one for each time step to be predicted.

2.2.3. Long Short-Term Memory (LSTM) Neural Networks

LSTM neural networks are an advanced version of the recurrent ANN architecture that allow the solution of the vanishing gradient problem affecting the learning of long-term dependencies [30,31,32]. In particular, an LSTM ANN is composed of a sequence of LSTM cells. Each LSTM cell comprises three gates, i.e., forget, input, and output gates. Before dealing with the rule of the three gates above, it is necessary to specify that an LSTM architecture succeeds in handling both short- and long-term dependencies among input data by kee** track of a hidden state and a cell state. The hidden state is used to track a short memory dependency, whereas the cell state tracks a long-term relationship among input data. Specifically, in an LSTM cell, both the hidden state and the cell state referring to the previous time-step are given as input to the forget gate, which is responsible for deciding what information to erase and what to propagate ahead.

The input gate is then exploited to measure the relevance of the new information borrowed from the last input data and to merge the new information with the previous hidden state for updating the cell state. Finally, the output gate is responsible for elaborating a new hidden state to be assumed as the input of the subsequent cell. In general, LSTM ANNs are ideal when working with input data that exhibit short- and/or long-term dependencies; this could be the case of epidemiological indicators, whose trend can be analyzed as a time-sequence.

The output layer of an LSTM ANN is typically defined by a FCNN made up of a set of neurons corresponding to the number of predicted time-steps. The present work implemented two LSTM ANNs by taking as input data regarding 21 epidemiological indicators (reported in Section 2.1) at a single time step.

2.3. Experimental Setup

For each Italian region and autonomous province, the present work exploited data concerning 21 epidemiological indicators (input variables) in order to obtain seven-day ahead forecasts of the R_t index (target variable) every day, by first deriving estimates through the W–L method proposed in Section 2.1.1. Input data were scaled in [0, 1] to avoid numerical problems and to speed up the convergence of ANNs training procedures. Some regions required excluding a set of training data that had caused strong uncertainty in the W–L R_t estimates. In particular, Aosta Valley, Molise, Campania, Basilicata, and Calabria required the elimination of epidemiological data in the period from 3 March 2020 to 1 September 2020, whereas the R_t estimates in the period from 25 May 2020 to 10 August 2020 were excluded for Apulia, Piedmont, Liguria, Friuli Venezia Giulia, A.P. Trento, A.P. Bolzano, Umbria, and Sardinia. For the remaining regions, the entire available dataset was exploited (from 3 March 2020 to 11 April 2021). Therefore, given a set of 21 epidemiological indicators at a certain day, the forecasting task consists of predicting the temporal trend of the effective reproduction number in the subsequent seven days. To do so, the aforementioned six ANNs were trained and compared on the entire dataset by means of the K-fold Cross Validation method. This approach consists of two subsequent steps: (i) a preliminary splitting of the entire dataset into K non-overlap** subsets; (ii) an iterative procedure that, for each iteration, assumes a subset as the cross validation set and the remaining K − 1 subsets as the training set. Notice that, for recurrent architectures such as LSTM networks, K-fold Cross Validation slightly changes in order to ensure time consistency of the training and cross validation sets in each iteration. Specifically, this requires splitting the dataset into K non-overlap** ordered subsets so that, for each iteration $j\in 1, ... , K$, the first j contiguous subsets are used for training, whereas the $j+1$ subset is assumed as the cross validation set, kee** temporal dependency unaltered. To assess the performance of the proposed ANN architectures at a regional context, three error metrics were computed on the validation set: Root Mean Squared Error (RMSE), Spatial RMSE (RMSE_S), and Temporal RMSE (RSME_T), during each iteration of the K-fold Cross Validation method. The aforementioned metrics are formally defined as follows:

$$RMSE=\sqrt{\frac{1}{m}{\displaystyle \sum }_{i=1}^m{(y_i-\widehat{y_i})}^2}$$

(3)

$$RMS{E}_S=\frac{1}{t_f}{\displaystyle \sum }_{t=1}^{t_f}{\left(\sqrt{\frac{1}{m}{\displaystyle \sum }_{i=1}^m{(y_i\left(t\right)-\widehat{y_i}\left(t\right))}^2}\right)}_t$$

(4)

$$RMS{E}_T=\frac{1}{m}{\displaystyle \sum }_{i=1}^m{\left(\sqrt{\frac{1}{t_f}{\displaystyle \sum }_{t=1}^{t_f }{(y_t\left(i\right)-\widehat{y_t}\left(i\right))}^2}\right)}_i$$

(5)

where m represents the number of samples included in the k-th cross validation set and t_f is the lead-time, i.e., the number of time steps to predict ahead (t_f = 7 in this work). Because the K-fold Cross Validation procedure involves K iterations, it follows that, for each architecture, a set of K values is obtained for each error metric. Thus, an average value was computed across the aforementioned K values for each of the three-error metrics reported in (3)–(5). In this way, the Final Score (FS) for the six architectures was computed as follows:

$$FS=\frac{\left(\frac{1}{K} {\sum }_{k=1}^{K}RMSE\left(k\right)\right)+\left(\frac{1}{K} {\sum }_{k=1}^{K}RMS{E}_V\left(k\right)\right)+\left(\frac{1}{K} {\sum }_{k=1}^{K}RMS{E}_H\left(k\right)\right) }{3}$$

(6)

The three networks with lower FS values were selected for each region and autonomous province. The K-Fold Cross Validation algorithm in this work was executed with K = 10. Thus, a total of 10 (runs) × 6 (architectures) × 21 (regions) = 1260 fitting were performed only in this comparison phase. After the selection of the three best architectures, a further training procedure was performed on 70% of the entire dataset (training set), whereas the validation and test procedures required 30% of the remaining data (15% for validation and the remaining 15% for testing). Furthermore, training and validation data were preliminarily shuffled before passing to the models, except for recurrent architectures, because this would lead to breaking temporal consistency. In particular, the three models with the lowest FS were exploited after the training procedure to forecast the R_t index in each region. The forecasts produced by the models were further averaged to compute an ensemble prediction that allows for inherent training patterns learned through different architectures. Moreover, in order to assess both the uncertainty of forecasts and the random weights initialization, the 95% Prediction Interval (PI) was computed for the best model, according to Algorithm 2. The present study adopted a K value equal to 10 for evaluating the prediction intervals. This leads to a total of 10 (iterations) × 21 (regions/autonomous provinces) = 210 fittings performed during the PI computation phase.

Algorithm2 95% Prediction Interval (PI) estimation

Input:X (input data), K (number of iterations), s_train (training set size), s_val (validation set size), s_test (test set size), modeltype (architecture type)

Output: ${\hat{y}}_{MEAN}$ (average predictions), ${\hat{y}}_{LB}$ (lower bound of PI), ${\hat{y}}_{UB}$ (upper bound of PI)

Initializations:

1:       train = X[s_train] 2:       validation = X[s_val] 3:       test = X[s_test] 4:       modelspredictions = [ ]

LOOP Process: 5:       for j in range(1, K, step = 1) do 6:                model = ModelGenerator(modeltype) 7:                model.fit(train) 8:                model.validate(validation) 9:                predictions = model.predict(test) 10:              modelspredictions.append(predictions) 11:      end for 12:      ${\hat{y}}_{MEAN}$ = mean(modelspredictions) 13:      ${\hat{y}}_{STD}$ = standard_deviation(modelspredictions) 14:      ${\hat{y}}_{LB} = {\hat{y}}_{MEAN} - (1.96 \times {\hat{y}}_{STD})$ 15:      ${\hat{y}}_{UB} = {\hat{y}}_{MEAN} + (1.96 \times {\hat{y}}_{STD})$ 16:      return ${\hat{y}}_{MEAN}$, ${\hat{y}}_{LB}$, ${\hat{y}}_{UB}$

Each model was trained by means of the Adam optimizer, and all the fitting parameters are listed in

Table 2. Fitting parameters of the models.

Model	1-Layer FCNN	2-Layers FCNN	1-Layer 1D-CNN	2-Layers 1D-CNN	1-Layer LSTM	2-Layers LSTM
Batch Size	4	4	2	2	4	4
Learning Rate	0.0001	0.0001	0.0001	0.0001	0.001	0.001
Decay Rate	0.9	0.9	0.9	0.9	0.9	0.9
Loss Early Stop** (patience)	50	50	50	50	50	50
Validation Loss Early Stop** (patience)	50	50	50	50	50	50
Epochs	4000	4000	4000	4000	4000	4000
Shuffle	True	True	True	True	False	False

FCNN = Fully Connected Neural Network, 1D-CNN = Mono-Dimensional Convolutional Neural Network, LSTM = Long Short–Term Memory.

. As regards the software adopted for model implementation, both the architecture and training/test control flows were written in Python (v3.7.3), based on the Keras API (v2.3.1) [33], and relying on TensorFlow backend (v2.2.0) [34]. The training phase was performed in a distributed fashion by means of the Python multiprocessing framework in order to speed-up the computation. Furthermore, several Python libraries were exploited throughout the experiments, i.e., Numpy (v.1.16.2), Pandas (v.1.0.3), Sklearn (v.0.22.2), Scipy (v.1.4.1), Statsmodels (v.0.11.1), and Matplotlib (v.3.0.3). The execution of the entire workflow, including the K-Fold Cross Validation procedure along with the training of the best three models and PI computation, took approximately one hour for each region and autonomous province.

2.4. The Rolling Approach

As detailed in the previous sections, starting from 21 epidemiological variables at a given input day, t, the proposed framework provides seven-day ahead forecasts (f(t + 1), …, f(t + 7)) of the effective reproduction number, R_t, with an increasing prediction error as the forecast horizon increases. As a consequence, a rolling approach could be exploited to ensure the most accurate R_t forecast for each input day. Specifically, because all the required input data are provided daily by ICPD, it is only possible to collect the one-day ahead prediction, i.e., the most accurate forecast, for each input day, discarding all the others. This is made possible because new data are published daily, which allows the re-application of the same procedure at the subsequent days by recomputing the previously discarded predictions with a higher accuracy. Clearly, this procedure, called the rolling approach, is no longer applicable if new input data are not available daily. In this case, all seven-day ahead forecasts obtained from the input day, t, must be collected. The aforementioned operating workflow is described in Figure 2.

3. Results

Figure 3 reports the comparison between the test-set R_t ground-truth in Lombardy and the predictions obtained through the three regional best models, which ranked in the top 3 after the K-Fold Cross Validation procedure (#1, #2, and #3, respectively), along with their average ensemble. Specifically, the R_t ground truth was estimated by means of the W–L procedure, as shown in Figure 1. Figure 3 highlights that the prediction error increases as the predicted time grows, leading to inaccurate forecasts for the last instants of the prediction window (from five to seven days ahead). On the contrary, the forecasts obtained for the remaining prediction window time instants (from one to four days ahead) proved to be very accurate. Moreover, it should be noted that considering the three regional best models in the predictive phase is crucial because, as shown in Figure 3, the #2 model (1-layer 1D-CNN) is able to predict the test-set R_t trend better than the #1 model (1-layer FCNN). For each Italian region/autonomous province, the RMSEs evaluated on the test set are reported in

Table 3. Regional RMSE on the test set for the top three architectures.

Region	Architecture	f(t + 1) RMSE	f(t + 2) RMSE	f(t + 3) RMSE	f(t + 4) RMSE	f(t + 5) RMSE	f(t + 6) RMSE	f(t + 7) RMSE
Abruzzo	1-L FCNN	0.071	0.113	0.144	0.172	0.221	0.249	0.252
	1-L 1D CNN	0.08	0.104	0.083	0.109	0.095	0.11	0.106
	2-Ls 1D CNN	0.107	0.113	0.119	0.114	0.12	0.121	0.117
	Ensemble	0.085	0.109	0.113	0.131	0.142	0.157	0.156
Aosta Valley	1-L FCNN	0.101	0.113	0.137	0.155	0.163	0.158	0.149
	1-L 1D CNN	0.137	0.128	0.151	0.161	0.167	0.163	0.162
	2-Ls 1D CNN	0.137	0.13	0.144	0.163	0.175	0.172	0.177
	Ensemble	0.119	0.121	0.142	0.159	0.168	0.164	0.162
Apulia	1-L FCNN	0.043	0.059	0.065	0.075	0.083	0.086	0.09
	2-Ls 1D CNN	0.098	0.094	0.098	0.098	0.107	0.105	0.111
	1-L 1D CNN	0.094	0.101	0.081	0.101	0.11	0.115	0.115
	Ensemble	0.075	0.081	0.079	0.089	0.098	0.1	0.104
Basilicata	1-L LSTM	0.136	0.134	0.132	0.137	0.144	0.166	0.179
	1-L FCNN	0.066	0.095	0.115	0.143	0.172	0.206	0.217
	1-L 1D CNN	0.066	0.091	0.105	0.124	0.138	0.165	0.178
	Ensemble	0.078	0.101	0.115	0.134	0.149	0.176	0.187
Calabria	1-L LSTM	0.134	0.159	0.177	0.192	0.209	0.222	0.225
	1-L FCNN	0.064	0.096	0.119	0.139	0.163	0.187	0.187
	2-Ls 1D CNN	0.12	0.12	0.11	0.107	0.107	0.113	0.105
	Ensemble	0.105	0.124	0.134	0.144	0.158	0.171	0.169
Campania	1-L FCNN	0.032	0.052	0.065	0.072	0.081	0.088	0.094
	2-Ls 1D CNN	0.032	0.055	0.067	0.072	0.08	0.091	0.093
	1-L 1D CNN	0.032	0.054	0.065	0.073	0.083	0.093	0.098
	Ensemble	0.031	0.053	0.066	0.072	0.081	0.091	0.095
Emilia-Romagna	1-L FCNN	0.027	0.038	0.043	0.049	0.06	0.069	0.077
	2-Ls 1D CNN	0.066	0.09	0.112	0.134	0.16	0.18	0.198
	1-L 1D CNN	0.025	0.039	0.043	0.057	0.073	0.084	0.093
	Ensemble	0.031	0.042	0.053	0.065	0.081	0.095	0.107
Friuli Venezia Giulia	1-L FCNN	0.072	0.108	0.149	0.178	0.209	0.229	0.237
	1-L 1D CNN	0.074	0.108	0.136	0.18	0.207	0.214	0.238
	2-Ls 1D CNN	0.15	0.179	0.2	0.246	0.285	0.286	0.309
	Ensemble	0.096	0.13	0.161	0.2	0.233	0.243	0.261
Lazio	1-L FCNN	0.026	0.04	0.052	0.062	0.074	0.085	0.092
	2-Ls 1D CNN	0.167	0.168	0.173	0.173	0.165	0.163	0.167
	1-L 1D CNN	0.115	0.139	0.127	0.108	0.123	0.112	0.124
	Ensemble	0.098	0.105	0.103	0.098	0.1	0.096	0.103
Liguria	1-L FCNN	0.061	0.083	0.104	0.118	0.131	0.148	0.162
	1-L 1D CNN	0.035	0.048	0.065	0.079	0.092	0.102	0.11
	2-Ls 1D CNN	0.039	0.067	0.077	0.089	0.101	0.115	0.126
	Ensemble	0.04	0.061	0.078	0.092	0.106	0.119	0.13
Lombardy	1-L FCNN	0.041	0.059	0.077	0.088	0.099	0.108	0.117
	1-L 1D CNN	0.035	0.048	0.061	0.069	0.083	0.096	0.106
	2-Ls 1D CNN	0.084	0.093	0.101	0.109	0.118	0.129	0.138
	Ensemble	0.046	0.062	0.077	0.085	0.097	0.108	0.118
Marche	1-L FCNN	0.067	0.101	0.137	0.156	0.181	0.207	0.22
	2-Ls 1D CNN	0.185	0.177	0.151	0.144	0.155	0.145	0.149
	2-Ls FCNN	0.045	0.064	0.078	0.1	0.117	0.128	0.137
	Ensemble	0.091	0.107	0.118	0.131	0.148	0.156	0.164
Molise	1-L LSTM	0.097	0.112	0.128	0.133	0.142	0.144	0.15
	1-L FCNN	0.091	0.098	0.121	0.12	0.123	0.118	0.136
	1-L 1D CNN	0.092	0.119	0.136	0.146	0.155	0.167	0.18
	Ensemble	0.093	0.108	0.127	0.131	0.139	0.141	0.154
Piedmont	1-L FCNN	0.061	0.083	0.092	0.118	0.12	0.136	0.151
	1-L 1D CNN	0.037	0.051	0.065	0.071	0.079	0.086	0.092
	2-Ls 1D CNN	0.069	0.087	0.101	0.115	0.131	0.14	0.144
	Ensemble	0.045	0.063	0.074	0.088	0.099	0.109	0.116
Sardinia	1-L FCNN	0.158	0.213	0.233	0.252	0.286	0.315	0.346
	1-L 1D CNN	0.07	0.093	0.122	0.147	0.175	0.199	0.217
	2-Ls 1D CNN	0.091	0.119	0.154	0.182	0.205	0.234	0.261
	Ensemble	0.08	0.123	0.156	0.182	0.212	0.24	0.266
Sicily	1-L FCNN	0.054	0.076	0.092	0.112	0.118	0.129	0.128
	2-Ls 1D CNN	0.188	0.227	0.231	0.232	0.226	0.219	0.224
	1-L 1D CNN	0.068	0.052	0.055	0.057	0.06	0.065	0.079
	Ensemble	0.098	0.108	0.115	0.119	0.117	0.116	0.108
A.P. Bolzano	1-L FCNN	0.08	0.137	0.173	0.184	0.207	0.215	0.232
	2-Ls 1D CNN	0.141	0.229	0.298	0.367	0.406	0.455	0.462
	2-Ls FCNN	0.172	0.182	0.191	0.197	0.205	0.211	0.217
	Ensemble	0.127	0.18	0.218	0.245	0.268	0.288	0.298
A.P. Trento	1-L FCNN	0.077	0.097	0.113	0.13	0.141	0.154	0.158
	2-Ls 1D CNN	0.21	0.238	0.264	0.294	0.309	0.332	0.353
	1-L 1D CNN	0.081	0.103	0.127	0.145	0.171	0.184	0.196
	Ensemble	0.111	0.136	0.159	0.181	0.199	0.218	0.23
Tuscany	1-L FCNN	0.048	0.053	0.057	0.063	0.072	0.076	0.08
	2-Ls 1D CNN	0.05	0.056	0.073	0.082	0.105	0.085	0.088
	1-L 1D CNN	0.048	0.051	0.055	0.072	0.067	0.076	0.08
	Ensemble	0.03	0.041	0.047	0.056	0.067	0.067	0.069
Umbria	1-L FCNN	0.078	0.128	0.16	0.178	0.224	0.234	0.249
	2-Ls 1D CNN	0.119	0.122	0.119	0.121	0.121	0.119	0.116
	1-L LSTM	0.365	0.391	0.411	0.415	0.435	0.433	0.434
	Ensemble	0.168	0.199	0.217	0.227	0.246	0.248	0.254
Veneto	1-L FCNN	0.066	0.085	0.106	0.134	0.155	0.175	0.185
	2-Ls 1D CNN	0.143	0.145	0.156	0.161	0.17	0.181	0.19
	1-L 1D CNN	0.087	0.094	0.1	0.108	0.125	0.135	0.141
	Ensemble	0.095	0.105	0.119	0.132	0.149	0.161	0.168

A.P. = Autonomous Province, RMSE = Root Mean Squared Error, 1-L FCNN = 1 Layer Fully Connected Neural Network, 1-L 1D CNN = 1 Layer Mono-Dimensional Convolutional Neural Network, 1-L LSTM = 1 Layer Long Short-Term Memory, 2-Ls FCNN = 2 Layers Fully Connected Neural Network, 2-Ls 1D CNN = 2 Layers Mono-Dimensional Convolutional Neural Network.

for each selected architecture and their average ensemble. Figure 4 reports the ground-truth R_t values and the rolling-approach predictions for all the Italian regions and autonomous provinces from 12 to 25 April 2021. The reported predictions were obtained on epidemiological data that were not included in the dataset used for model training, validation, and evaluation procedures, with the aim of evaluating the performance of the proposed framework on new data available daily.

Specifically, for each Italian region or autonomous province, a forecasting phase from 12 to 25 April 2021 exploited the best model (#1) resulting from the K-Fold Cross Validation procedure as well as the average ensemble of the three best-ranked models after the aforementioned procedure. In this phase, the predictions obtained from the models were compared with the W–L R_t estimates (ground truth) and the weekly ISS R_t estimates, as reported in Figure 4. Moreover, Figure 4 also reports the 95% PI of the best model computed as described by Algorithm 2 in Section 2.3. As shown in Figure 4, the proposed framework provides accurate predictions in all the Italian regions and autonomous provinces. In particular, in most regions, all three models proved to be very accurate because their predictions turned out to be very close to the W–L estimates. On the contrary, in some regions, not all models provided accurate forecasts, meaning that some R_t patterns were only learned by specific architectures.

4. Discussion

The present work aimed to define a framework fully based on ANNs for predicting the COVID-19 effective reproduction number R_t in all the Italian regions and autonomous provinces. This indicator is crucial because its value denotes the status of an epidemic (i.e., spread or containment phases), thus allowing policymakers to act in a timely and targeted manner to tackle new outbreaks. In particular, out of the six proposed models, the proposed framework identified three ANNs for each Italian region and autonomous province, which turned out to be optimal for the seven-day ahead prediction of the R_t trend, starting from a daily updated set of 21 epidemiological variables given at a generic day, t. As reported in Section 3, the proposed framework succeeds in the predictive task nationwide, with a prediction error that increases as the prediction time window grows. However, as shown in Table 3, the predictions accuracy varies from territory to territory, meaning that different R_t trends were observed in the Italian country due to the heterogeneous spread of the pandemic. Focusing on Lombardy, which is the region where COVID-19 was more widespread, the predictions obtained through the proposed approach proved to be very accurate, as shown in Figure 3, with a minimum RMSE (on the R_t scale) ranging from 0.035 at day t + 1 to 0.106 at day t + 7. However, because Italian epidemiological data are communicated daily by the ICPD [21], a rolling approach was proposed in order to obtain the most accurate predictions possible at any day to be predicted (Section 2.4), thus improving the prediction accuracy of the models by updating each forecast as soon as new epidemiological data become available. Figure 4 highlights the high level of accuracy of the forecasts produced by the suggested rolling approach for each Italian region and autonomous province. Specifically, it should be noted that, for some regions, not all three models selected for the inference phase produced accurate predictions, meaning that particular trends were learned during the training phase by different architectures. For this reason, the proposed approach does not only rely on a single predictive model, but it also exploits three different ANN architectures in order to maximize the possibility of predicting heterogeneous and complex R_t trends. Notice that the proposed approach allows the forecast of the R_t trend at a daily temporal resolution, whereas official R_t data are currently communicated with a weekly temporal resolution. Forecasts at a finer (daily) temporal resolution could be crucial in a pandemic context, because they could help track the spread of the disease more effectively by rapidly identifying new outbreaks or sudden growth trends. Furthermore, it should be pointed out that traditional statistical methods for estimating the R_t index, typically based on Bayesian approaches, require the observation of epidemiological indicators in a time window whose dimension defines the accuracy of the estimates: the greater the time window, the more accurate the obtained estimates (Section 2.1.1). The proposed framework aims to overcome this limitation by proposing a multivariate approach that exploits the information of different variables at a given day to produce accurate forecasts for the subsequent 7 days. Moreover, it is important to consider that, since 3 November 2020, different containment measures have been applied in Italian regions and autonomous provinces by identifying three pandemic risk scenarios related to increasing critical levels (yellow, orange, and red) [5]. Each risk scenario is based on a weekly evaluation of a set of epidemiological parameters, including the effective reproduction number, R_t, and the number of regional hospitalizations. Therefore, it is important to highlight that the proposed framework could be exploited to predict not only the daily effective reproduction number, but also other epidemiological variables, such as the number of COVID-19 hospitalizations, by simply redefining the target variable in the model training, validation, and testing procedures. In this way, accurate predictions regarding several epidemiological indicators could be obtained in order to facilitate the timely assignment of each territory to the proper risk level in a cost-effective manner. This approach presents some limitations. First, it should be noted that, as expected, an increasing prediction error was observed as the prediction time window increased. For this reason, when possible, the rolling approach should be adopted in order to improve prediction accuracy. When no updated data are available, predictions cannot be updated, leading to a greater prediction error, which varies from territory to territory, as reported in Table 3. Second, consider that the training procedures were conducted by exploiting publicly available data that are regularly revised and updated, thus including an error bias that could be reduced by exploiting more reliable official data. Last, because no official data were available daily regarding the regional effective reproduction number, R_t, the setup of the supervised learning problem required the estimation of the daily ground-truth by means of the Wallinga–Lipsitch algorithm, as described in Section 2.1.1.

5. Conclusions

The present work proposes a framework to forecast the daily $R_t$ index 7 days ahead by exploiting a set of 21 epidemiological indicators at a given day, in each Italian region and autonomous province. The obtained predictions were characterized by a finer temporal resolution (daily) compared with the official $R_t$ data, which were reported weekly by the Italian National Institute of Health. A finer temporal resolution represents a benefit because the $R_t$ is one of the most representative indicators of the infection status. The ANN predictions resulting from the proposed study highlight the fact that accurate forecasts can be achieved for each Italian region and autonomous province, also exploiting the rolling approach procedure that allows the update of the obtained forecasts on a daily basis as soon as new input data become available. For instance, the predictions obtained through the proposed approach in Lombardy, which is a region where COVID-19 was more widespread, proved to be very accurate, with a minimum RMSE (on the R_t scale) ranging from 0.035 at day t + 1 to 0.106 at day t + 7. The proposed framework could be easily rearranged for predicting different epidemiological indicators concerning not only the COVID-19 pandemic, but also other diseases that could potentially appear and spread worldwide. Clearly, this is crucial for containing the disease spread, by imposing timely and appropriate containment measures in order to tackle new outbreaks. As a consequence, the saturation of the health system and the strong impact on the socio-economic fabric could be mitigated because the proposed framework can represent a valid and accurate decision support tool for policy makers.