Modeling the Geospatial Evolution of COVID-19 using Spatio-temporal Convolutional Sequence-to-sequence Neural Networks

Data concerning diseases is commonly available by administrative regions. For many infectious diseases, statistics on the number of cases by countries, municipalities, or other administrative regions are widely available. Thus, the most common approaches to map and visualize the risk of disease is through the use of choropleth maps, where those pre-defined spatial regions are colored in accordance with recorded incidence rates. In the wake of the COVID-19 pandemic, many such analyses have been performed and made available [2, 3, 4]. However, such maps suffer from the intrinsic limitation caused by the arbitrary boundaries imposed by the pre-defined administrative regions considered and exhibit undesirable discontinuities at these boundaries. Geostatistics provides the framework to model spatial variability and to account for data with varying spatial support delivering continuous maps of the risk of disease, overcoming the limitations of choropleth maps. Additionally, geostatistical risk maps can integrate data uncertainty resulting from population size, which is very important in this case, given that the population density varies widely across the country.

A number of methods that estimate the spatial risk of diseases from known datapoints have been proposed [5], including distance-based approaches and kernel density estimation (KDE) methods. These and other methodologies may be used to convert point data to continuous risk maps. However, these methods have been designed to model static or slow-varying distributions, and cannot be trivially adapted to dynamic modeling of diseases. Combined with time-series prediction, KDE has had some success in visualizing disease evolution [6] but depends critically on the accuracy of the time-series prediction method.

Direct estimation of the geospatial risk of COVID-19 in Switzerland using high-resolution Spatio-temporal data analysis [7] based on the modified space–time density-based spatial clustering of application with noise (MST-DBSCAN) algorithm [8] has been proposed. However, high Spatio-temporal resolution health data are rarely available, due to limited information and other limiting factors, such as personal data protection policies.

Our modeling of the geospatial risk model, described in Section 3 is based on geostatistical block sequential simulation [9], assuming a Poisson model for rare diseases as proposed by Waller and Gotway [10] to estimate the local distributions functions [11]. Given the actual municipality-level data, with varying support, this method makes it possible to integrate this uncertainty into the spatial model to generate the incidence rate maps and to assess its spatial uncertainty in any region in the country.

The block sequential simulation we use assumed that the disease propagation is isotropic, at the resolution involved in the simulation. In practice, this may not be true in some scales, since streets, roads, city organization, and other characteristics of built space may privilege some specific directions of transmission. However, this assumption makes the problem tractable and our experiments with mobility information (provided by a cellphone operator and tested as additional inputs) did not contribute to improving the precision of the models, contrary to our initial expectations. The method also assumes that the incidence rate collected at the municipality level is a good proxy for the incidence rate in the municipality since it is used in the kriging process. This is an assumption of the direct block sequential simulation (Block-DSS) method that cannot be avoided. However, we weigh the location of the centroid by the population density of each municipality, in order to decrease the possible negative effects of this assumption.

Compartmental models have been extensively used in epidemiological studies. The first proposals to group populations in different compartments and to use differential equations to model the transitions between compartments, deriving the dynamics of the processes, date back more than a century [12, 13]. Compartmental models use continuous differential equations to model the evolution of an epidemic in a population of individuals, where each individual is assigned a state (a compartment) that determines its stage in the epidemic process [13, 14]. Since then, many articles have been published on the use of compartmental models to predict and simulate the dynamics of epidemics. In the last year alone, hundreds of articles have been published on the application of compartmental models to study the dynamics of the COVID-19 pandemic in many countries and locations, such as China, Italy, and India [15, 16, 17]. A complete description of the many applications of the models, as well as the variants proposed, is outside the scope of this work.

Our own application of compartmental models to predict the dynamics of the pandemic in Portugal is based on the Susceptible-Infected-Recovered-Dead (SIRD) model, which uses four different compartments. Despite the model’s simplicity, many assumptions and approximations are required to make it fit real-world data, including the consideration of time-varying parameters, adjustments for small-sample, high-variance, data, and estimation of unknown values in the observed data series. We used, with adaptations described in Section 4.3, the methods proposed in recent work that used data from the COVID-19 pandemic in New York City, Madrid, and Stockholm, among other cities, together with various states, countries, and regions to estimate a SIRD model [18].

