Forecasting COVID-19 confirmed cases in South Korea using Spatio-Temporal Graph Neural Networks

Abstract

Since the outbreak of the coronavirus disease 2019 (COVID-19) pandemic, a lot of efforts have been made in the field of data science to help combat against this disease. Among them, forecasting the number of cases of infection is a crucial problem to predict the development of the pandemic. Many deep learning-based models can be applied to solve this type of time series problem. In this research, we would like to take a step forward to incorporate spatial data (geography) with time series data to forecast the cases of region-level infection simultaneously. Specifically, we model a single spatio-temporal graph, in which nodes represent the geographic regions, spatial edges represent the distance between each pair of regions, and temporal edges indicate the node features through time. We evaluate this approach in COVID-19 in a Korean dataset, and we show a decrease of approximately 10% in both RMSE and MAE, and a significant boost to the training speed compared to the baseline models. Moreover, the training efficiency allows this approach to be extended for a large-scale spatio-temporal dataset.

1. Introduction

1.1. An overview of the study

On February 6th, 2021, COVID-19 has caused more than 105M confirmed cases and over 2.3M deaths globally [1]. To prepare, understand, and control the spread of the disease, researchers worldwide have come together in a collaborative effort to model and forecast COVID-19 [2]. In the state of pandemic, the ability to accurately forecast caseload is extremely important to help inform policymakers on how to provision limited healthcare resources, rapidly control outbreaks, and ensure the safety of the general public [2].

Forecasting COVID-19 cases is a time series prediction problem in which given 𝑀 observations, we need to predict the most likely number of cases in 𝐻 following time steps.

\(\hat{v}_{t+1}, \ldots, \hat{v}_{t+H}=\underset{v_{t+1}, \ldots, v_{t+H}}{\operatorname{argmax}} \log P\left(v_{t+1}, \ldots, v_{t+H} \mid v_{t-M+1}, \ldots, v_{t}\right)\)       (1)

where 𝑣_t is an observation of the number of confirmed cases at time step 𝑡.

Time series learning approach is popular in forecasting the next number of cases, for example, applying curve-fitting [3], Autoregression (AR) [4], or deep learning [5-7] on time series data. These approaches rely on the information from the time series. They can learn some useful insights or hidden patterns from the time series and use those to forecast next items in the sequence.

However, the development of these time series as well as the epidemic in general is affected by many other factors, such as people movement, weather, geography, etc. Therefore, in this work, we aim to incorporate information from spatial data (geographic location) with temporal data (time series) to make prediction about the number of cases for all locations at the same time. It is very natural and suitable to model our data using the structure of a graph where nodes present regions, edges present the geographic distance among regions, and the temporal features are the features of node through time.

This work is inspired by spatio-temporal graph neural network which is successfully applied in traffic field, for example STGCN [5], DCRNN [6], and the work of Amol Kapoor et al. [2] in COVID-19 forecasting. However, we cannot employ complicated models from traffic field, in which the data is recorded by minutes instead of days like ours. Therefore, we aim to improve the performance from the work of Amol Kapoor et al. [2], with a simpler model than that in traffic [8,9].

We evaluate the effectiveness of the approach when predicting the daily infection cases using a real COVID-19 dataset in 17 provinces of Korea. Based on the previous works, we propose a spatio-temporal graph neural network that can learn dynamic features from both spatial and temporal domains. Subsequently, we use this model to forecast COVID-19 daily new cases for all provinces at once.

Through this work, our contributions are:

• We propose a graph neural networks-based model for COVID-19 forecasting which can incorporate both spatial and temporal data. The architecture allows the model to process all time series simultaneously. With its training efficiency, the model can be extended for large scale dataset.

• Our proposed model is superior to several models, but not as complicated as the one used in traffic field.

• A case study when we apply our proposed method on a real COVID-19 dataset in Korea.

