Statistical Inference Using GLEaM Model with Spatial Heterogeneity and Correlation between Regions

Compartmental model over multiple regions

In GLEaM ([1]) and other works ([26, 27]) involving transportation, the whole area is first divided into subdivisions. For example, the GLEaM model divides the total area of 220 countries into over 3300 subpopulations centered around major airports and [26] divided Milwaukee County and Dane County in the state of Wisconsin into several regions. In this paper, we consider abstract subdivisions in the whole area, which will be referred as “regions” hereinafter until further specifications in the later experiment sections. We denote by n the number of regions.

In our model, we use continuous time t ∈ [0, T ], where it is assumed that the evolution of the epidemic lasts within a period of T time units. The unit of time is fixed as one day throughout this paper. Note that when we introduce the transportation model in below the traveling matrix is assumed to be constant within each day, and also the observed data is collected on a daily basis. Thus we will use notation of discrete time (days) from 1, …, T hereinafter, however, the evolution dynamic itself is modeled over continuous time.

For each region, we consider the following epidemiology compartments adapted from the SEIR model:

• S_k(t): Susceptible.

• E_k(t): Exposed and infectious.

• H_k(t): Hospitalized.

• R_k(t): Removed.

The subscript k of the states denotes that they belong to the k-th region, and dependence on the continuous time t is addressed through expressing the states as functions of t ∈ [0, T ]. Removed individuals are those who have recovered or died.

At time t, we use N_k(t) = S_k(t)+E_k(t)+H_k(t)+R_k(t) to denote the total population in the k-th region. N_k(t) is allowed to be time-varying since the inter-region mobility is taken into consideration, especially for the days before the implementation of travel restrictions. However, since the traveling volume is not comparable to the total population in a region, the fluctuation of the total population in a region is not obvious. In this paper, we assume that keeps constant over time, which means that we consider a closed system in this paper, where exported/imported cases are not considered. However, it is worth noting that we do allow the transportation of active virus carriers between regions within our system. We remark in advance that this assumption is reasonable for the real-world data sets considered in this paper. From January to February 2020, strict international travel restrictions were imposed in China. While for data in Europe, from May to August 2020, the local spread of the epidemic has reached a relatively high level, and the imported cases were not comparable to the indigenous cases. We also denote (N_a)_k(t) = S_k(t) + E_k(t) + R_k(t) as the total population that are permitted to move in the k-th region, excluding the hospitalized ones.

Transportation between regions

Transportation plays an essential role in the spread of COVID-Actually, [31, 32] indicated that the travel restrictions were remarkably important in mitigating the transmission of COVID-19, especially in the early stage of the pandemic. Recently, as detailed in [33], the Omicron variant had spread to 110 countries and had become dominant in many of them by 22 December 2021, only one month after its first report from South Africa on 24 November 2021. This motivates us to also take transportation into consideration in this paper. Motivated by the GLEaM model [1], we introduce the transportation between regions. Specifically, we denote (w_l)_ij as the traveling volume from region i to j on the l-th day (l = 1, …, T). Then the traveling matrix W_l on the l-th day can be written as

Then, inspired by [34], we consider the transition dynamic of ξ(t) = {S₁(t), E₁(t), H₁(t), R₁(t), …, S_n(t), E_n(t), H_n(t), R_n(t)} ∈ {0, 1, 2, …, N }⁴ⁿ over t ∈ [0, T ] described below:

• Transmission in the k-th region: A case from E_k(t) chooses an individual from N_k(t) randomly at Poisson rate λ_k(t), and the individual chosen is infected if it is of state S_k(t). Note that this is different from traditional SEIR models, since we assume that for the COVID-19 case, the pre-symptomatic patients from E_k(t) can be contagious.

• Hospitalization in the k-th region: Each individual in E_k(t) will be hospitalized with Poisson rate δ.

• Recovery or death in the k-th region: Each individual in H_k(t) will transfer into R_k(t) with Poisson rate γ_k. The rate γ_k owns spatial heterogeneity due to the uneven distribution of medical resources.

• Transportation between regions: At the end of the l-th day, all the individuals in region i except the ones in H_i(l) have the same probability to travel from the i region to the j region, and the total traveling volume from the i region to the j region is (w_l)_ij. We assume that there are no transmissions happening during the transportation between regions. If we denote as the number of people transported from {S_i(l), E_i(l), R_i(l)} to {S_j(l), E_j(l), R_j(l)} at the end of the l-th day, then follows a multinomial distribution. Specifically, with .

As a consequence, N_k(t) is a piece-wise constant function of t which only changes at the end of each day. Specifically, for any time t ≥ 0,

Differential equation with spatial heterogeneity