ARMA models have been widely used in the modeling of time series and they can be used to predict the evolution of a univariate stationary stochastic process using two polynomials, one for the autoregression and a second one for the moving average. Proposed by Peter Whittle [19] in 1951, it has been extensively used in time series analysis and represents a solid baseline against which other, more sophisticated, models can be tested.

The ARMA model suffers from a strong limitation since it computes the evolution of the incidence rate at each cell not taking into account the values of any other cells. To circumvent this limitation we also tested a VAR model [20], which computes the incidence rate at each cell using all the values from the other cells. Vector autoregression is commonly used to capture the relationship between different time-varying quantities, as it is a generalization of the (univariate) autoregressive model that does not suffer from the same limitations.

Recurrent neural network (RNN) architectures, based on Long Short-Term Memory (LSTM) units [21] or gated recurrent units (GRU) [22] have been extensively employed in forecasting tasks involving time-series data. However, these models consider the input data as independent sequences of vectors and do not explore the spatial context available in Spatio-temporal data.

To address this limitation, spatial convolution operations were combined with LSTMs, thereby exploiting the abilities of convolutional neural networks (CNNs) and RNNs to effectively model spatial and temporal information, respectively. In the Convolutional LSTM (ConvLSTM) approach [23], all input data structures are 3D tensors, with the first dimension corresponding to either the number of measurements or the number of feature maps and the last two dimensions representing the spatial dimensions (i.e., width and height). By replacing the product operation in the original LSTM with the convolution operation, the future states of a certain cell depend on the inputs and past states of itself and its local neighbors.

Many different architectures have been proposed using ConvLSTM building blocks, and applied to different problems, such as weather prediction from satellite image sequences [24], air quality forecasting [25], and video frame prediction [26].

Inspired by the aforementioned ConvLSTM architecture, Wang et al. [27] proposed another recurrent model, PredRNN, which captures spatial and temporal features in a unified memory pool. In PredRNN, the states of an adapted LSTM cell can travel both vertically between layers and horizontally across states. The authors introduced a new cell state in each LSTM unit, that flows in a zig-zag direction, first upwards across layers and then forwards over time. This extension to standard LSTM units, called Spatio-Temporal Long Short-Term Memory (STLSTM), allows simultaneous flow of both spatial and temporal memory and the standard temporal memory, present in traditional ConvLSTMs. The authors defined the PredRNN as a multi-layer architecture employing ST-LSTMs, having tasked the model with predicting 10 future observations from three distinct datasets, given the previous 10 observations. Results showed that PredRNN utilizing the ST-LSTM architecture significantly outperformed every other baseline, including previous state-of-the-art models. In general, the proposed model was able to more accurately predict the trajectories and digits, particularly in sequences that displayed overlapping sections, when compared to the baseline models.

In the architecture that served as the basis for this work, STConvS2S [28], the authors used a standard encoder-decoder structure comprised of multiple stacked ConvLSTM layers, with each decoder layer initializing its hidden states from the output of the corresponding encoder layer. The final predictions of the model are given by the concatenation of the hidden states from the decoder network followed by a \(1 \times 1\) convolution. This method was adapted to this problem by introducing a number of improvements, including local weights and learnable inputs, described in Section 4.4.

We use the incidence rate, the number of confirmed positive tests in a region, in a period of 14 consecutive days, divided by the population of the region, as a proxy for assessing the risk of infection. In particular, we consider the incidence rate per municipality \(z_\alpha (t)\), defined bywhere \(c_\alpha (t)\) is the number of confirmed positive tests in the 14 days preceding a given day t, at each municipality, \(\alpha\), and \(n_\alpha\) is the respective population size. We will use \(z_u\) to denote the incidence rate at a specific location u (or cell) in the country. In the sequence, we may drop the explicit dependence on t to simplify the notation, when not required.

The daily numbers of new cases by the municipality are reported regularly by the Portuguese DGS, and they have been converted to incidence rates by dividing by the population count of each municipality, as given by official statistics, and scaled to 100,000 inhabitants (thus corresponding to the number of new cases in the most recent 14 days per 100,000 inhabitants). This will be the scale used in the maps and reports in this article.