The remaining of this paper is organized as follows. We briefly introduce related works in the rest of this section, including time series forecasting and graph neural networks. In Section 2, we first present our approach, from an overview to details of each component. We then present about experimental setup, including dataset, data preprocessing, and baseline models. The results are shown in Section 3 and discussed in Section 4. Finally, we end up this paper with conclusions in Section 5.

1.2. COVID-19 and time series forecasting

In time series analysis, one of the most consolidated classical statistics-based methods is autoregressive integrated moving average (ARIMA) and its variants [6]. Some studies have been done to apply ARIMA to COVID-19 time series data [10-11]. However, this method is constrained by the stationary assumption of time series and is incapable of dealing with the spatial and temporal dependencies. It is also infeasible to apply a large model for multiple time series. Nowadays, classic statistical methods have been overwhelmed by machine learning models, especially neural networks as higher prediction accuracy and more complicated data representation can be obtained by these models. In this work, we will not consider ARIMA as a baseline model because our effort is not to optimize models for each time series but a model performing well for all time series.

In past few years, recurrent neural network (RNN) [7] shows its effectiveness in dealing with sequence data, such as connected handwriting recognition [8] or speech recognition [14-15]. Long short-term memory (LSTM) [16-17], an artificial RNN architecture, was proposed to process long sequence dependencies and to avoid vanishing or exploding gradients problem. LSTMs are the efficient algorithms to build a time series sequential model and can be applied for COVID-19 data [18-19].

Sequence to Sequence (Seq2Seq) [9] was originally invented for neural machine translation. Attention mechanism from Bahdanau's paper [10] and Luong's paper [11] helped to boost the performance of the Seq2Seq model. Because it can deal with sequence data, it can be adapted to process COVID-19 time series data [2].

1.3. Graph neural networks

Recently, graph neural network has drawn more attention because of the great expressive power of graphs [12]. There are two approaches have been proposed to generalize convolution neural networks (CNN) for graph structure [5]. Those are to apply convolution operation on certain grid forms of rearranged vertices or in spectral domain with graph Fourier transforms [13]. Several studies follow the latter approach to reduce the computational complexity of graph convolution using approximation [25-26].

Spatio-temporal graph neural network has been applied successfully in traffic, such as STGCN [5], DCRNN [6], etc. In traffic prediction on road graphs, the problem is to predict the most likely traffic measurements (e.g., traffic flow or speed) in the following time steps for 𝑛 road segments. By defining the traffic network on a graph and focusing on structured traffic time series, they can predict for short and long terms with relatively low errors.

The work related directly to this study is from Amol Kapoor et al. [2]. They try to incorporate mobility data into the graph. Mobility data, which is generated by GPS of portable devices under permission of people, to some extents, can describe the movement of people into and out of regions through time. This is a crucial information for modeling the disease related directly to human transmission. In this work, we cannot have access to the mobility data of people, hence, we replace it by the location of regions, which is similar to traffic field when using location of sensors.

1.4. Graph convolution network

Convolution on Graphs operator ∗_𝒢 in the spectral domain can be regarded as the multiplication of a signal 𝑥∈ℝⁿ (a scalar for every node) with a kernel Θ.

Θ_∗𝒢𝑥=Θ(𝐿)𝑥=Θ(𝑈Λ𝑈^T)𝑥=𝑈Θ(Λ)𝑈^T𝑥       (2)

where graph Fourier basis 𝑈∈ℝ^n×n is the matrix of eigenvectors of the normalized graph Laplacian \(L=I_{n}-D^{-\frac{1}{2}} W D^{-\frac{1}{2}}=U \Lambda U^{T} \in \mathbb{R}^{n \times n}\) (l_n is an identity matrix, 𝐷∈ℝ^n×n is the diagonal degree matrix with \(D_{i i}=\sum_{j} W_{i j}\), 𝑊 is the adjacency matrix) with a diagonal matrix of its eigenvalues Λ∈ℝ^n×n; filter Θ(Λ) can be interpreted as a function of the eigenvalues of Λ, and is also a diagonal matrix [5]. By this definition, a graph signal 𝑥 is filtered by a kernel Θ with multiplication between Θ and graph Fourier transform 𝑈^T𝑥 [14].