Following [35, 36], we consider the mean-field differential equation system defined by (2.4), which is a deterministic version of the stochastic model above: where is the deterministic counterpart of S_k(t) and so on. Note that since the traveling volume is a piece-wise constant function of t, evolves as a piece-wise linear function of t, and coincides with N_k(t) when t takes integer values. and are the accumulated confirmed and removed cases determined by the ordinary differential equation (ODE) system (2.4) in the k-th region at time t respectively. (or only takes the increment of (or into account, meanwhile guaranteeing that each case is only accounted for once.

Furthermore, we assume that the accumulated confirmed and removed cases are available from data on a daily basis, which are denoted as and , respectively. We also assume that R_k(0) is 0, while E_k(0) and H_k(0) are left to be inferred for each k = 1, …, n. For inference of parameters, we further denote as the newly confirmed cases on the i-th day (k = 1, …, n, i = 2, …, T) determined by (2.4), and (ΔC_a)_k(i) as the newly confirmed cases computed from data, namely (ΔC_a)_k(i) = (C_a)_k(i) − (C_a)_k(i − 1). The same convention holds for the definitions of and (ΔR_a)_k(i). Note that the data are random in nature.

Note that the model and the inference algorithm described below can be applied to estimate spatially heterogeneous transmission parameters as long as , and are available. The availability of and is required in many works that use the SEIR model to estimate transmission rates of the epidemic [37, 38, 39], and the transportation network is also used in GLEaM [1] and its applications. However, in contrast to the works based on GLEaM [12, 4, 5, 6, 15], we allow parameters to possess both spatial and temporal heterogeneity, and further utilize the correlation between regions in the inference of parameters. We remark that the spatial and temporal heterogeneity is reflected in the fact that the transmission parameters {λ_k} are allowed to vary in both space and time in our model. The temporal heterogeneity is introduced in more detail for the real-world data in Section 4.2.

Step 1: Estimate

We first estimate by maximizing over γ_k for each k.

By assuming that the newly removed cases in one day follow a Poisson distribution whose mean equals to the product of γ_k and the accumulated hospitalized cases (which is the difference between the accumulated confirmed cases and the accumulated removed cases, and thus is observable) the day before, we have that

Then, we estimate for each k sep-arately.

Step 2: Estimate and the marginal posterior distributions of

Next, we estimate the remaining parameters as follows.

First, note that λ_k is the transmission rate, which reflects the infectivity of the virus, the prevention and control measures, and the rates of close contacts in the k-th region.

Notice that the ODE system (2.4) is the mean-field version of our stochastic model, and (especially are determined by the parameters, thus by the Markov property, {(ΔC_a)_k(i)} are all independent for k = 1, …, n, i = 1, …, T conditioned on the parameters . Furthermore we suppose that . Thus, the likelihood of can be written as

Now, we denote the posterior distribution of as . Then by Bayesian formula, where is the prior distribution of to be determined and Z = p {(ΔC_a)_k(i)}_k,i is a constant irrelevant to . We further denote , then

To fit the realistic evolution of the epidemic more precisely, is estimated as arg min with reasonable prior distribution . Then, MCMC sampling scheme starting from is applied to get the posterior distribution for . This process might possess higher computational efficiency than choosing the initial point for MCMC randomly or empirically.

Choosing prior distribution

The remaining problem is to choose the prior distribution . Presuming that the transmission rates in the regions owning more similarities are closer, is designed to combine the information of correlations between regions. In particular, given a matrix A which characterizes the pairwise similarities between the regions, and if we denote λ = (λ₁, …, λ_n)^T ∈ ℝⁿ, where D = diag{d₁, …, d_n} is the degree matrix of A with is a constant depending on σ and A. Here, a small σ is chosen for to be a probability measure without imposing much restriction on λ. Then,

The parameter estimation procedure could be extended to the case when {λ_k} are time-varying by modifying (2.11), which we will introduce in more detail in Section 4.2 for the real-world data in China and Europe.

A is constructed from affinity matrix W by further addressing the correlation between regions with more similarities. To treat different data sets and W with a unified approach, we assume that max_i,j W_ij = 1 (W can be re-scaled entry-wise if necessary).

For a given W = (W_ij)_{i,j∈{1,…,n}} whose choice is detailed later, the next step of attaining A is to divide the n regions into d groups (, and ∀i ≠ j, D_i ∩ D_j = ∅) where the regions in the same groups have more similarities. Then, for given β ∈ (0, 1) and a given penalty factor µ > 0, a_ij is constructed as follows:

We remark that β in (2.10) is taken to be 0.1 for all the experiments in this paper. By constructing A as in (2.10), correlations for regions in the same groups are further strengthened, whose transmission parameters are imposed with stronger restrictions.

Now we specify the choice of W for data sets that will be analyzed later in this paper. For simulated data and real-world data in China, in which cases the transportation data are available, we construct W from the traveling volume matrices . Specifically, , where . Nevertheless, for real-world data in Europe, where we are not aware of traveling data publicly available that are sufficient for the proposed model, W is just the all 1 adjacency matrix. We remark that the affinity matrix W may also be obtained by ways other than using the transportation data, as long as it reflects the similarities between districts.

From the definition of a_ij in (2.10), the inference of can be equivalently written as follows it can be seen that by choosing the prior distribution as in (2.8), a l₂ regularization term is imposed for better generalization.

Finally, after choosing determined by A and σ, the optimization process arg min is accomplished through Matlab function fmincon. In addition, the posterior distribution for could be obtained by classical MCMC sampling scheme starting from as mentioned above.

We summarize the procedure described in Section 2.3 in Algorithm 1 below.

Algorithm 1

GLEaM model with Graph Laplacian Prior Distribution

• Download figure
• Open in new tab

Experimental settings

The results in Section 3 are for 100 replicas. In each replica, three random trajectories are sampled independently according to the stochastic model with prefixed parameters, part of which are treated as the ground truth training, validation, and testing trajectory, respectively (see more details in Appendix B.1.1). For each model, we first fit the parameters using the training trajectory and then predict the testing trajectory using the estimated parameters. Note that we detail the choice of hyper-parameters for the model proposed in Appendix B.2. In particular, the penalty factor µ > 0 is chosen by cross-validation and chosen as the value minimizing the validation error, since for the simulated data, the validation error are identically distributed as the testing error. The comparison of trajectory prediction is from one typical realization, for which we compare the ground truth training and testing trajectories with the fitted training and predicted testing trajectories for all the models. Additionally, parameter estimation and quantitative evaluations are compared with mean and standard deviation over all 100 replicas. The detailed computations of training, validation, and testing errors can be found in Appendix C.

Sensitivity analysis of σ

Note that parameter inference with the proposed model involves the parameter σ, as shown in (2.11). Therefore, the sensitivity analysis for the parameter σ in (2.11) is performed for the four provinces case. Specifically, the results for σ varying from 10⁻⁶ to 10⁰ are presented and compared. Details can be found in Sections 3.2.4 and 3.2.5. Similar results are obtained for other data sets, and details are omitted.

Mismatched partition of regions

We also note that the graph Laplacian penalty of the proposed model depends on the partitioning the regions into several groups, as described in Section 2.3. Since the graph knowledge is usually not fully known, it is a question whether our methods can still perform well without accurate prior knowledge. For the thirty provinces case, we report the results of the proposed model with a mismatch between the partition of the regions and the ground truth division, the details of which can be found in Section 3.3.

3.2.1 Data Description

In this simulated study, we let n = 4, T = 20, and set the threshold T_th separating training and testing data to be 10 (more detailed can be seen in Appendix B.1.1). The other prefixed parameters are listed below:

• For k ∈ {1, 2, 3, 4}, N_k(0) = 10⁶, E_k(0) = 30, H_k(0) = 10.

• For l ∈ {1, …, T }, i, j ∈ {1, 2, 3, 4} and i ≠j, (W_l)_ij = 5 × 10³.

• λ₁ = 0.5, λ₂ = 0.47, λ₃ = 0.4, λ₄ = 0.37, δ = γ₁ = … = γ₄ = 0.14.

The four provinces are divided into two groups, with the first group consisting of province 1 and 2 and the second group consisting of province 3 and 4. The similarities within groups are reflected in the settings that the values of {λ_k} are closer for provinces in the same group.

3.2.2 Models to Compare

We first specify the models to be compared below. The last one is the proposed model, and the first four models serve as baselines with different settings.

1. The model with uniform prior distribution, without heterogeneity or migration.

2. The model with uniform prior distribution, without heterogeneity but with migration.

3. The model with uniform prior distribution, with heterogeneity but without migration.

4. The model with uniform prior distribution, with both heterogeneity and migration.

5. The model with prior distribution based on graph Laplacian, with both heterogeneity and migration.

For better illustration and comparison between the models in the experiment results, the Models 1-5 are summarized in Table 2 below.

First, the models with uniform prior distributions themselves (Models 1-4) are compared according to whether two key assumptions exist in the model:

1. Whether the transmission rates {λ_k} are allowed to vary over regions.

2. Whether there exists transportation between regions.

Then, the model with prior distribution based on graph Laplacian (Model 5) is compared with those using uniform distributions as prior distributions (Models 1-4). The former one utilizes the correlation between subpopulations by adding a l₂ regularization term for the model. In contrast, only lower and upper bounds are imposed on parameters without other prior information being used in the latter ones.

For Model 5, we conduct the sensitivity analysis for parameter σ in (2.8) and present results for σ= 10⁻⁶, 10⁻³ and 10⁰ respectively. Then for Model 5, by (2.11), are estimated by the optimization problem in (3.1): where β is taken to be 0.1.

We remark that in Models 1-4, the estimation of still follows a similar two-step procedure as in Model 5, and the first step of obtaining remains formally the same. The difference lies in the optimization object of . First, the l₂ regularization term becomes prior knowledge of the parameters’ upper and lower bounds. Second, for models without heterogeneity of parameters, are forced to be the same in (2.4). For models without transportation between regions, terms involving W_t disappear in (2.4). Additionally, for the other data sets considered in the following sections, the parameter estimation methods for the baseline models are similar and thus will not be repeated.

Finally, note that the model in [12] is similar to Models 1 and 2, since they all assume a spatially homogeneous transmission parameter. However, [12] assumed that the epidemic was seeded from one seed region while Models 1 and 2 do not make such assumption. Moreover, [12] focused more on the spread of the epidemic from the seed region at the early stage of the pandemic, and only the introduction dates in the other regions were utilized for the estimation of transmission parameters. In comparison, the estimation of transmission rate in Models 1 and 2 exploits the data in all regions in the whole process.

3.2.3 Results of Trajectory Prediction

First, we remark that in Model 5, µ is chosen to minimize the mean of validation errors over 100 replicas. The weighted (simply averaged) validation errors, and and , are defined as in Appendix C. We remark that the superscript ^[Val] refers to when the error is computed on validation data, and the subscripts _(w) and _(s) denote that the errors are the weighted and simple average of relative errors over time respectively.

The mean of weighted validation errors and are shown in Figure 1 below, and the mean of simply averaged validation errors and are shown in Figure A1. µ is chosen to be 10^2.7, which minimizes the mean of validation errors, as can be seen from Figure 1 and Figure A1.

• Download figure
• Open in new tab

Figure 1:

Mean of weighted validation errors and against µ (computed by (A5)) over 100 replicas for simulated data (four provinces). The horizontal lines show the values of the mean of the validation errors when µ = 0.

Then, The trajectories of a typical realization are plotted in Figure 2 and the differences between the true and fitted trajectories are shown in Figure 3. As shown in these two figures, heterogeneity helps improve the prediction of testing data more than transportation, while introducing migration without heterogeneity of parameters worsens the estimate as can also be noticed from Table 4. More explanations can be found in Section 3.2.5.

• Download figure
• Open in new tab

Figure 2:

True and fitted trajectories of Models 1-5 in the four provinces. The vertical lines show the threshold T_th = 10.

• Download figure
• Open in new tab

Figure 3:

Absolute errors of fitted trajectories of Models 1-5 in the four provinces. The vertical lines show the threshold T_th = 10.

Additionally, Model 5 with prior distrbution based on graph Laplacian lowers the absolute errors of predicted trajectories compared with Model 4.

As shown in Figures 2 and 3, Models 1 and 2 have slightly better generalization accuracy than Model 5 for Province 3. On the one hand, for data in this replica, the estimated λ₃ is 0.3860 using Model 5 and 0.4640 using Models 1 or 2 (recalling that the ground truth λ₃ is 0.4). On the other hand, due to the randomness of the generated testing data, the sampled newly confirmed cases are much more than the deterministic ones in Province 3 obtained by running (2.4) with the ground truth parameters. Hence, although all these estimates of λ₃ are biased from the ground truth 0.4, estimates using Models 1 and 2, which are biased up, lead to less absolute errors.

3.2.4 Results of Parameter Estimation

The mean and standard deviation of {λ_i} estimated by the five models for four provinces case are reported in Table 3. It can be observed that models allowing heterogeneity estimate parameters more accurately, and Model 5 that integrates the correlation leads to slightly better estimate for λ₂. We can see that compared to Model 4, the estimates of smaller λ_k’s (such as λ₂, λ₃, λ₄) become larger, and the estimate of λ₁ which has the largest value becomes smaller, since the graph Laplacian penalty tends to make closer to each other.

Moreover, we performed a sensitivity analysis for the hyper-parameter σ to check that the results are robust to σ. The last two rows of Table 3 show the parameter estimation results for Model 5 with σ= 10⁻³ and 10⁰ respectively (more values of σ ∈ [10⁻⁸, 10⁰] are tested and the results are similar as well). We observe that the variation of parameters estimated by Model 5 with σ varying from 10⁻⁶ to 10⁰ does not exceed 1%. Therefore, the parameter estimation results are not sensitive to the choice of σ as long as σ is not too large.

3.2.5 Further Model Evaluation

The training and testing errors, , as defined in Appendix C (in the equation (A5)), are listed below in Table 4 with mean and standard deviation. We remind the readers that the superscripts ^[Tr] and ^[Te] represent the errors are computed on training and testing data respectively, and the subscript _(w) denotes that the error is weighted average of the daily relative errors over time. It can be seen from Table 4 that the presence of both heterogeneity and transportation helps reduce the training and testing errors by comparing the first four models. By comparing Model 4 and Model 5, it can be seen that using the graph Laplacian prior could further reduce the testing errors.

In addition, it can be seen from Table 4 that the errors increase greatly after transportation is included while heterogeneity remains absent. A possible explanation for this might be that without heterogeneity of parameters and transportation between provinces, the estimated λ_k’s are lower than the true λ_k’s for group 1, which leads to that the estimated newly confirmed cases are fewer than the true ones in group 1. For the same reason, the estimated newly confirmed cases are higher than the true ones in group 2. When the transportation is considered, more confirmed cases in group 1 are transferred to group 2 than the cases transported in the opposite direction. As a result, when the transmission parameters do not have heterogeneity, migration between provinces will worsen the prediction performance compared to the case without migration.

Furthermore, the last two rows of Table 4 report the training and testing errors for Model 5 with the same µ = 10^2.7 while σ= 10⁻³ and σ= 10⁰ respectively. As a consequence of the robustness of the parameter estimation regarding σ, the errors of Model 5 are also robust to σ. The similar analysis is also performed for the other data sets and the similar results can be obtained which we do not report repetitively. Hereinafter, the results are presented with σ= 10⁻⁶.

The plots of the mean of weighted and simply averaged testing errors and against varying µ are shown in Figure 4 and Figure A2 respectively. Recall that the subscripts (w) and (s) denote the weighted and simple average respectively. Note that Model 5 with µ = 10^2.7, which minimizes validation errors as shown in Figure 1 and Figure A3, achieves the minimal values of testing errors and (also and ). This is because validation and testing errors have the same distribution in this case. We also note that the testing error of Model 4 is nearly the same as that of Model 5 with small µ, since Model 4 is basically the same as Model 5 with µ = 0.

• Download figure
• Open in new tab

Figure 4:

Mean of weighted testing errors and against µ (computed by (A5)) over 100 replicas for simulated data (four provinces). The horizontal lines show the values of the mean of the testing errors when µ = 0.

3.3.1 Data Description

We set the total number of regions n = 30, the total days considered T = 30, and set the threshold that separates the training and testing set as T_th = 10. The other prefixed parameters are listed as below:

• For k ∈ {1, 2, …, n}, N_k(0) = 10⁶, E_k(0) = 30, H_k(0) = 10.

• For l ∈ {1, …, T }, i, j ∈ {1, 2, …, n} and i ≠j, (W_l)_ij = 5 × 10³.

• For provinces in Group 1:

  λ₄ = 0.307, λ₅ = 0.295, λ₆ = 0.303, λ₁₃ = 0.299, λ₁₆ = 0.296, λ₂₁ = 0.296, λ₂₅ = 0.294, λ₂₈ = 0.293, λ₂₉ = 0.305.

  For provinces in Group 2:

  λ₁ = 0.416, λ₂ = 0.410, λ₇ = 0.405, λ₈ = 0.393, λ₉ = 0.400, λ₁₁ = 0.402, λ₁₂ = 0.406, λ₁₅ = 0.407, λ₁₈ = 0.395, λ₁₉ = 0.410, λ₂₂ = 0.396, λ₂₃ = 0.419.

  For provinces in Group 3:

  λ₃ = 0.493, λ₁₀ = 0.512, λ₁₄ = 0.514, λ₁₇ = 0.503, λ₂₀ = 0.507, λ₂₄ = 0.498, λ₂₆ = 0.488, λ₂₇ = 0.495, λ₃₀ = 0.482.

Note that the thirty provinces are randomly assigned to three groups, and we denote this ground truth partition as P :

• D₁: Provinces 4, 5, 6, 13, 16, 21, 25, 28, 29.

• D₂: Provinces 1, 2, 7, 8, 9, 11, 12, 15, 18, 19, 22, 23.

• D₃: Provinces 3, 10, 14, 17, 20, 24, 26, 27, 30.

The transmission rates of the provinces in the same groups are closer as listed above. Specifically, the transmission rates of provinces in D₁, D₂, and D₃ are randomly sampled from normal distribution with mean 0.3, 0.4, and 0.5 respectively and all with standard deviation 0.01.

3.3.2 Models to Compare

We also compare the Models 1-5 listed in Table 2 as detailed in Section 3.2.2.

Note that for Model 5, the parameter σ in (2.8) is still fixed to be 10⁻⁶. When we assume that the ground truth group division P as listed above is known, in model performance evaluation as well as the optimization of µ for Model 5, by (2.11), are estimated by the optimization problem in (3.2): where β = 0.1.

Moreover, for Model 5, we also report results when there are mismatches in partition of the provinces, since in real-world applications the graph information may not be fully known. Specifically, since most of the results reported in the following sections are for Provinces 4 (in D₁), 1 (in D₂), and 3 (in D₃), we consider the other two partitions and , which are different from P in the assignments of Provinces 4, 1, and 3, and another three provinces respectively. Comparison between the results using the ground truth partition and the results of these two different kinds of mismatches may reflect the potentially different impact of these mismatches on the results. For Model 5 with partitions P′ or is estimated by (3.2) with P = {D₁, D₂, D₃} replaced by or .

The partition deviates from P in the assignment of Provinces 4 (in D₁), 1 (in D₂), and 3 (in D₃):

• : Provinces 3, 5, 6, 13, 16, 21, 25, 28, 29.

• : Provinces 4, 2, 7, 8, 9, 11, 12, 15, 18, 19, 22, 23.

• : Provinces 1, 10, 14, 17, 20, 24, 26, 27, 30.

The partition deviates from P in the assignment of Provinces 5 (in D₁), 2 (in D₂), and 10 (in D₃):

• : Provinces 4, 10, 6, 13, 16, 21, 25, 28, 29.

• : Provinces 1, 5, 7, 8, 9, 11, 12, 15, 18, 19, 22, 23.

• : Provinces 3, 2, 14, 17, 20, 24, 26, 27, 30.

For mismatched partitions P′ and P′′, the provinces assigned to the wrong groups are marked in bold.

We remark that the parameter inference for Models 1-4 is the same as detailed in Section 3.2.2.

3.3.3 Results of Trajectory Prediction

Note that the validation and testing errors have the same distribution for the simulated data. Therefore, for Model 5 (with the ground truth partition P), we choose µ = 10^1.8, which minimizes the mean of the validation errors as marked with red pentagrams in Figure 5 and Figure A3. The following results for Model 5 are all obtained with µ = 10^1.8.

• Download figure
• Open in new tab

Figure 5:

Mean of weighted validation errors and against µ (computed by (A5)) over 100 replicas for simulated data (thirty provinces). The horizontal lines show the values of the mean of the validation errors when µ = 0.

Then, we choose three provinces from the total thirty provinces (one province from each group), and plot the prediction of the trajectories and also the absolute errors from one specific replica as in Figures 6 and 7. It can be observed from Figure 7 that Model 5 (the green lines) achieve the most accurate prediction results due to the graph Laplacian regularization. Model 4 (the purple lines) also have good generalization performance, thanks to the heterogeneity and transportation involved in this model.

• Download figure
• Open in new tab

Figure 6:

True and fitted trajectories of Models 1-5 for three provinces chosen from the total thirty provinces. The vertical lines show the threshold T_th = 10.

• Download figure
• Open in new tab

Figure 7:

Absolute errors fitted trajectories of Models 1-5 for three provinces chosen from the total thirty ones. The vertical lines show the threshold T_th = 10.

We note that Model 3 also behaves much worse than the best models, especially in Provinces 4 (in D₁) and 3 (in D₃). In fact, the estimates for the λ_k’s in these three provinces using Model 3 are 0.2908, 0.4198, and 0.4639, respectively, which are not much different from the ground truth values. Nevertheless, Model 3 does not involve transportation between provinces. This leads to the confirmed cases of Province 3 in D₃ (which are much more than the cases of provinces in the other two groups) not output to other provinces; hence the predicted newly confirmed cases are much more than the truth. A similar situation happens in the case of Province 4 (in D₁). Namely, Model 3 does not take the imported cases from the provinces in the other two groups into account, causing the predicted newly confirmed cases to be much less than the truth.

3.3.4 Results of Parameter Estimation

The mean and standard deviation of λ_k’s in Provinces 4 (in D₁), 1 (in D₂), and 3 (in D₃) are reported in Table 5. We can still see that Models 3-5 estimate the transmission parameters more accurately than Models 1 and 2 especially for Provinces 4 and 3. In addition, the estimates by Model 5 of λ_k are closer than those by Model 4, due to the existence of the regularization term.

The last three rows show the parameter estimates using Model 5 with three partitions P, P′, and P′′ respectively. Recall that P is the ground truth underlying graph structure, and P′ and P′′ are mismatched partitions introduced in Section 3.3.2. It can be seen that compared with the estimated parameters using Model 5 with P′ that differs from P in the grouping of Provinces 4, 1, and 3, the estimates using Model 5 with P′′ that differs from P in the grouping of Provinces 5, 2, and 10 are closer to the estimates using Model 5 with the ground truth partition P. Hence, the results imply that for the parameter inference of simulated data generated from a certain partition, the incorrect division causes more discrepancy in the mismatched regions when comparing the estimates using the ground truth partition.

3.3.5 Further Model Evaluation

The mean and standard deviation of the weighted training and testing errors (formulas can be found in Appendix C, specifically in (A5)) are presented in Table 6. We see that Model 5 with µ = 10^1.8 achieves the minimum testing errors among all the models, and hence have the best generalization performance.

The last three rows show the comparison of errors for Model 5 with three partitions P, P′, and P′′ respectively. We may see that the testing errors of Model 5 with the incorrect partitions are basically the same and larger than those with the ground truth partitions. Nevertheless, these errors are still smaller than the other baseline models. Therefore, the incorrect division might worsen the generalization performance of Model 5; however, when the division does not deviate from the ground truth much, it still behaves better than Models 1-4 with uniform prior.

Moreover, as shown in Table 6, both Model 2, which allows transportation but not heterogeneity, and Model 3, that introduces heterogeneity but no transportation, behave worse than the baseline Model 1. Model 2 having more significant testing errors is due to the same reason as the case of four provinces, as has been analyzed in Section 3.2.5. The reason for Model 3 not being able to predict well is explained in detail in Section 3.3.3.

However, we also note that Model 4, which involves both transportation and heterogeneity of parameters, has greater testing errors than Model 1 on average. This phenomenon may be partly explained by the fact that, on the one hand, there may be more infected cases imported from regions with higher λ_k’s to those with lower transmission rates. On the other hand, although the transmission parameters are the same and all-around 0.42 in Model 1, E₀ and H₀ are still allowed to be spatially heterogeneous, which makes the trajectories not that far away from the ground truth. Meanwhile, though the λ_k’s are allowed to be spatially heterogeneous in Model 4, it often underestimates λ_k’s for provinces in D₁ and overestimates λ_k’s for provinces in D₃, making the trajectories deviate from the ground truth.

Figures 8 and A4 plot the change of weighted and simply averaged testing errors and regarding µ (recall that the subscripts (w) and (s) denote weighted and simple average respectively). We can still observe that Model 5 with the minimizers of validation errors also has the minimal testing errors, since the errors have the same distribution for the simulated data.

• Download figure
• Open in new tab

Figure 8:

Mean of weighted testing errors and against µ (computed by (A5)) over 100 replicas for simulated data (thirty provinces). The horizontal lines show the values of the mean of the testing errors when µ = 0.

Real-world Data in China

Data in China include the number of confirmed cases and removed cases (consisting of recovered cases and fatalities) in China’s major provinces and municipalities from January 21st to March 28th (2020). These publicly available data are downloaded from the websites of the National Health Commission of the People’s Republic of China [28] and Chinese Center for Disease Control and Prevention [29]. The corresponding total population in each province or municipality is from [41].

We also utilize two transportation data sets in China extracted from Baidu Qianxi [42], lasting from January 10th to February 10th in 2020. The first data set, called as migration indexes, reflects the number of people moving out from the provinces. The second data set reflects the percentage of people moving to the corresponding destinations from their starting points. The traveling volume is estimated by combining these two data sets.

Real-world Data in Europe

Data in Europe include the number of confirmed and removed cases (consisting of recovered cases and fatalities) in the following 11 countries from May 1st to August 31st in 2020: Denmark, Finland, Norway, Austria, Germany, Switzerland, Italy, Spain, Belgium, France, and Ireland. The data are downloaded from [43, 13], where the data are collected from European Center for Disease Prevention and Control [30]. The data of total population in each country is obtained from [44]. During the study, we are not able to obtain suitable transportation data between these European countries needed by our model, and thus we take simplification by incorporating only spatial and temporal heterogeneity and no transportation in this experiment, details in Section 4.4.

We remark that the data in China are from the initial outbreak of the epidemic. On the contrary, the data in Europe are from when the epidemic has been ongoing about six months. The choice of the aforementioned time period of interest is to illustrate that the proposed model can predict both the case when the epidemic breaks out from some seed region and the case when the epidemic has already established sizable local spread within different regions.

4.3.1 Data Description

After the selection process described in Appendix D, n = 21 provinces or municipalities are taken into consideration, which are Anhui, Beijing, Fujian, Gansu, Guangdong, Guangxi, Hebei, Henan, Hubei, Hunan, Jiangsu, Jiangxi, Liaoning, Ningxia, Shandong, Shan-Xi, Shanxi, Shanghai, Sichuan, Zhejiang, and Chongqing.

We set T = 31, from January 10th to February 10th 2020. However, we remark that only the data after January 21st are available. For measuring the fitting and generalization performance, we set the threshold T_th = 27 (February 5st). The data from T_th,1 = 11 (January 21st) to the T_th-th day (February 5st) are treated as the training set, and the data from T_th + 1 (February 6th) to the t-th day (February 10th) are treated as the testing set. Furthermore, for cross-validation of parameter µ of the proposed model, we choose another threshold T_th,2 = 17 (January 27th) and another ending data T₂ = 20 (January 30th). The parameters are first fitted on the data from T_th,1 (January 21st) to T_th,2 (January 27th), and then the validation errors are computed on the data from T_th,2 (January 27th) to T₂ (January 30th). µ is chosen according to validation errors.

In addition, from observation of data from January 10th to February 10th, the transmission rates in the selected provinces or municipalities change after some specific time point. Thus, we allow {λ_k} to be time-varying for this data set with appropriate changes to the model described in Section 2, which are detailed in Section 4.2. Recall that T_k denotes the day when λ_k changes in the k-th region. For real-world data in China, noting that multiple containment measures, including domestic travel restrictions, contact tracing and epidemiological investigation, and even city lock down were implemented since late January/ early February, we choose T_k = 24 (February 2nd) for Anhui, Beijing, Henan and Hubei, and T_k = 20 (January 29th) for the other provinces from the initial observation of the data.

Furthermore, the data from Baidu Qianxi might not be the exact traveling volumes between municipalities and provinces. We assume that the actual traveling volume from one starting point to one destination is proportional to Baidu migration index (which reflects numbers of people departed from the stating point) and the percentage of population traveling from this origin to the destination. The corresponding scaling parameter α also needs to be inferred from the data for all the models.

4.3.2 Models to Compare

As in Section 3.2.2, the same Models 1-5 are compared for COVID-19 data in China. Recall that Model 5 is the model proposed in this paper and the other four are the baseline models for comparison.

Note that for Model 5 with prior distribution based on graph Laplacian, the provinces are chosen to be divided into two groups, which are Hubei and other provinces except Hubei, due to the overwhelmed medical resources in Wuhan at the early stage of the pandemic [45]. For Model 5, by (4.3), the time-varying extension of (2.11) in Section 4.2, and scaling parameter α are estimated by the optimization problem in (4.4): where β = 0.1, and and are the transmission rates in the k-th province before and after the T_k-th day respectively, and σ is chosen to be 10⁻⁶. , where is the averaged traveling volume between the i-th and j-th province with rescaling, treated as the proximity between the two provinces.

For Models 1-4, if the model does not have heterogeneity of transmission parameters, then are forced to be the same and so do , while the transmission rates in the two periods are allowed to be different; if the model does not allow transportation between provinces, then terms involvingW_t disappear in (2.4) as described in Section 3.2.2.

4.3.3 Results of Trajectory Prediction

We first remark that for Model 5, three values of µ are chosen, and part of the results below are present from these three choices to perform the sensitivity analysis of µ. The first value is µ = 10^2.4, which minimizes the relative validation errors (, and ), marked with black pentagrams in Figure 9 below and Figure A5. The other two values are µ = 10² and µ = 10^2.7, whose validation errors are slightly larger, marked with red and green pentagrams in Figure 9 and Figure A5 respectively. Since for the real-world data in China, the validation error does not necessarily have the same distribution as the testing error, some of the results from various choices of µ are presented for better comparison.

• Download figure
• Open in new tab

Figure 9:

Weighted validation errors and against µ (computed by (A5)) for real-world data in China. The horizontal lines show the values of the validation errors when µ = 0.

Figure 10 shows the true and predicted trajectories for newly confirmed cases in Hubei. In Figure 10,

• Download figure
• Open in new tab

Figure 10:

True and fitted trajectories in Hubei. The red line shows the true trajectory, the black lines show the predicted deterministic trajectories with inferred using (4.4) with µ = 10^2.4, 10², 10^2.7 respectively), black and red scatter plots show 100 stochastic trajectories with Θ^∗ and sampled from the posterior distribution of respectively ( inferred using (4.4) with µ = 10^2.4). The vertical lines show the threshold T_th = 27 for the training data and testing data.

• The red line shows the true trajectory.

• The black lines show the predicted deterministic trajectories obtained by running (2.4) with Θ^∗ inferred using (4.4) with 10², 10^2.4, 10^2.7 respectively.

• The black scatter plots show 100 stochastic trajectories with Θ^∗ inferred using (4.4) with 10^2.4

• The red scatter plots show 100 stochastic trajectories sampled from the posterior distribution of , which is estimated by 5 × 10⁵ MCMC iterations staring from inferred using (4.4) with 10^2.4

The deterministic trajectories obtained by the other four models and the corresponding absolute errors are also plotted for better comparison. It can be seen that all the predicted deterministic trajectories generally capture the trend of the true trajectory. Moreover, sampling trajectories from the posterior distribution of generates more randomness than sampling trajectories with Θ^∗. Besides, it can be seen from Figure 10 that the performances of all the five models (including Model 5 with the three choices of µ) are similar for Hubei. This may be explained by the fact that the cases in Hubei outnumber those in other provinces or municipalities, which makes all the models tend to fit the trajectory of Hubei best.

Similarly, Figure 11 shows the predicted trajectories in Henan. As can be seen in Figure 11, unlike Hubei, heterogeneity of parameters helps the model fit and predict the trajectory of Henan in a more significant manner. When the transmission parameters in the two periods are forced to be the same in all provinces, the increment of the newly confirmed cases of the predicted trajectory is faster than the trend shown in the true trajectory in the first period. Meanwhile, the decrement of the newly confirmed cases of the fitted deterministic trajectory in the second period is also slower than the trend reflected in the true trajectory. Thus, heterogeneity helps improve the performance of fitting and predicting. Additionally, although µ = 10^2.4 minimizes the validation errors, it can be seen from Figure 11 that the deterministic trajectory obtained by Model 5 with µ = 10² achieves better performance in prediction than the one obtained by Model 5 with µ = 10^2.4. Therefore, the result suggests that the choice of µ is not necessarily limited to the minimizer of the validation error. Instead, we may also compare trajectories with µ that have slightly larger validation error for possibly better generalization performance.

• Download figure
• Open in new tab

Figure 11:

True and fitted trajectories in Henan. The red line shows the true trajectory, the black lines show the predicted deterministic trajectories with inferred using (4.4) with µ = 10^2.4, 10², 10^2.7 respectively), black and red scatter plots show 100 stochastic trajectories with Θ^∗ and sampled from the posterior distribution of respectively ( inferred using (4.4) with µ = 10^2.4). The vertical lines show the threshold T_th = 27 for the training data and testing data.

Moreover, from Figure 12 which shows the fitted trajectories in Anhui, it can still be seen that heterogeneity helps the models predict the testing trajectory better. Furthermore, introducing correlation helps Model 5 (with µ = 10^2.4) predict the trajectory marginally better than the other models.

• Download figure
• Open in new tab

Figure 12:

True and fitted trajectories in Anhui. The red line shows the true trajectory, the black lines show the predicted deterministic trajectories with inferred using (4.4) with µ = 10^2.4, 10², 10^2.7 respectively), black and red scatter plots show 100 stochastic trajectories with Θ^∗ and sampled from the posterior distribution of respectively ( inferred using (4.4) with µ = 10^2.4). The vertical lines show the threshold T_th = 27 for the training data and testing data.