Since not all cases are detected, this incidence rate may underestimate the real incidence rate. However, it is reasonable to assume that the ratio between the actual number of newly infected persons and the number of detected cases is approximately equal in the different regions of the country, enabling us to use the incidence rate computed in this way as a proxy to the real incidence rate. In practice, we may be underestimating the incidence rate by a significant factor, since it is known [29] that a relevant fraction of the COVID-19 cases are asymptomatic and remain undetected. This has an impact on the value of the parameters of the dynamic models, since important characteristics of the process such as the group immunity threshold, depending on the actual rates of infection and recovery. In practice, we believe this underestimation does not have a strong impact on the accuracy of the models (other than a systematic underestimation), since all model parameters are estimated from actual data and updated as time goes by in order to reflect the dynamics of the pandemic. Therefore, the only significant impact of the fact that not all cases are detected is the underestimation of the incidence rate, by a factor that can be considered approximately constant, if one assumes that there are no significant differences in testing strategy from region to region. This assumption can be justified by the relative homogeneity of the health system in a country like Portugal.

We consider the country divided into a regular grid of square cells, which, in this work, have a dimension of 2 km by 2 km (see Figure 1 for an example of the discretization used).

The methodology can be trivially adapted to a different grid. We compute the incidence rate maps by estimating \(z_u\) for each cell in the country, using geostatistical stochastic simulation, namely block-DSS [30], which is based on a stationary spatial covariance model (i.e., stationary variogram model). This modeling technique allows the daily update of COVID-19 incidence rate maps at high-resolution together with the associated spatial uncertainty. The geostatistical incidence rate maps are continuous (up to the scale of the discretization) and their spatial distribution follows a given spatial pattern as revealed by a variogram model inferred from the data. The incidence rate maps do not show any sharp discontinuities and are not influenced by administrative boundaries, such as municipality borders. Since incidence rates refer to different population sizes, these must be weighted when calculating the experimental variogram such that municipalities with large populations have greater weighting.

Below, we briefly describe this geostatistical simulation method. A full description of the method can be found in previous work by the authors [1]. We assume that the geometric centroid of each municipality \(\alpha\) has coordinates \((x_\alpha ,y_\alpha)\). Data are available for the 278 municipalities in mainland Portugal, as illustrated in Figure 2.

For clarity, in this section, we will drop the explicit dependency on t. The daily update of incidence rates \(z_\alpha\) known for each municipality is assigned to each centroid (weighted by population density) and are used as experimental data in the geostatistical simulation. This is a reasonable approximation given the relatively small area of each municipality when compared with the area of the country.

Within a geostatistical framework \(z_\alpha\) is interpreted as a realization of a random variable \(Z_\alpha\) corresponding to the true, and unknown, incidence rate in municipality \(\alpha\). The expected value of \(E[Z_\alpha ]\) is an approximation of the incidence rate. Z depends on the population size of each municipality. For example, a given municipality with a small population size \(n_\alpha\) (i.e., a small denominator) will have high variance and consequently higher uncertainty in the incidence rate estimate. We use the Poisson kriging model [11] to define the risk varianceBlock-DSS [30] is a stochastic sequential simulation method based on the Poisson kriging model and an adequate modeling technique to deal with data with varying spatial support (i.e., municipalities with varying size and shape) and to make spatial predictions with a change of support, e.g., from the area (block) to point support. In this case, the scale related to the map cells can be referred to as point support, because it denotes a small area when compared with the municipality areas. The following sequence of steps summarizes the block-DSS algorithm [30]:

Block-DSS generates alternative incidence rate maps, designated realizations, at each run, as the random path changes and consequently the conditioning data when simulating a location \({u}\). A given set of realizations provides the value of the uncertainty of incidence rate distribution in a given cell. We have simulated daily sets of one hundred realizations of incidence rates for mainland Portugal.

We computed the incidence values for each cell in the territory, for each day in the period under consideration. The final dataset, used to assess the performance of the models, was generated from data in 278 municipalities, for the period between March 1st, 2020 and February 28th, 2021. The territory was discretized into 40,608 square cells (in a \(288 \times 141\) grid), each cell corresponding to a surface of 2 km by 2 km. The resulting dataset consists of information (incidence and uncertainty estimates) for each cell during the first 365 days of the pandemic. The first 184 days (March 1st until August 31st) were used as the training set, to train or define the parameters of the models, and the following 181 days were used to test the prediction ability of the models.