Evaluating Eq. (2) is computationally expensive because of 𝒪(𝑛²) multiplications with eigenvector matrix 𝑈 (graph Fourier basis) and expense of eigen-decomposition of 𝐿 for large graph. Thus, two approximation strategies were proposed to reduce the complexity.Chebyshev Polynomials Approximation To avoid the problem, the kernel Θ can be well approximated by a truncated expansion up to 𝐾^th order in terms of Chebyshev polynomial 𝑇_k(𝑥), as Θ(Λ)≈\(\sum_{k=0}^{K} \theta_{k} T_{k}(\widetilde{\Theta})\) with rescaled \(\widetilde{\Theta}=\frac{2 \Lambda}{\lambda_{\max }}-I_{n}\) (𝜆_max is the largest eigenvalue of 𝐿), and 𝜃∈ℝ^k+1 is a vector of Chebyshev coefficients [15]. The graph convolution equation now becomes:

\(\Theta_{* G} x=\Theta(L) x \approx \sum_{k=0}^{K} \theta_{k} T_{k}(\tilde{L}) x\)       (3)

where 𝑇_k(\(\tilde{L}\))∈ℝ^n×n is the Chebyshev polynomial of order 𝑘 evaluated at the scaled Laplacian \(\tilde{L}=\frac{2 L}{\lambda_{\max }}-I_{n}\). Now, the equation is 𝐾-localized as it is a 𝐾^th order polynomial in the Laplacian. 𝐾 denotes the kernel size of graph convolution, which determines the maximum steps that nodes are away from the central node. By using this 𝐾-localized convolutions, evaluating Eq. (2) costs 𝒪(𝐾|ℰ|) (i.e., linear in the number of edges) as Eq. (3) shows [16].

1^st-order Approximation We now can limit 𝐾=1 in Eq. (3) to make the layer-wise convolution operation linear w.r.t 𝐿 and therefore linear on the graph Laplacian spectrum. As a result, a deeper neural network-based model stacking these first-order approximation graph convolution layers can still recover spatial information without being limited to the explicit parameterization given by the Chebyshev polynomials [17]. Because neural networks are expected to adapt parameters to the change of scale during training, we can further approximate 𝜆_max=2. Under these approximations, the Eq. (3) is simply expressed as:

\(\Theta_{*_{G}} x \approx \theta_{0} x+\theta_{1}\left(\frac{2}{\lambda_{\max }} L-I_{n}\right) x \approx \theta_{0} x-\theta_{1}\left(D^{-\frac{1}{2}} W D^{-\frac{1}{2}}\right) x\)       (4)

where 𝜃₀ and 𝜃₁ are two shared parameters of the kernel. To constrain the number of parameters and alleviate numerical instabilities, 𝜃₀ and 𝜃₁ can be replaced by a single parameter 𝜃=𝜃₀=−𝜃₁; 𝑊 and 𝐷 are renormalized by \(\widetilde{W}\)=𝑊+𝐼_n and \(\widetilde{D}_{i i}=\sum_{j} \widetilde{W}_{i j}\). This leads to the alternative expression of the graph convolution:

\(\Theta_{*_{G}} x \approx \theta\left(I_{n}+D^{-\frac{1}{2}} W D^{-\frac{1}{2}}\right) x \approx \theta\left(\widetilde{D}^{-\frac{1}{2} \widetilde{W} \widetilde{D}^{-\frac{1}{2}}}\right) x\)       (5)

Stacking 1^st-order approximation of graph convolution layers effectively exploits the information from the 𝐾-order neighborhood of a node, where 𝐾 is the number of successive filtering operations or convolutional layers in the model. Furthermore, the layer-wise linear formulation allows us to construct deeper models with a fixed computational budget, because the order of the approximation is limited to one [5], [17].

Generalization of Graph convolutions The graph convolution operator ∗_𝒢 defined on a signal 𝑥∈ℝⁿ can be generalized to a multi-dimensional signal 𝑋∈ℝ^n×c_i with 𝐶_i channels (i.e. a 𝐶_i-dimensional feature vector for every node) as follows:

\(y_{j}=\sum_{i=1}^{C_{i}} \Theta_{i, j}(L) x_{i} \in \mathbb{R}^{n}, 1 \leq j \leq C_{o}\)       (6)

with the 𝐶_i×𝐶_o vectors of Chebyshev coefficients Θ_i,j∈ℝ^K (𝐶_i, 𝐶_o are the size of input and output of the feature maps, respectively). The graph convolution for 2-dimensional variables (2D) is denoted as Θ_∗𝒢𝑋 with Θ∈ℝ^K×C_i×C_o [5].

2. Materials and Methods

In this section, we present our proposed method and experimental setup. For the former part, we show an overview of our method and explain each component in detail in the following subsections. For the latter one, we present how we conduct experiments to evaluate our methods and other baseline models on a real dataset.

2.1. Proposed method

The general architecture of our approach is described in Figure 1. The input is time series of all regions/provinces, and the output is predictions for next time steps of corresponding provinces. In particular, we employ 1D convolution network (1D CNN or Conv1D) and graph convolution to extract temporal and spatial features, respectively. After convolution operators, the output is processed and passed through a fully connected layer to produce final prediction. In the following subsections, we will discuss more details about graph data structure modeling and prominent components of our model.

Figure 1. A visualization of our approach's architecture.

One special point of this design is to process and predict all time series simultaneously instead of building a model for each time series or iterating through each time series.

For clarity, only process for center node is presented, and the same mechanism is applied for all other nodes. Conv1D is used to extract temporal features. At each graph convolution layer, the output of the first Conv1D (𝐻₀, represented in purple) is concatenated and propagated to the next layer. In a graph convolution layer, a node aggregates information from its immediate neighbors (nodes in same color region) to process in the next layer. The output of the last convolution is passed through a fully connected layer to produce final prediction.

2.1.1. Modeling the graph

Before explaining about components of our model, we present how to model the data structure, which is the input to our model.

In this problem, we have multiple time series data which indicate the number of daily infection cases in each location (province). In other words, 𝑣_t in the Eq. (1) can be regarded as a vector (or a matrix with one dimension is 1) containing infection cases for all provinces at day 𝑡. It is natural to model locations as nodes in a spatial graph where an individual node can have connections to others. Time series data can be incorporated as temporal features or temporal edges of these nodes.

In order to model spatial and temporal dependencies, we create a graph with two types of edge. In the spatial domain, edges represent distance between each pair of locations based on coordinates. In the temporal domain, edges simply represent connection to previous days. Temporal edges can be interpreted as sets of connected temporal features through time. The graph-structured data is visualized in Figure 2. Note that, graph notion in Figure 1 is the graph representation in hidden layers, and that in Figure 2 is spatial graph data in each time step.

Figure 2. Graph data structure showing spatial and temporal edges (black and red, respectively).

Each node (province) has spatial edges to others in the same time step (same bounding box) and has temporal edges to itself in previous days (previous bounding box). At time step 𝑡, 𝑣_t indicates the number of cases in all provinces. For clarity, only temporal edges to the center node are shown. In practice, every node in the graph has direct temporal edges to itself in 𝑀−1 previous days. When processing, a graph at a time step will be passed through all steps in Figure 1.

Spatial dependencies are built from the coordinates of provinces. We first computed a matrix of distance which contain pairwise distance from provinces' coordinates. It is then used to compute adjacency matrix 𝑊 in Eq. (5), as same as the work in [6].

where 𝑑𝑖𝑠𝑡(𝑎_i, 𝑎_j) is the distance between two nodes 𝑎_i and 𝑎_j, 𝜎 is the standard deviation of distances and 𝒦 is a threshold for sparsity. We set 𝒦=0.1 for our dataset.

The spatial graph of provinces in Korea and the unnormalized adjacency matrix is shown in Figure 3 and 4, respectively. The normalized adjacency matrix is used for graph convolution and is computed as from Eq. (5).