4.3.4 Results of Parameter Estimation

To illustrate the estimate of the posterior of parameters, Figure 13 shows the estimated posterior distributions of and in Hubei using Model 5, in which the red lines show the inferred and with µ₁ = 10^−3.5, and Table 7 shows the inferred and with the mean and standard deviation of their estimated posterior distributions from MCMC.

• Download figure
• Open in new tab

Figure 13:

Estimated posterior distribution of λ in Hubei with µ₁ = 10^2.4. The vertical red lines represent the values of corresponding θ^∗.

By comparing the transmission rates in Hubei in the two periods, it can be seen that the transmission rate decreases greatly after February 2nd, which indicates that the containment measures adopted in Hubei against the COVID-19 outbreak are effective.

4.3.5 Further Model Evaluation

Table 8 shows the training errors and testing errors of the five models computed as detailed in Appendix C. As can be seen from Table 8, both heterogeneity of parameters and transportation between provinces help reduce the training errors and the testing errors as the simulated data case. The utilization of graph Laplacian prior helps further reduce the testing errors. Note that the testing errors of Model 4 are close to those of Model 3, which might be explained as after the traveling restrictions are imposed on January 23rd [46], transportation no longer has much influence on the spread of the epidemic.

Figures 14 and A6 plot weighted and simply averaged testing errors against varying µ respectively. We can see from both Table 8 and Figure 14 that though µ = 10^2.4 achieves the minimal validation error, it has slightly larger testing error than µ = 10².

