Non-stationary Spatio-Temporal Modeling of COVID-19 Progression in The U.S.

2.1. Model 1: Piecewise Stationary SIR Model

The assumptions underlying the transmission and recovery rates of the standard SIR model are stringent. Both environmental factors and changes in population behavior can lead to time varying behavior and this has been the case for Covid-19; see, e.g., discussion in (Giordano et al., 2020). Variants of the SIR model with time varying parameters have been proposed in the literature (Liu and Stechlinski, 2012).

Let β(t) and γ(t) denote the transmission and recovery rate at time t. A time varying analogue of the system in (3) is:

The solutions for β(t) and (t) are given by

Given time series data on {I(t),R(t)),t ∈ [1,T]}, the corresponding B_t =(β(t),γ(t))‣ can be calculated for each t =1, ···, T − 1 (Chen, Lu and Chang, 2020).

For our application, we assume that the transmission and recovery rates are piecewise constant over time, reflecting the fact that their temporal evolution is impacted by intervention strategies and environmental factors. Hence, it becomes important to identify the break points wherein these rates exhibit changes. Specifically, suppose there exist m₀ change points such that the rate vector B_t =(β(t),γ(t))′ exhibits a change from B^(j)=(β^(j),γ^(j))′ to B^(j+1)=(β^(j+1),γ^(j+1))‣ at time point t_j, while it remains fixed until the next break point. Then, model (4) takes the following form where . In each segment [t_j−1,t_j), all model parameters are assumed to be fixed; hence, these break points divide the time series data into stationary segments, as reflected in model (6).

The goal then becomes to utilize data on infected and recovered individuals over time to detect both the change points t_js, as well as the model parameters β_j and γ_j within each fixed segment. The key steps in our proposed change point detection strategy are outlined in Algorithm 1.

Next, we elaborate on the key steps of the Algorithm. First, few notations are defined.

Notation. Denote the indicator function of a subset A as l_A. R denotes the set of real numbers. For any vector v ∈ ℝ_p, we use ‖υ‖₁, ‖υ‖₂ and ‖υ‖_∞ to denote and max_1≤i≤p |υ_i|, respectively. The transpose of a matrix A is denoted by A′.

Block Fused Lasso: the objective is to partition the observations into blocks of size b_n, wherein the model parameter B_t remains fixed within each block and select only those blocks for which the corresponding change in the parameter vector is much larger than the others. Specifically, let n = T − 1 be the number of the times points for the response data Y_t and define a sequence of time points 1 = r₀ <r₁ < ··· <r_kn = n + 1 for block segmentation, such that r_i − r_i−1 = b_n for i =1,…,k_n − 1, , where is the total number of blocks. For ease of presentation, it is further assumed that n is divisible by b_n such that r_i − r_i−1 = b_n for all i =1,…,k_n. By partitioning the observations into blocks of size b_n and fixing the model parameters within each block, we set θ₁ = B⁽¹⁾ and for i =2, 3,…,k_n. Note that θ_i ≠ 0 for i ≥ 2 implies that θ_i has at least one non-zero entry and hence a change in the parameters. Next, we formulate the following linear regression model in terms of : where , , and E ∈ ℝ²ⁿ.

A simple estimate of parameters Θ(k_n) can be obtained by using an ℓ₁-penalized least squares regression of the form which uses a fused lasso penalty to control the number of change points in the model. This penalty term encourages the parameters across consecutive time blocks to be similar or even identical; hence, only large changes are registered, thus aiding in identifying the change points. Further, a hard-thresholding procedure is added to cluster the jumps into two sets: large and small ones, so that those redundant change points with small changes in the estimated parameters can be removed. We only declare that there is a change point at the end point of a block, when associated with large jump of the model parameters.

Hard Thresholding: is based on a data-driven procedure for selecting the threshold η. The idea is to combine the K-means clustering method (Hartigan and Wong, 1979) with the BIC criterion (Schwarz et al., 1978) to cluster the changes in the parameter matrix into two subgroups. The main steps are:

• Step 1 (initial state): Denote the jumps for each block by and let v₁ = 0. Denote the set of selected blocks with large jumps as J (initially, this is an empty set) and set BIC^old = ∞