To illustrate the data, Figure 3 depicts an example of the computation of the incidence rates, for a specific day in the period.

Since the method also derives confidence intervals for the incidence rate in each cell, we also make these available (see example in Figure 4).

The value of the incidence rate at each cell represents a good approximation, given existing data limitations, to the real incidence rate at that location. Both the incidence risk maps and the confidence intervals are made available as supplemental data.

The baseline ARMA(p, q) model is a simple, cell-level, autoregressive-moving-average model [19]. In this model, the incidence rate for each cell u is computed usingwhere \(k_u\), \(\varphi _u(i)\), and \(\theta _u(i)\) are calculated, for each cell u, using linear regression, in order to minimize the observed errors, \(\epsilon _u(t)\).

We used the ARIMA package in Python to compute an independent ARMA model for each cell u. The resulting model was used to predict the temporal evolution of the incidence rate \(z_u(t)\) for each cell u at time t.

We used a VAR model [20], which computes the incidence rate at each cell using all the values from the other cells, in accordance withwhere \(z_u\) denotes the incidence rate and A is an \(N \times N\) time-invariant matrix, calculated, together with \(k_u\), using linear regression, in order to minimize the observed errors, \(\epsilon _u(t)\). We used the VAR package in Python to compute a VAR model for the data and a prediction of \(z_u(t)\), for each cell u at time t.

This approach uses recently proposed techniques for Spatio-temporal modeling using neural networks, resulting in a model derived from the STConvS2S architecture [28], coupled with new normalization-activation layers and learnable parameters intrinsic to each location.

The STConvS2S architecture is a sequence-to-sequence model comprised exclusively of convolutional layers, which are suited to capture spatial features by design (see Figure 5). Convolutions are performed with factorized 3D kernels \(K = t \times d \times d\), where t is the size of the temporal kernel and d is the size of the spatial kernel.

The first component is a temporal block that extracts temporal features through the temporal kernels (i.e., \(t \times 1 \times 1\)), while the second component is a spatial block that receives the output of the previous block in order to extract spatial features through the spatial kernels (i.e., \(1\, \times \, d\, \times \, d\)). Both the temporal and spatial blocks are comprised of successive convolutional blocks with normalization-activation operations, with the temporal block using causal convolutions to maintain temporal coherency during prediction (i.e., predictions for a time-step t make no use of future information from time-steps \(t + 1\) onward).

The feature maps generated throughout the architecture are adequately padded in order to maintain the initial dimensions, and the intermediate feature maps are increased by a factor of 2 in each layer, with the last convolutional layer in both the spatial and temporal block reducing the number back to the initial amount.

Instead of the traditional batch normalization and activation function operations following each convolution, we apply an extension of the normalization-activation layers recently proposed [31]. The chosen operation is an adapted version of the EvoNormB0 layer, which was the best performing batch-based version reported, and explicitly models Spatio-temporal scenarios, considering sequences of two-dimensional inputs when calculating both the batch and instance variance (i.e., the values associated to each time-step in the input sequence are considered separately when computing the variance). Let \(v_1\), \(\theta\), and \(\psi\) denote learnable parameter vectors, and let \(s_{b,d,w,h}\) and \(s_{d,w,h}\) represent the variance of a mini-batch and the variance of a single instance, respectively. This extension, denoted here as B03D, is defined as follows:

Another technique used in our model was the extension of the input data with learnable local features and weights, intrinsic to each spatial location, as originally proposed for the area of wind forecasting, where changes in location imply changes in the model parameters [32]. This technique seeks to improve Spatio-temporal forecasting scenarios by combining both (i) global location-invariant features, and (ii) location-specific features.

Typical CNNs are mostly translation invariant. Translating an input x and convolving with a filter k will yield the same result as translating the feature map resulting from a convolution between x and k. Therefore, standard CNNs treat each spatial location equally and thus learn global patterns. This is often insufficient in many Spatio-temporal forecasting scenarios, where the behavior in specific locations should be guided by their own intrinsic local features.