• Download figure
• Open in new tab

Figure 14:

Weighted testing errors and against µ (computed by (A5)) for real-world data in China. The horizontal lines show the values of the testing errors when µ = 0.

4.4.1 Data Description

11 countries in Europe are considered in the experiments, which are Denmark, Finland, Norway, Austria, Germany, Switzerland, Italy, Spain, Belgium, France, and Ireland. Some countries such as Iceland and Luxembourg are excluded since their total population is too small, and some other countries like Sweden, Holland, and England are ignored due to that the removed cases in these countries are zero (from the data sources [43]) and thus the two-step procedure for parameter estimation in this paper may not be applied.

Since the original data oscillate much, especially over the weekends, the data are first smoothed by taking a simple moving average (which treats the unweighted average of data in 7 days around the i-th day as the new data point for the i-th day). After the smoothing, T = 116.

Additionally, from the observation of newly confirmed cases from May 1st to August 1st, the transmission rates change after some threshold in the selected countries in Europe. Thus, we still allow {λ_k} to be time-varying for this data set as described in Section 4.2. T_k = 50 is chosen to be the day when λ_k changes for all countries except Norway, Italy and France. In addition, T_Nowway = 75, T_Italy = 85 and T_{F rance} = 60. Furthermore, for evaluating fitting and generalization performance, we set T_th = 95. For choosing µ (the details can be found in Appendix B.2), we set T_th,2 = 18 and T₂ = 36. Namely, the validation errors are computed by comparing the data from T_th,2-th day to T₂-th day with the predicted trajectories using parameters fitted on the data from the 1st day to T_th,2-th day. Then, the choice of µ is based on validation errors.