• Step 2 (recursion state): Apply K-means clustering to the jump vector with two centers. Denote the sub-vector with a smaller center as the small subgroup, V_S, and the other sub-vector as the large subgroup, V_L. Add the corresponding blocks in the large subgroup into J. Compute the BIC by using the estimated parameters after setting for each block i ∉ J and denote it by BIC^new. Compute the difference BIC^diff = BIC^new − BIC^old and update BIC^old = BIC^new. Repeat this step until BIC^diff ≥ 0.

Block clustering: the Gap statistic (Tibshirani, Walther and Hastie, 2001) is applied to determine the number of clusters of the candidate change points. The basic idea is to run a clustering method (here, K-means is selected) over a grid of possible number of clusters, and to pick the optimal one by comparing the changes in within-cluster dispersion with that expected under an appropriate reference null distribution (for more details, see Section 3 in Tibshirani, Walther and Hastie (2001)).

Exhaustive search: Define and , where C_i’s are the subsets of candidates blocks by block clustering procedure. Denote the subset of corresponding block indices by J_i. Define the following local coefficient parameter estimates: for all i =1,…,m₀, where J₀ = {1} and .

Now, given a subset C_i, we apply the exhaustive search method for each time point s in the interval (l_i,u_i) to the data set truncated by the two end points in time, min(C_i) − b_n and max(C_i)+ b_n, i.e. only consider the data within the interval [min(C_i)−b_n, max(C_i)+b_n). Specifically, define the final estimated change point as for i =1,…,m₀, where and are the local coefficient parameter estimates based on the first step by block fused lasso. Denote the set of final estimated change points from (10) by .

Remark: An alternative approach for detection of break points is to run a full exhaustive search procedure for both single and multiple change point problems. Such procedures are computationally expensive, and not scalable for large data sets. Simple fused lasso is (block of size one) another method, which although computationally fast, it leads to over-estimating the number of break points; hence, it requires additional “screening” steps to remove redundant break points found using the fused lasso algorithm (Harchaoui and Lévy-Leduc, 2010; Safikhani and Shojaie, 2020). Such screening steps usually include tuning several hyper-parameters. This task not only slows down the detection method, but is also not robust. The proposed approach (block fused lasso coupled with hard-thresholding) aims to solve this issue by choosing appropriate block sizes, while it only needs a single tuning parameter (the threshold) to be estimated. Estimation of Infection and Recovery rates: Once the locations of break points are obtained, one can estimate the model parameters by running a separate regression for each identified stationary segment of the time series data. The work of Safikhani and Shojaie (2020) shows that this strategy yields consistent model parameter estimates.

Some Properties of the Detection Procedure: The estimated change points obtained by the block fused lasso step includes all points with non-zero , which leads to overestimating the number of the true change points in the model. Hence, the reason for employing a hard-thresholding step to “thin out” redundant change points exhibiting small changes in the estimated coefficients. Nevertheless, the block fused lasso parameter estimates enjoy a prediction consistency property, which implies that the prediction error, i.e. , converges to zero with high probability as n → +∞. This result is stated and proved in Theorem A.1 in the Appendix A under some mild conditions on the behaviour of the tail distributions of error terms.

After the hard-thresholding step, those candidate change points located far from any true change points will be eliminated when the block size is large enough (b_n = c log n for a large enough constant c> 0). On the other hand, there may be more than one selected change points remaining in the set in b_n-neighborhoods of each true change point. The exhaustive search step picks one estimated break point in each such b_n-neighborhoods. Therefore, by selecting b_n = c log n for a large enough constant c> 0, one can conclude that the proposed detection algorithm locates the true break points with an error bounded by the order . Similar bounds have been found in the literature for similar processes (see for example Chan, Yau and Zhang (2014) for auto-regressive models). For more details, see appendix A.

2.2. Model 2: Piecewise Stationary SIR Model with Spatial Heterogeneity

The standard SIR model and its piecewise stationary counterpart do not account for any influence due to inter-region mobility and travel activity. To that end, we extend the previous model for each region, by considering the influence exerted by its q neighboring regions (cities, counties or states) (see also Srivastava and Prasanna (2020)). Then, for each time point t =1, ···, T − 1, the model becomes where N^j denotes the total population in neighboring region j; α is a spatial effect parameter; ω_j’s are spatial weights such that . The latter two parameters capture inter-region mobility patterns.