Figure 3. Korea spatial graph map build from coordinates of each province.

Figure 4. Unnormalized adjacency matrix for spatial graph of provinces in Korea.

\(\hat{A}=\widetilde{D}^{-\frac{1}{2}} \widetilde{W} \widetilde{D}^{-\frac{1}{2}}\)       (8)

To build temporal dependencies, time series data is grouped by 𝑀 days (containing 𝑀−1 previous days) for each node in each time step. 1D CNN is used to extract temporal features from these 𝑀-step time series.

2.1.2. CNN for temporal features extraction

Although RNN-based models are widely used for time series analysis, they still suffer from time-consuming and complicated gate mechanism. On the contrary, CNNs have the advantage of simple structure, fast training, and no dependency constraints to previous steps [5]. In this model, we employed 1D temporal convolution layer (1D CNN) with a kernel size of 𝐾_t.

For each node in the graph 𝒢, the input of CNN can be regarded as a length-𝑀 sequence with 𝐶_i channels 𝑌∈ℝ^M×C_i (𝐶_i=1 for our input time series data). In this scenario, the 1D CNN explores 𝐾_t neighbors in the input series without padding, which shortens the length of sequences by 𝐾_t−1 units each time. Therefore, the convolution kernel Γ∈ℝ^K_t×C_iC_o maps the input 𝑌 to an output 𝑍∈ℝ^{(M-(K_t-1))×C_o}. We apply a non-linear activation 𝜎 (ReLU in our case) to the output as follows.

𝜎(Γ∗_𝒯𝑌)=𝜎(𝑍)       (9)

The CNN can be extended to 3D input by utilizing the same convolution kernel Γ to every node 𝒴_i∈ℝ^M×C_i (i.e., provinces) in 𝒢, noted as Γ∗_𝒯𝒴 with 𝒴∈ℝ^M×n×C_i where 𝑛 is the number of nodes.

2.1.3. Graph convolution network for spatial features extraction

In our model, we use Graph convolution network proposed by Kipf and Welling [17] for extracting spatial features. Feature of each node contains multiple channels for each time step. Hence, we utilize generalization of graph convolutions to aggregate information from neighbor nodes for each channel and in each time step, as shown in Figure 1. Stacking two layers of graph convolution allows a node to aggregate information from its 2-hop neighbor nodes.

In particular, the input of case prediction contains 𝑀 time steps of province graphs as Figure 2 shows. Each time step 𝑣_t can be regarded as a matrix 𝑋∈ℝ^n×C_i, i.e., a matrix of 𝑛 𝐶_i-dimensional vectors for 𝑛 nodes (in this case, 𝐶_i=1 for the input data). For each time step 𝑡 of 𝑀, the same kernel Θ of graph convolution operation is applied on 𝑋_t∈ℝ^n×C_i in parallel. Hence, the graph convolution can be further extended for 3D signals, denoted as Θ_∗𝒢𝒳 with 𝒳∈ℝ^M×n×C_i.

In our model, graph convolution is not applied directly to the input 𝒳 but to the processed data \(\tilde{x}\)∈ ℝ^{M-(K_t-1)×n×C_o} which is the output of the 𝐾_t-kernel CNN. However, the calculation is unchanged as the meaning of each dimension is similar.

2.1.4. Skip-connection model

Inspired by the work of Kapoor et al. [2], we utilize skip-connections between layers to avoid diluting the self-node feature state. In particular, the output of each graph convolution layer is concatenated with the output of the first temporal convolution as supplemental information from temporal node features. Our proposed model can be formulated as:

𝑯₀=𝜎(Γ_∗𝒯𝒳)=𝜎(Γ_∗𝒯(𝑣_t|𝑣_t-1|…|𝑣_t-M+1))       (10)

𝑯_l+1=𝜎(Θ_∗𝒢𝑯_l)|𝑯₀       (11)

𝑯_o=𝜎(Γ_∗𝒯𝑯_s)       (12)