4.4.2 Models to Compare

For COVID-19 data in Europe, only three models detailed below are compared, which allow heterogeneity of parameters but do not include transportation. The last model is the proposed one in this paper, and the first two are baseline models. The transportation data needed in this model are not available to the best of our knowledge. However, the impact of transportation is expected to be less significant due to travel restrictions [46].

• 1’. Model with uniform prior distribution, without heterogeneity or migration.

• 2’. Model with uniform prior distribution, with heterogeneity but without migration.

• 3’. Model with prior distribution based on graph Laplacian, with heterogeneity but no migration.

  For Model 3’, 11 countries are partitioned into d = 4 groups according to geographical locations:

  • D₁ (Northern Europe): Denmark, Finland, Norway;

  • D₂ (Central Europe): Austria, Germany, Switzerland;

  • D₃ (Southern Europe): Italy, Spain;

  • D₄ (Western Europe): Belgium, France, Ireland.

We denote this partition as P. Since this might not be the unique appropriate partition for the n = 11 countries, we also present results from another partition denoted as P′ in the following sections. P′ also groups the countries that are geographically close together:

• : Finland, Norway;

• : Denmark, Austria, Germany;

• : Spain, Belgium, France, Ireland;

• : Switzerland, Italy.

The results for Model 3’ in Section 4.4.5 are presented with both the partitions P and P′.