Model (11) can be succinctly written as where and is the weighted average of spatial effect over the neighboring regions at time t.

Note that assuming that the influence of the spatial effect component Z is small, we can use Model 1 for each region of interest to estimate both the change points and the corresponding transmission and recovery rates. Subsequently, the piecewise constant SIR model, augmented with an overall spatial effect, is used as follows: suppose that there are m₀ change points being detected, denote the final estimated change points by Model 1 as . We now formulate the following linear regression model in terms of : where , , and E ∈ ℝ²ⁿ. The model parameters are then estimated by least squares.

2.3. Model 3: A Stochastic Piecewise Stationary SIR Model with Spatial Heterogeneity

The standard homogeneous SIR model, previously discussed, is deterministic; hence, the output of the model is fully determined by the parameter values of the transmission and recovery rates and the initial conditions. Its stochastic counterpart (Bartlett, 1949; Bailey, 1953), possesses some inherent randomness. The same set of parameter values and initial conditions will lead to an ensemble of different outputs; for an in depth discussion and survey of deterministic and stochastic epidemic models, see Daley and Gani (2001).

The transition probabilities associated with the Markov process underlying the stochastic version of the SIR model are given by where S(t),I(t) ∈{0, 1, 2, ···, N} are discrete random variables, denoting the number of susceptible and infected individuals at time t.

Using the approximation result in Greenwood and Gordillo (2009) (see also Allen (2017); Buckingham-Jeffery, Isham and House (2018)), the general stochastic epidemic model can be approximated by the stochastic differential equation: where the random variables S(t) and I(t) are continuous, and W =(W₁,W₂)′ is a vector of two independent Wiener processes, i.e., .

Under the assumption that S(t) ≈ N and given (16), the stochastic SIR model can be written as: where the increments ΔW₁(t) and ΔW₂(t) are two independent normal random variables, i.e., .

It can be seen that the resulting regression model, based on the discrete analogue of (18), will have an error term exhibiting temporal correlation, driven by I(t) and R(t). An examination of the temporal correlation patterns in the residuals of model 2 from COVID-19 data (see left two panels in Figures 8 and 11 in Appendices B and C) supports this finding. To that end, we model the error process as a Vector Auto-Regressive (VAR) with lag p. The corresponding temporal correlation plots for the residuals after inclusion of the VAR(p) component are depicted in the right two panels in Figures 8 and 11 in Appendices B and C and clearly show the importance of considering such an error structure.

The piecewise stationary SIR model with spatial heterogeneity and a VAR(p) error process is given by where u_t is white noise with mean 0 and variance σ², ϕ₁;…;ϕ_p are the corresponding autoregressive parameters.

The estimation of the model parameters is accomplished in the following three steps: Step 1: Fit model 1 for each region of interest to obtain the change points; Step 2: Obtain the transmission and recovery rates and spatial effect as in Model 2. Step 3: Compute the residuals from Step 2 and fit a VAR model to them, see Lütkepohl (2005).

3.1. Simulation Scenarios

We consider three different simulation settings. In each simulation scenario, we consider three different block size settings. In all scenarios, the block sizes are selected to be b_n =4, 8, 12.

Simulation Scenario A (SIR model with piecewise constant rates): the transmission and recovery rates are chosen to be piecewise constant. In this scenario, we set the number of time points T = 250, m₀ = 2, the change point t₁ = 100 and t₂ = 200. We choose β⁽¹⁾ =0.10, β⁽²⁾ =0.05, β⁽³⁾ =0.01,γ⁽¹⁾ =0.04, γ⁽²⁾ = 0.06, γ⁽³⁾ =0.04. Results are based on data generated from the SIR model in (4) with β(t) ∼ and . The coefficients for the SIR model are depicted in Figure 1.

• Download figure
• Open in new tab

Fig 1:

True transmission rate and recovery rate in Scenario A, B, C. The left plot provides the true transmission rate β(t) and recovery rate γ(t) in the SIR model in Scenario A; the middle plot provides the true transmission rate β(t) and recovery rate γ(t) in the SIR model in Scenario B and C; the right plot provides the true transmission rate β^s(t) and recovery rate γ^s(t) for generating the spatial component in Scenario B and C.

Simulation Scenario B (SIR model with piecewise constant rates and spatial effect): the transmission and recovery rates are chosen to be piecewise constant. In this scenario, we set the number of time points T = 200, m₀ = 1, the change point . We choose α = 1, β⁽¹⁾ =0.10, β⁽²⁾ =0.05, γ⁽¹⁾ =0.04, γ⁽²⁾ =0.04, , , , . We construct a sequence of T − 1 values of decreasing from to and a sequence of T − 1 values of increasing from to on the linear scale and then generate the spatial effect data from SIR model in (4). By plugging in the spatial effect data, we generate the response variable Y_t from (11) with additional white noise error term from . The SIR model’s coefficients for the response variable and the spatial effect variable are depicted in Figure 1. Simulation Scenario C (piecewise constant SIR model with spatial effect and VAR(1) structure in error term): the data are generated based on (19) with piecewise constant transmission and recovery rates. Similar to scenario B, we set the number of time points T = 200, m₀ = 1, the change point , β⁽¹⁾ =0.10, β⁽²⁾ = 0.05, γ⁽¹⁾ =0.04 and γ⁽²⁾ =0.04. For the spatial effect, we set α = 1, , , , . The sequence of T − 1 values of and are on the linear scale, same as above in scenario B. We first generate the spatial effect data from SIR model in (4) with parameter and and generate the error term by VAR(1) model with the covariance matrix of the noise process Σε =0.1I₂. By plugging in the spatial effect data and error term data, we generate the dataset of the response variable Y_t from (19). The SIR model’s coefficients for the response variable and the spatial effect variable are depicted in Figure 1. The autoregressive coefficient matrix has entries 0.8, 0, 0.2, 0.7 from top left to bottom right.

3.2. Simulation Results

The mean and standard deviation of the location of the selected change point, relative to the the number of time points T – i.e., – for all simulation scenarios are summarized in Table 1. The results clearly indicate that, in the piecewise constant setting, our procedure accurately detects the location of change points. The results of the estimated transmission rate , recovery rate and spatial effect in Table 2 also suggest that our procedure produces accurate estimates of the parameters, under the various settings considered.

In Scenario C, we generate additional 20 days worth of data for prediction testing. The results of the out-of-sample mean relative prediction error (MRPE) for I(t) and R(t) are provided in Table 3. The results indicate that adding the spatial effect can significantly improve the prediction, when the spatial component influences the individual data series.

The estimated transition matrices of the VAR component are close to the true values, regardless of block size selection. For example, for b_n = 4, the entries of estimated autoregressive coefficients matrix which are averaged among 100 replicates are 0.7577, -0.0057, 0.2361, 0.6663, with standard deviations 0.1125, 0.0461, 0.3582, 0.1147, respectively, while the entries of true autoregressive coefficients matrices are 0.8, 0, 0.2, 0.7 from top left to bottom right. Moreover, for b_n = 8, the entries of estimated autoregressive coefficients matrix are 0.7652 -0.0057, 0.2088, 0.6767, with standard deviations 0.0533, 0.0464, 0.0459, 0.0573, respectively. These results confirm the good performance of the proposed algorithm in finite sample while also indicates its robustness with respect to block size selection.

4.1. Data Description