\(\hat{\boldsymbol{H}}_{o}\)=𝑯_o.reshape(𝑯_o.shape([0], 𝑯_o.shape[1], −1)       (13)

𝑷=𝐹𝐶(\(\hat{\boldsymbol{H}}_{o}\))       (14)

where 𝒳 represents for the input data containing of temporal features through 𝑀 time steps. The shape of tensor 𝒳 is (batch_size, n_nodes, n_timesteps, n_features) in which batch_size is batch size, n_nodes is the number of nodes (provinces) in the graph, n_timesteps is the number of time steps (𝑀 for the input), n_features is the number of features of the data (1 for the input). 𝑯 represents for hidden state at layer 𝑙 with shape of (batch_size, n_nodes, n_timesteps, n_features). Note that, n_timesteps and n_features are different after each layer of our model. Γ_∗𝒯 and Θ∗_𝒢 are convolution operators for extracting temporal features and spatial features (on graph), respectively. | is concatenation and 𝜎 is non-linear function (ReLU for our case).

According to Equations 10 – 14, the first hidden state 𝑯₀ is the output of a CNN on temporal input 𝒳. This hidden state is then concatenated with the output of each graph convolution to create skip connections. After 𝑠 graph convolution, the hidden state 𝑯_s is passed through another CNN to extract temporal information and produce hidden state 𝑯_o. The number of time steps is decreased by 2∗(𝐾_t−1) (𝐾_t is the kernel of CNN) because there are two CNN layers. The hidden state 𝑯_o is reshaped to flatten the time and feature dimensions. The final prediction 𝑷 having the shape of (batch_size, n_nodes, n_predicted_timesteps) is the output of a fully connected layer (𝐹𝐶) applied on the last hidden state \(\widehat{\boldsymbol{H}}_{o}\) (n_predicted_timesteps is the number of time steps ahead to forecast). A visualization of these equations is presented in Figure 1.

2.2. Experimental setup

In this section, we present preparations to conduct experiments, including dataset, data preprocessing, and architecture of models.

2.2.1. Dataset

We verify our model on a real COVID-19 in Korea dataset which is reported by Statistic Korea using API provided by the public data portal site [31-32]. It includes the number of daily infection cases from April 1^st, 2020 to January 12^th, 2021 in all 17 provinces of Korea. We trained all models on 70% of data (until October 20^th) and evaluate the performance on the other 30%.

In addition, we also utilize the geographical location (coordinates) of all provinces in Korea extracted from COVID-19 in Korea dataset to construct the spatial dependency among nodes in our graph. This dataset is available on Kaggle [33].

2.2.2. Data preprocessing

In our time series dataset, we first transform our data into stationary time series. Although many deep learning models do not need to assume a stationary input data, the empirical results of this COVID-19 in Korea dataset suggest transforming it into a stationary form. In this case, the trend can be removed from the series, and added back to forecasts later to return the predictions to the original scale and calculate comparable metrics.

A standard way to discard non-stationarity in the data is by differencing. The difference value is the result of subtraction between the observation of the current time step 𝑡 and the previous observation 𝑡−1. The transformation gives us a different series, in other words, the changes to the observations from one time step to the next.

Subsequently, we re-scale the stationary data using standard scaler of scikit-learn, which has the formulation

\(z=\frac{x-\mu}{\sigma}\)       (15)

where 𝜇 is the mean of the training samples, and 𝜎 is the standard deviation of the training samples.

2.2.3. Baselines

For our experiments, we try to forecast one time step ahead using 𝑀=7 previous observations.

Persistence model (PM) considers that the number of cases at 𝑡+1 is equal to that at 𝑡. It is also called the naive predictor.

LSTM baseline model contains a stack of two LSTM layers with 8 units each and a final dense layer. The LSTM layers learn sequential information from the input through recurrent neural network. The fully connected layer takes the final output from the second LSTM layer and produce a vector of size 17, which is equal to the number of time series (or provinces) in the data.

Seq2Seq model has an encoder-decoder architecture. The encoder comprises of a fully connected layer (with 4 units) and a GRU layer [18] that can learn from a sequential input and return the full output sequence and the last state. The decoder has an inverse structure of the encoder. It takes the output sequence and last state from the encoder and applies Bahdanau attention [10] for each decoding step to make predictions for target sequences. The number of units for both GRU layers is 4 and the decoder's final dense layer is 1 because we are dealing with a regression problem. We use a single Seq2Seq for all time series in our data.

For both LSTM and Seq2Seq models, we use MSE loss, Adam optimizer [30] with a learning rate of 0.01, and a droprate of 0.2 for training. We implement these models using Tensorflow framework.

As mentioned earlier in Section 1.2, we will not consider ARIMA as a baseline model. We will also not consider STGCN [5] here because this model is more complicated and requires much more data.

2.2.4. Proposed method setup

For our proposed model, thanks to the parallel and combination of extended CNN and graph convolution, our model is able to process multi-channel input and arbitrary batch size. The first 1D CNN turns 1-channel time series into 32-channel graph features. The next two graph convolution layers have 16 output channels each. The second 1D CNN also has 16 output channels and the final dense layer has one unit which is equal to number of time step ahead to predict. We train the model for 30 epochs with a batch size of 16, Adam optimizer [30] with a learning rate of 0.001, MSE loss, and a droprate of 0.2. We implement this model using PyTorch framework for flexibility.

We also implement STGNN model from the work of Kapoor et al. [2]. There are two main differences between our model and STGNN. Firstly, instead of 1D CNN layers in Figure 1, STGNN simply uses multilayer perceptron (MLP) for the first and the last layer with 4 and 1 unit, respectively. Secondly, STGNN utilizes implementation of graph convolution layer from Kifp and Welling [17], which can only process 1-channel signal and one graph at a time. For two graph convolution layers, the graph output features (before concatenation) are 1- channel signal with 4 output units each. The model is trained for 30 epochs with Adam optimizer [30], a learning rate of 0.001, MSE loss, and a droprate of 0.2.

A summary of several hyperparameters is presented in Table 1. The hyperparameters were carefully chosen for each model. Note that, arbitrary batch size is only for Seq2Seq and our model. In LSTM model, we use stateful LSTM which requires to define batch size in advance. Thus, the batch size of 1 guarantees the model to work with any size of training and test set. Meanwhile, as mentioned above, STGNN can only process 1 graph at a time.

Table 1. Summary of hyperparameters.

We run all the experiments for 5 times and report the mean and standard deviation of the metrics for each model. The metrics we use in this study are root mean squared error (RMSE) and mean absolute error (MAE). The experiments are conducted on a computer using a NVIDIA GeForce GTX 1050 graphic card. All source code is available on Github [34].

3. Results

In this section, we provide the results of the experiments for evaluation. Concretely, we report the metrics and running time, and visualize the prediction.

3.1. Metrics evaluation

Table 2 shows the performance of our proposed model and baselines in terms of RMSE and MAE on the test dataset. It shows the mean and standard deviation of the metrics through 5 runs. It can be seen that graph neural network models (GNN) outperform baselines, reducing the errors by around 10%. Between two GNN models, our model produces the lower RMSE and MAE error (17.8492 and 8.2180 in comparison to 18.1490 and 8.2579).

Table 2. Mean ± Std of RMSE and MAE through 5 runs of each model.

For a detailed view of the errors for each province, we report in Table 3 and 4. GNN models are able to reduce the errors in all time series prediction. Although the performance of 2 GNN models are relatively equivalent, our proposed model is likely to reduce the errors in some extreme cases in which the errors are extremely high.

Table 3. Mean ± Std of RMSE through 5 runs for each province prediction of each model.

Table 4. Mean ± Std of MAE through 5 runs for each province prediction of each model.

For visualization purpose, we select to plot time series of three provinces which have the 25-percentile, medium and 75-percentile of the errors. Figures 5 – 7 show the ground truth of daily number of infection cases and models' forecasts on the test set in Daegu, Jeonbuk, and Ulsan, respectively. According to the figures, Seq2Seq model cannot perform well when it even makes negative prediction. This is obvious because one single Seq2Seq model cannot deal with different patterns which are featured for each time series. Both LSTM and Seq2Seq predictions are about 1 day lagging the ground truth, thus, their performances are similar or worse than the PM model. Meanwhile, our model makes more stable predictions. It can follow the trend of the ground truth when the numbers are stable but lag a bit behind when the ground truth fluctuated wildly. Our model's numbers do not change as significantly as those of the others.

Figure 5. The ground truth and predictions (mean and 95% confidence interval) of the number of cases in Daegu.

Figure 6. The ground truth and predictions (mean and 95% confidence interval) of the number of cases in Jeonbuk.

Figure 7. The ground truth and predictions (mean and 95% confidence interval) of the number of cases in Ulsan.

The last thing to note is that in December 2020, Korea witnessed the third wave of infection. It was more serious than ever before when the number of daily cases continuously peaked. That makes it harder for the models to provide accurate predictions.

3.2. Training efficiency

Turning to traiing time in Table 5, it respectively takes LSTM and Seq2Seq approximately 4.5 and 5.3 seconds to finish training an epoch. For our model, only about 0.04 seconds is needed for an epoch with a batch size of 16, accelerating about 100 times.

Table 5. Elapsed time for training of each model.

4. Discussions

As shown in Table 2, our model and STGNN were able to reduce the errors comparing to other baselines. This means that, to some extents, the spatial data does impact on the behavior of the temporal data. Because other baseline models only relied on the time series data. Meanwhile, our model and STGNN also modeled distance relationship among provinces as spatial dependencies among these time series.

As a result, spatio-temporal graph is an effective method to combine temporal and spatial data together. In this research, we can utilize time series data of all provinces and the distance relationships among these provinces to forecast the following time steps. For further study, we can use different data to construct temporal and spatial dependencies, such as transportation, temperature or human activities which can have impacts on the spread of the disease. Moreover, this work can be also adapted for other spatio-temporal data if we can define spatial and temporal domains.

Between two GNN models, our model is superior to previous STGNN as it further reduces the errors Table 2. This probably comes from the temporal convolution layer (1D CNN or Conv1D). Particularly, we employ a CNN layer instead of a MLP because a CNN layer can explore surrounding time steps better than a normal fully connected layer. Another advantage of our model over STGNN is the processing capacity. Because of the parallel and combination of extended CNN and graph convolution, our model can process multi-channel input and arbitrary batch size. Therefore, we can also incorporate more input information by increasing the number of channels of the data.

As illustrated in Table 5, our model accelerates the training speed by approximately 100 times comparing to RNN based models due to couple of reasons. Although there are some differences in the implementation, we are also beneficial from the model design. First, we utilize convolution instead of recurrent neural network which is simpler and faster, as mentioned in subsection 2.1.2. The model design also allows us to train and predict all time series simultaneously, which save us from iterations. The training efficiency shows the potential scalability of our model for large scale dataset.

5. Conclusions

Forecasting daily COVID-19 infection cases for each region or province is one of crucial problems to help the combat with the pandemic. In this study, we developed a graph neural network-based approach for this forecasting problem in which we integrated both spatial and temporal data, rather than only time series as other RNN-based models. Specifically, we utilized both time series data of all provinces and the distance relationships among these provinces to construct spatial and temporal dependencies of the input data. We then applied our proposed model to learn both temporal and spatial features and forecast the following time step. The model's design allows us to process all time series simultaneously. We evaluated the performance of our model on the real COVID-19 dataset in Korea and showed that our model outperformed other RNN-based models, reducing the errors by 10% and cutting down the running time significantly. For future work, further processing can take more data into account, such as transportation, temperature or human activities which can have impacts on the spread of the disease. Moreover, the training efficiency opens potentials for our model to be extended for large scale spatio-temporal data or adapted for other fields.