Besides, by the same time-varying extension (4.3) detailed in Section 4.2, for Model 3’, are estimated by the optimization problem in (4.5): or by (4.5) with the partition P = {D₁, D₂, D₃, D₄} replaced by another partition . Here, σ is still taken to be 10⁻⁶, and β is still 0.1.

4.4.3 Results of Trajectory Prediction

We still remark that for Model 3’ (with P), as in the case of real-world data in China, part of the results reported are from three choices of µ for careful consideration, since the validation error may not have the same distribution as the testing error. Note that Figures 15 and A7 plot the weighted and simply averaged relative validation errors, and respectively. One choice is µ = 10^3.9, which minimizes the validation errors (computed with P), marked with black pentagrams in Figures 15 and A7. The other two choices are µ = 10^3.2 and µ = 10⁵, which correspond to slightly larger validation errors and are marked in red and green pentagrams respectively in Figures 15 and A7.

• Download figure
• Open in new tab

Figure 15:

Weighted validation errors and against µ (computed by (A5)) for real-world data in Europe. The horizontal lines show the values of the validation errors when µ = 0.

Figure 16, 17 and 18 present the true and predicted trajectories in Austria, Germany, and Italy, respectively. As can be seen from Figure 16, 17 and 18, heterogeneity of transmission parameters in Model 2’ (green lines) and Model 3’ (black lines) helps improve the performance of fitting and generalization. Furthermore, utilization of correlation between countries in Model 3’ further improves the prediction of the trajectories. We also note that for Model 3’, compared to the results for µ = 10^3.9, which minimizes the validation errors, µ = 10⁵ predicts the trajectories of Austria and Italy more accurately. Hence, the results for Model 3’ are sensitive to µ, and it might be better to compare the trajectories from multiple choices of µ.