The COVID-19 data used in this study are obtained from The New York Times (2020). The curated data and code used in the analysis are available at the authors’ GitHub repository¹. The analysis is performed both at the state and county level and the raw data include both cases and deaths, as reported by state and local health department and organized by the NY Times. However, due to lack of complete information on recovered individuals (which are an important covariate in the models considered but the daily number of recovered cases are only reported at the national level (Wang et al., 2020b)), we calculate the number of recovered cases for each region (state/county) as follows: death in the region, multiplied by the nationwide cumulative recovered cases and divided by the nationwide deaths. Specifically, here we assume that the recovery versus deceased ratio for each state/county is fixed, and is approximated by the nationwide recovery to-death ratio. As coronavirus infections keep multiplied and laboratory testing is limited, reporting only confirmed cases and deaths leads to inconsistent data collection and undercounts of the disease’s impact. On April 14, 2020, CDC advised states to count both confirmed and probable cases and deaths. As more states and localities have started to include probable cases and deaths in their counts, in this study, we rely on the combined cases, which include both confirmed and probable cases. The populations of states and counties are obtained from the U.S. Census Bureau (2019). Further, to decide which neighboring states/counties to include in Model 2, their distance to the target state/county of interest is used. The latter is obtained from National Bureau of Economic Research (2010). Here, we define the regions within 500 miles for states and 100 miles for counties/cities as neighboring regions in model 2 and 3. For those areas with a large number of neighboring regions, such as New York state, New York City and King county, we only consider the top five regions with the smallest distances.

We analyzed the daily count of COVID-19 cases at the state-level from March 1, 2020, to August 18, 2020, and at the county-level from the first day after March 1, 2020, when the region records at least one positive COVID-19 case to August 18, 2020. In particular, for the five states presented next (NY, OR, FL, CA, TX), New York city and King county, n = 170; for Miami-Dade county, n = 160; for Orleans parish, n = 161; for Jefferson parish, n = 162; for Charleston, Greenville, Richland and Horry counties, n = 165, 156, 154, 156, respectively.

Most of the states together with New York City, King and Miami-Dade counties were selected due to being severely affected for a period of time during the course of COVID-19. The remaining counties illustrate interesting patterns gleaned from the proposed models.

Let I(t) and R(t) denote the number of infected and recovered individuals (cases) on day t. Day 1 refers to the first day after March for which the region records at least one positive COVID-19 case. Figure 2 depicts the actual case numbers I(t) and R(t) in the five states considered.

• Download figure
• Open in new tab

Fig 2:

Number of infected cases (dark orange) and recovered cases (dark green) in several states. The first vertical black dotted line indicates the “stay-at-home” order start date for each state, while the second vertical black dotted line indicates the “reopening” start date.

4.2. Results for US States

The various models considered are applied on scaled versions (divided by their standard deviations) of the predictors matrix X_t and the response vector Y_t. Further, the block size is set to b_n = 7, i.e. partition the observations into blocks of 7 days for all regions except Texas, wherein a block size b_n = 5 is used. We choose different block sizes since, in some regions, there are more than one break points of the transmission rate very close to each other which violates the assumptions of our model. To mitigate this issue, we select a smaller block size to detect most of the change points. Another situation is that the transmission rate in some regions exhibit many small noise disturbances. In this case, we use a larger block size to smooth out the parameter estimates and detect the change point more accurately. Finally, the tuning parameter λ_n is selected via a cross-validated grid search (Hastie, Tibshirani and Friedman, 2009).

For Model 2, we consider four different types of weights: equal weights (Model 2.1), distance-based weights (Model 2.2), similarity-based weights (Models 2.3 and 2.4). In Model 2.1 and 2.2, the neighboring regions are selected based on distance. For states, a threshold of 500 miles is used and the resulting neighbors are displayed in Table 4. When the spatial resolution is high (county level aggregated data), neighboring regions may exhibit similar patterns in terms of the evolution of transmission and recovery rates. Therefore, constructing weight matrices based on distance is meaningful as it is a common practice in spatial statistics (Cressie, 2015). Such spatial smoothing through proper weight matrices is specially helpful in increasing the statistical power through increasing the sample sizes, hence yielding to more accurate predictions.

However, the evolution of COVID-19 may exhibit different patterns across neighboring states due to the coarse spatial resolution. Thus, defining weight matrices based on distance may not be ideal. Hence, Models 2.3 and 2.4 select regions based on similarities of infected/recovered cases. Similarity between the region of interest and the j-th potential similar region is defined as

Model 2.3 uses the top five regions with the smallest similarity score, while Model 2.4 uses all states in the country. For states, the resulting neighbors for Model 2.3 are displayed in Table 5.

In the equal weight setting (Model 2.1), ω_j =1/q for any j =1,…,q. In both the distance-based weight and similarity-based weight settings, power distance weights are used, wherein the weight is a function of the distance/similarity to the neighboring region: where d_j is the distance score and s_j is the similarity score for the j-th region. Under the constraint that , we obtain the normalized weights as:

Before applying Model 3, we compare the in-sample and out-of-sample MRPEs (defined in Section 3) in all four variants of Model 2 and select the parameter values estimated by the best-performing model. In subsequent analysis, the results from Model 2.3 are reported, since it proved to be the best performing one.

As expected, change points detected for state data are related to “stay-at-home” orders, or phased reopening dates issued by state governments. We define the reopening date as the time when either the “stay-at-home” order expired or state governments explicitly lifted orders and allowed (select or even all) businesses to reopen (Mervosh et al., 2020). The “stay-at-home” and reopening dates for each state are shown in Table 6. In Model 1, a change point is detected from March to April for four states: New York, Oregon, Florida, and California. These change points coincide with the onset of “stay-at-home” orders and correspond to a significant decrease in the transmission rate. The first change point detected for Florida is around two weeks after the state’s Governor signed a statewide “stay-at-home” order, which is consistent with the fact that COVID-19 symptoms develop 2 days to 2 weeks following exposure to the virus. The first change points in New York and Oregon are detected around 20 days after the statewide “stay-at-home” order, which is longer than Florida. This delayed change point is probably because New York was a hot spot at that time since the virus was in circulation before the start of our observation period (March 1), and thus it took a much longer time to suppress the transmission rate by the statewide “stay-at-home” order. On the other hand, the death toll in New York, which contributes to the recovery rate, continued to rise from the end of March to the beginning of April, even though the “stay-at-home” policy was in effect. Although Texas exhibits a clear downward trend of transmission in April, as gleaned from the estimated transmission rate β(t) shown in the Figure 7, Model 1 fails to detect the change point for Texas in April. This is mainly because its transmission rate keeps decreasing smoothly during the month of May and violates the model’s piecewise constant assumptions.

Remark: The daily new cases data in California and Texas are highly skewed. This could be another explanation on why the detected change points by our algorithm are not close enough to the “stay-at-home” order start date. To remedy this issue, we considered applying log transformation to the data for these two states to deal with the skewness. The new break points detected by our method applied to the transformed data are much closer to the “stay-at-home” order start date. Specifically, for California, the detected change point occurs at April 8th -almost two weeks after the “stay-at-home” order start date-while the detected change point for Texas occurs at April 10th (almost a week after the “stay-at-home” order start date). A more detailed summary of these findings are stated in Section 5 of the supplementary materials.

As can be seen from Figure 7, in addition to the downward trend after lockdowns have been put in place, Florida and Texas have a clear upward trend after state reopenings began in May. Two change points in June and July are detected for Texas and Florida, which relate to their reopening. Note that in Florida, there is another jump in estimated parameters in June 21st. However, since this jump is close to the already selected change point in July, the algorithm does not pick this time point as a new break point since it violates the modeling assumption that true break points can not be close to each other. Such assumptions are quite natural and common in the anomaly detection literature (Safikhani and Shojaie, 2020).

Note that the restriction in either phase 2 reopening plan in Florida or the phase 3 reopening plan in Texas is quite similar in terms of restaurants, bars, and entertainment businesses. Starting June 5, restaurants and bars in Florida could increase their indoor seating to 50% capacity. Movie theaters, concert venues, arcades, and other entertainment businesses can also open at 50% capacity. Starting June 3, all businesses in Texas could expand their occupancy to 50% with certain exceptions. Bars may increase their capacity to 50% as long as patrons are seated.

Assuming that Florida had not begun the phase 2 reopening plan on June 5, our model predicts 50,585 infected cases by June 12 while the actual number of infected cases is 52,269 (3.3% higher). Similarly, by June 19, our model predicts 57,895 infected cases, while the actual number of infected cases is 69,571 (20.2% higher). Similarly, suppose that Texas had not begun the phase 3 reopening plan on June 3, then by June 10, our model predicts 69,054 infected cases, while the actual number of infected cases is 69,317 (around 0.4% higher). By June 17, our model predicts 82,120 infected cases, while the actual number of infected cases is 85,823 (4.5% higher).