To allow the learning of such features, we used two complementary techniques, Learnable Inputs (LI) and Local Weights (LW). The LIs correspond to a set of n trainable parameters with the same spatial dimensions as the original input, concatenated with the input before processing with a convolutional layer. These parameters allow the network to learn local features regarding every spatial location, which can complement the learned global patterns or be ignored by the convolutional kernels in situations where they are not beneficial. The LWs further reinforce the individualized learning of different spatial locations, through a locally-connected layer (i.e., a convolutional layer with a different filter at each input region, allowing different spatial locations to be weighed according to their relevance) of weights over the input, resulting in m trainable weights with the same spatial dimensions as the input. Similar to LIs, these weights are afterward concatenated with the original input. In our experiments, both the LIs and LWs are concatenated with the inputs on the channel dimension (sharing the LIs across the different input time-steps and using different LWs in each time-step), prior to every convolutional layer. LIs are implemented with a locally connected layer that receives as input a constant unitary tensor with the same spatial dimensions as the original input, and the LWs are implemented with a separate locally connected layer that receives the original tensor as input.

The ARMA and VAR baseline models used as input data the entire past sequence for each cell up to time t. They were then used to predict the T next days sequentially. We considered only the last value in the output sequence, corresponding to the time-step \(t+T\), and the process was repeated by sliding the window one day into the future and concatenating the ground-truth value corresponding to the time-step \(t+1\) with the existing data.

For the SIRD compartmental models, no training is necessary. As explained in Section 4.3.5, the time-varying parameters for each municipality were determined and used to solve the equations forward in time. For each parameter, we calculate each day in the historical dataset and then use extrapolation to project its values forward T days. By solving the equations for each municipality \(\alpha\) we get a prediction of the future number of new cases for that municipality, \(z_{\alpha }\). We then applied the block-DSS algorithm to determine the incidence rate for each cell in the territory, as described in Section 4.3.5. The application of the block-DSS algorithm used the same spatial covariance matrix as the one used to build the gold standard dataset.

For the STConvS2S architecture, the model received as input a sequence of the previous T contiguous days, and generated a sequence corresponding to the next T days. The final prediction is thus given by considering only the last time-step in the output sequence. This process was repeated by sliding the input window one day into the future, as was done for the ARMA and VAR models. Furthermore, in order to incrementally update the parameters as new data becomes available, we apply a simple online learning procedure, fine-tuning the model with each new input sequence for 5 epochs, after the prediction was made. This enables the model to quickly adapt to sudden changes.

Regarding hyperparameter choices, for the ARMA model, we consider \(p = 7\) and \(q = 1\). For the VAR model, we used \(p=4\). For the SIRD model, we used the pseudo-count parameters \(K=10,\!000\) in the prediction seven days ahead and \(K=100,\!000\) in the prediction 10 and 14 days ahead. These values were selected as the ones that obtained better results. As in the original proposal, the STConvS2S architecture used 3 convolutional layers for both the temporal and spatial block, plus the final convolution to reduce the channel dimensionality back to one. The initial number of convolutional filters was set to 32, and each filter was of dimensionality 5 \(\times\) 5. In order to select the most optimal Local Inputs and Learnable Weights configuration, we varied the filter size k of the LWs between 1 \(\times\) 1 \(\times\) T, 2 \(\times\) 2 \(\times\) T, and 1 \(\times\) 1 \(\times\) 1 (i.e, in this last case, considering a direct element-wise weighing unique to each spatial location and time-step in the input sequence), where T corresponds to the number of time-steps in the input sequence. Furthermore, we varied the number of LIs and LWs \(n \in \lbrace 1, 2, 3\rbrace\). The parameters were finally set as k = 1 \(\times\) 1 \(\times\) 1 and \(n = 2\). Training relied on the AdaMod optimizer using a learning rate of \(10^{-3}\) and batch size of five. Furthermore, for the STConvS2S model, in order to fit in the GPU memory and to avoid losing information over individual regions, the same architecture was trained and tested on three different subsets of the input region, corresponding to the upper third, middle third, and lower third of Portugal. The final prediction results were then concatenated, resulting in a prediction encompassing all of mainland Portugal.

The predictions made by each method were compared with the reference incidence rate values, obtained as described in Section 3. We computed two commonly used figures of merit in order to compare the predictions of the models with the reference data, the RMSE, and the sMAPE.

For a given day t, the root mean square error (RMSE) between the predicted value and the reference value, averaged over all cells, is given bywhere \(z_u(t)\) is, as before, the incidence rate in cell u at time t and \(\widehat{z_u}(t)\) is the value predicted by the model.