• Download figure
• Open in new tab

Figure 16:

True and fitted trajectories in Austria. The red line shows the true trajectory, the black lines show the predicted deterministic trajectories with inferred using (4.4) with µ = 10^3.9, 10^3.2, 10⁵ respectively), black and red scatter plots show 100 stochastic trajectories with Θ^∗ and sampled from the posterior distribution of respectively ( inferred using (4.4) with µ = 10^3.9). The vertical lines show the threshold T_th = 102 for the training data and testing data.

• Download figure
• Open in new tab

Figure 17:

True and fitted trajectories in Germany. The red line shows the true trajectory, the black lines show the predicted deterministic trajectories with inferred using (4.4) with µ = 10^3.9, 10^3.2, 10⁵ respectively), black and red scatter plots show 100 stochastic trajectories with Θ^∗ and sampled from the posterior distribution of respectively ( inferred using (4.4) with µ = 10^3.9). The vertical lines show the threshold T_th = 102 for the training data and testing data.

• Download figure
• Open in new tab

Figure 18:

True and fitted trajectories in Italy. The red line shows the true trajectory, the black lines show the predicted deterministic trajectories with inferred using (4.4) with µ = 10^3.9, 10^3.2, 10⁵ respectively), black and red scatter plots show 100 stochastic trajectories with Θ^∗ and sampled from the posterior distribution of respectively ( inferred using (4.4) with µ = 10^3.9). The vertical lines show the threshold T_th = 102 for the training data and testing data.