In July, many states paused plans to reopen, amid rising infected case counts. These pausing actions effectively slowed the spread of COVID-19, as can be clearly seen in the downward trend of the transmission rate in July in both Florida and Texas. Hence, a third change point is being detected in Texas and Florida, mainly related to this pausing.

The fitted number of infected cases and recovered cases are defined as for all t =2,…,T.

The observed and fitted number of infected cases are displayed in Figure 3. In summary, the piecewise constant SIR model (Model 1) significantly improves the performance in prediction of the number of infected cases compared with the baseline SIR model with fixed rates. Model 2.3 further improves the fit of the data for the California and Texas, due to the addition of a spatial smoothing effect. Finally, Model 3 provides further improvements, in particular for New York and Florida.

To determine the significance level of the spatial effect in model 2, we provide the estimate, p-value, and 95% confidence intervals for the parameter α in Table 15. We find that the influence of the infected or recovered cases in adjacent states is statistically significant (p-value < 0.05) for all states.

• Download figure
• Open in new tab

Fig 3:

Observed (black) and fitted (red) number of infected cases estimated by three models in five states.

Next, we assess the prediction performance for the three proposed models and the baseline SIR one. The out-of-sample mean relative prediction error (MRPE) is used as the performance, given by where is the predicted number of infected cases at time t; and I(t) is the observed number of infected cases at time t. The MRPE of R(t) can be obtained by respectively replacing the and I(t) with and R(t). We set the last two weeks of our observation period as the testing period and use the remaining time points for training the model. Note that predicted number of infected and recovered cases are defined as for all t =2,…,T, where the and are predicted by using the estimated parameters from the training period of the models.

The results of out-of-sample MRPE of I(t) and R(t) in the five states are reported in Tables 7 and 8. The calculated MRPEs of I(t) show that the model 3 which includes spatial effects and VAR temporal component outperforms the other models. Spatial smoothing itself (model 2) reduces the prediction error significantly in some states. For example, in Florida, the spatial smoothing reduced the MRPE by 80% when using model 2.2 (distanced-based weight) while in Oregon and Texas, the spatial smoothing reduced the MRPE (I) by more than 60% when using model 2.4 and model 2.3, respectively (similarity-based weight). In the state of New York, the prediction error reduction using spatial effects was not significant mainly due to the fact that its neighbors had different COVID-19 progression compared to New York. New York state was the hotspot at the beginning of the pandemic while the neighboring states had a surge in the number of cases much later. Finally, the reduction in MRPEs using model 3 (as reported in Tables 7 and 8) justifies the presence of the VAR component in the modeling framework. Additional results related to the VAR component (including the estimated auto-regressive parameters) are reported in Section 4 of the supplementary materials.

In Figure 4, we provide the predicted values in the last 2 weeks from all three different models in three states: New York, Florida and California. In both the daily infected cases ΔI(t) and the daily recovered cases ΔR(t), the predictions by model 2 and model 3 perform better than model 1 in all three states.

• Download figure
• Open in new tab

Fig 4:

Prediction of response variable y_h in three models. Black line: true value; orange: model 1; green: model 2; blue: model 3.

4.3. Results for Selected US Counties

We focus on the following counties in this analysis: New York City, King County in Washington, Miami-Dade County in Florida, Orleans and Jefferson parishes in Louisiana and Charleston, Greenville, Richland, and Horry in South Caroline. Note that The New York Times (2020) reports a single value for New York City which includes the cases in New York, Kings, Queens, Bronx and Richmond Counties.

The block size b_n is set to 7 for New York City, King and the South Carolina counties. For Miami-Dade, a smaller block size b_n = 5 is used, while for Louisiana counties a larger one b_n = 10. Further, for determining their neighbors displayed in Table 9, a threshold of 100 miles is used. Model 2.3 uses the top five counties in the corresponding state with the smallest similarity score, while Model 2.4 uses all counties in the given state. The resulting neighbors for Model 2.3 are displayed in Table 10.

Since the parameters obtained from (5) vary greatly, the 7-day moving average of transmission β(t) and recovery γ(t) rates is provided in Figure 9 and 10. As can be seen from those figures, our estimated piecewise parameters are very close to the real parameters in the 7-day moving average.