Since absolute values of the incidence rate vary greatly by location, it is also relevant to compute the second figure of merit, the symmetric mean average error (sMAPE):

Figures 7 and 8 show the evolution of these two figures of merit for the four methods used (day 0 corresponds to September 7th, 2020). The superior performance of the STConvS2S method is clear, when compared with both the baseline ARMA and VAR models, and the SIRD model. As always, the incidence rate is reported as the number of new cases in the last 14 days per 100,000 inhabitants, and the value of the RMSE uses the same scale.

Table 1 shows the values for the RMSE and the sMAPE of the predictions, averaged over all days in the period September 7th, 2020 to February 28th, 2021. The table includes the four methods studied and also a reference baseline, the ConvLSTM [23]. This table shows the clear superiority of the neural network-based Spatio-temporal convolution methods over the alternatives and, in particular, the superiority of the STConvS2S algorithm, not only in the quality of the predictions but in the stability of these forecasts as the prediction horizon is moved forward.

Figures 9 and 10 illustrates the quality of the predictions for week seventeenth of the predicted period, the week after Christmas (December 28th, 2020–January 3rd, 2021). This was a particularly relevant week since it corresponded to a sharp inflection of the tendency and the start of the third wave. The predictions were made with the data available until December 21st, 2020, for the first day in the week. The window was then adjusted by one day for each successive day. The results clearly show the superior predictive ability of the STConvS2S model which, in this particular week, only exhibited significant error in a sparsely populated region in the south of Portugal, Alentejo, while the other models exhibited a much more significant error in different parts of the country. Other weeks exhibit similar behavior, although there is significant variation in the model errors over time, as shown in Figures 7 and 8.

The results obtained in these tests have shown conclusively the superior predictive ability of the modified STConvS2S method, when compared with the other methods described in Section 4. We attribute this superiority to the ability of this model to take into consideration the incidence rate at neighboring cells, when predicting the evolution of that rate in a given cell. Although, conceptually, the VAR model could also use this information, it has no access to the position of each cell and, even though it beats the basic ARMA method, it does not reach the level of precision attained by the STConvS2S model. STConvS2S also performs significantly better than the baseline Spatio-temporal neural network ConvLSTM. We conjecture that the use of local weights and learnable inputs, coupled with the modified normalization-activation layers, made the STConvS2S more able to adapt to dynamics that change with the characteristics of the territory and, consequently, more precise.

Furthermore, the predictions of the STConvS2S and the ConvLSTM degrade more slowly as the prediction horizon is moved forward, while the predictions of the other methods, namely ARMA and SIRD, degrade rapidly as we ask the models to make predictions further into the future. These results also show that precise predictions for more than two weeks ahead, in this specific problem, are difficult to make, possibly because the dynamics of the phenomenon change in ways that are hard to learn from past data. Still, the Spatio-temporal neural networks (ConvLSTM and STConvS2S) are clearly the ones that degrade more gracefully, although at the cost of significant use of computational resources that, ultimately, make long-term predictions to error prone and too expensive.

The SIRD model suffers from the same limitations as the ARMA model, in that it cannot use information from the neighboring municipalities to infer the dynamics of the pandemic. In principle, this model should have been the one better tuned to the reference data, since it uses exactly the same method to infer cell level incidence rates from municipality level rates as the approach used to obtain the reference data: the Block-DSS algorithm. We attribute the lower precision of this method to two factors: the inherent limitations of the time-varying SIRD models to make predictions far in the future due to the changing pandemic dynamics; and the inability of the SIRD model to take into consideration the geographical relations between municipalities. We conjecture that a modified SIRD model that uses information from neighboring municipalities may improve significantly these results.

The adherence between the predictions and the reference data obtained with the STConvS2S method is dependent on a set of assumptions, analyzed in Section 2.1, which are used by the kriging method used to generate the reference data from the official infection numbers. The models also assume that the number of detected cases is a reasonable and stable proxy for the actual number of infections in a region. We believe that the reference data, which we make publicly available as part of additional material, is relevant for the study of the evolution of the pandemic in Portugal and accurately reflects the underlying incidence rate (and, implicitly, infection risk). However, the methodology used to generate the reference data may be more or less precise in different geographies, where phenomena that violate one or more of the assumptions analyzed in Section 2.1 may occur more frequently.