4.4.4 Results of Parameter Estimation

Same as before, Figure 19 shows the estimated posterior distributions of and using Model 3’, in which the red lines show the inferred and with µ = 10^3.9. Table 9 shows the inferred and and the mean and standard deviation of their estimated posterior distribution from MCMC. The estimated is larger than the estimated , which is consistent with the trend that the newly confirmed cases in Italy first decrease and then increase.

• Download figure
• Open in new tab

Figure 19:

Estimated posterior distributions of λ in Italy with µ = 10^3.9. The vertical red lines represent the values of corresponding θ^∗.

4.4.5 Further Model Evaluation

Table 10 shows training and testing errors of the three models. It can be seen from Table 10 that Models 2’ and 3’ have smaller both training and testing errors than Model 1’, due to the introduction of heterogeneity of transmission parameters. We can also see that compared with Model 2’, Model 3’ that adds Graph Laplacian prior has better generalization performance. Especially, Model 3’ with µ = 10⁵ greatly reduces the testing errors, although it does not have the minimal validation errors.

Recall that P and P′ are two different partitions of countries in Europe introduced in Section 4.4.2, both based on geographical locations. By comparing the corresponding results in Table 10 for Model 3’ with partitions P and P′ respectively, it can be seen that the testing errors are almost the same. This is different from the implications of the simulated data. A possible explanation for this might be that the regions in the real world have more complicated correlations than the simulated regions where clustering information is uniquely and artificially prefixed. Although reasonable groupings may not always be unique for real-world cases, the proposed model could still predict the trajectories more accurately than the baseline models.

Figure 20 below and Figure A8 plot the weighted and simply averaged testing errors (whose definitions are in Appendix C) against varying µ respectively, from which one may see that by imposing regularization properly through choosing a moderate µ helps improve the generalization performance. By comparing

• Download figure
• Open in new tab

Figure 20:

Weighted testing errors and against µ (computed by (A5)) for real-world data in Europe. The horizontal lines show the values of the testing errors when µ = 0.

Figure 20 (and Figure A8) and Figure 15 (and Figure A7), we see that the change in validation and testing errors concerning increasing µ are not the same. This further verifies the indications from the previous findings that it would be better to examine and compare the results from multiple choices of µ with relatively low validation errors for Model 3’ to choose a better µ for the specific data sets.