The statewide and countywide policy start dates and the detected change points are shown in Tables 11 and 12. A change point is detected from March to April for the first three counties. These change points coincide with the onset of “stay-at-home” orders and indicate a significant decrease in the transmission rate. A change point in June is detected in King county and Miami-Dade county, which may be related to the reopening plan after June as both change points correspond to a substantial increase in the transmission rate. Similar to Texas, a third change point is being detected in Miami-Dade county in July due to pausing in the reopening plans.

Next, using the proposed set of models we discuss two cases that provide interesting insights into policy questions.

Case I: Neighboring Parishes in Louisiana. A change point is detected in April for the Orleans and Jefferson parishes in Louisiana. This is primarily due to the fact that the virus circulated widely due to Mardi Gras activities in the end of February, that led to a high transmission rate in the March/early April period. However, the estimated transmission rate is much higher in Orleans, compared to Jefferson, which translates to a much higher number of infections in the former county. After reopening in early June, the transmission rate exhibited an increase (and hence a change point), only for Jefferson county. A possible explanation is that the high transmission rate in March in Orleans county led to a much larger number of infections, thus possibly conferring immunity to a larger portion of the population, while that was not the case in Jefferson county. Hence, in Jefferson’s case, the reopening led to a subsequent increase in the transmission rate. Note that seroprevalence survey data from the CDC² show a much higher percentage of immunity for Louisiana at the time of this writing, but probably not uniformly distributed across the state.

Case II: Neighboring Counties in South Carolina. We consider the following four neighboring counties (harleston, Greenville, Richland, and Horry), where mask wearing policies exhibited some variation in their timing, while Greenville county did not impose one. We can see that the transmission rates in these four counties share a very similar pattern with an upward trend in May and a downward trend from late June to August. Although all four counties show a significant decline in transmission rates in July, the mask mandate seems to lead to a reduction in transmission rates. Specifically, in Charleston at the change point in July the estimated transmission rate has a jump of 0.04, which is larger than the one in Greenville’s June change point of 0.02. The difference in jump size may relate to mandating the facial covering. What is difficult to delineate is whether the timing of the imposition of the mask mandate in three out of the four neighboring counties should have occurred earlier, so as the transmission rate not to exhibit such a large jump, especially in the most populous Charleston county.

Figure 11 shows the ACF plot of model 2’s residual time series for three counties. As can be seen from these plots, in Orleans parish, there is a significant temporal dependence among the time-lagged residual time series processes. This temporal dependence has been successfully captured by including the VAR(p) model for the residuals, as seen in the ACF plot of the new residuals, i.e. .

A more intuitive result is shown in Figure 5, where the fitted values of infected cases by model 2 and model 3 are much closer to the observed value than the predicted value by model 1. Specifically, Model 2.3 improves the fit of the data in all regions, due to the addition of a spatial smoothing effect. Model 3 further improves the fit of the data in New York City and Orleans parish.

• Download figure
• Open in new tab

Fig 5:

Observed (black) and fitted (red) number of infected cases estimated by three models in parishes/counties/cities.

We also provide the out-of-sample MRPE of I(t) and R(t) in Tables 13 and 14. The MRPE of I(t) results show that the piecewise constant model (model 1) performs the best in New York City, Orleans parish, while the piecewise constant model with spatial effect (model 2) performs the best in King County, Miami-Dade County, Jefferson parish and Richland. Adding the VAR(p) (model 3) performs the best in Charleston and Horry. The MRPE of R(t) results show that the piecewise constant model with spatial effect (model 2) performs the best in King County, Miami-Dade County, Greenville, Richland and Horry, while adding the VAR(p) (model 3) performs the best in New York City, Orleans, Jefferson and Charleston.

Finally, Figure 6 shows the predicted values in the last 2 weeks from all three different models. Here, we provide the prediction plots in three counties: Orleans, Charleston and Richland. In terms of the infected cases, the predictions by model 2 and model 3 perform much better than model 1 in all three counties.

• Download figure
• Open in new tab

Fig 6:

Prediction for the response variable y_h in three models. Black line: true value; orange: model 1; green: model 2; blue: model 3.