BayesSMILES: Bayesian Segmentation Modeling for Longitudinal Epidemiological Studies

3.1 Data notations

During a pandemic such as COVID-19, the most accessible and complete data are the daily reported numbers on confirmed cases and deaths. Suppose N is the total population size in a given region. Let C = (C₁, …, C_T) and D = (D₁, …, D_T) be the sequences of cumulative confirmed case and death numbers observed at T successive equally spaced points in time (e.g. day), where C_t and D_t ∈ ℕ for t = 1, …, T. For a region of which recovery cases are closely monitored day by day, we use E = (E₁, …, E_T) to denote the sequence of cumulative recovery case numbers. Thus, due to the compositional nature of the basic SIR model, the three trajectories can be constructed as S = (S₁, …, S_T) with S_t = N − C_t, R = (R₁, …, R_T) with R_t = D_t + E_t, and I = (I₁, …, I_T) with I_t = N − S_t − R_t = C_t − D_t − E_t. For a region of which recovery cases do not exist or under-reported, we consider both R and I as missing data and reconstruct these two sequences according to the last equation of (1) with a pre-specified constant removal rate γ. Specifically we set I₁ = C₁ and R₁ = 0, and generate R_t = R_t−1 + rI_t−1 and I_t = I_t−1 +(C_t − C_t−1) − (R_t − R_t−1) from t = 2 to T sequentially, where. : [0, ∞) → N denotes the ceiling function. For the choice of γ, we suggest estimating its value from publicly available reports in some region where confirmed, death, and recovery cases are all well-documented, or from prior epidemic studies due to the same under-reporting issue in actual data. Lastly, let Y₀ be the initial value and be the lag one difference of a sequence Y, where and each following entry , i.e. the difference between two adjacent observations. Note that Y could be any of those time-series data introduced above, i.e. C, D, E, S, I, and R.

3.2 Modeling epidemic dynamics via a modified stochastic SIR model

An SIR model has three trajectories, one for each compartment. The compositional nature of the three trajectories, i.e. S(t) + I(t) + R(t) = N, implies that we only need two of them to solve the ODEs as shown in (1). As mentioned previously, assuming S(t) ≈ N for all t results in dI(t)/dt = (β − γ)I(t) and further leads to an exact solution: I(t) = I(0) exp((β − γ)t). For modeling daily reported infectious data I, we utilize its discrete-time version, with a time-varying rate β_t to account for the transmissibility changes of the disease. For simplicity’s sake, we assume a constant removal rate γ. Based on (2), we introduce a probabilistic model, which approximately mimics the dynamics of the deterministic standard SIR model as shown in (1) during the onset phase of a pandemic. Specifically, we suppose the infectious population size at time t is sampled from a Poisson model, where α_t = I₀ exp ((β_t − γ)t)/N is a redefined time-varying transmissibility parameter that depends on the initial infectious population size I₀, the transmission rate β_t, the removal rate γ, and any latent factors (e.g. public health intervention, social behavior, virus mutation, healthcare quality, etc.) that may affect the disease transmissibility. This model automatically accounts for measurement errors and uncertainties associated with the counts. Note that (3) can be generalized to a negative binomial (NB) model, i.e. I_t|α_t ∼ NB(Nα_t, ϕ_I) if needed, where ϕ_I is a dispersion parameter to account for over-dispersion that might be observed in I. Here we use NB(µ, ϕ), µ, ϕ > 0 to denote an NB distribution with expectation µ and variance µ + µ²/ϕ.

3.3 Detecting multiple change-point via a Poisson segmented regression model

Our multiple change-point detection builds upon the above modified stochastic time-variant SIR model (3) with stationary transmissibility α_t in a certain segment. Particularly, we assume that β_t only changes at certain time points. Identifying those change points is of significant importance, as it not only enables us to characterize the temporal dynamics of the pandemic but also helps policymakers evaluate the effectiveness of the past and ongoing mitigation and intervention strategies.

In this paper, the change points are defined as those discrete time points that significantly alter the transmission rate β_t between two adjacent segments, while assuming the removal rate γ is constant across all time points. We introduce a latent binary vector γ = (γ₁, …, γ_T), γ_t ∈ {0, 1}, with γ_t = 1 if time point t is a change point and γ_j = 0 otherwise. We restrict γ₁ = 1, although the first time point is not a change point. Those points with γ_j = 0 can be partitioned into segments bounded by two adjacent change points. Thus, we use another vector z = (z₁, …, z_T), z_i ∈ {1, …, K} to reparameterize γ, where z_t = k if time point t is in segment k, that is, between the k and (k + 1)-th change points. The total number of change points excluding the first time point is K − 1. Mathematically, z is the cumulative sum of γ, i.e. , while γ is the lag one difference of z, i.e. γ_t = z_t − z_t−1 with γ₁ = 1. Figure 1 shows the representation of γ and z for a simulated time-series dataset (T = 10) with two change points.

• Download figure
• Open in new tab

Figure 1:

An example of time-series data (T = 10) with two change points (K = 3) and its associated parameterizations in terms of γ and z, respectively. Red triangles are change points, while black circles are not. Note that the first time point is not a change point, although it is defined as a “change point” for convenience’s sake.

To infer γ or z given the sequence of infectious population size I, we adopt a Poisson segmented regression framework, which can be written as, where x_t is a p-dimensional column vector of covariates that includes a scalar of one for the intercept, time t, and p − 2 explanatory variables observed at time t if applicable. Those explanatory variables could contain the number of tests, weather information, or other temporal information that are important to adjust for during a longitudinal epidemiological study. b_k is a p-dimensional column vector of segment-specified coefficients that includes an intercept representing the proportion of infectious people at logarithmic scale, i.e. b_k1 = log(I₀/N), in segment k, and a slope accounting for the segment-varying transmission rate, i.e. b_k2 = β_k − γ. For simplicity’s sake, we assume that ϵ_t is an independent and identically distributed random variable from a zero-mean normal distribution with segment-specific variance . To ensure the identifiability of all model parameters, we try a set of considerably small values for and employ a robust cross validation method called Pareto-smoothed importance sampling leave-one-out (PSIS-LOO) cross validation to determine the best choice (Vehtari et al., 2017).

Let α_k be the sequence of all α_t’s in segment k, i.e. , where we denote as the time point corresponding to the k-th change point and as the number of time points in segment k. We can rewrite the second stage in (4) as , where is an identity matrix and X_k can be explicitly written as We assume a zero-mean multivariate normal distribution for b_k, that is, b_k ∼ MN(0_p, H), where 0_p is an p-by-1 all zero column vector and is set to be a diagonal variance-covariance matrix. For a non-informative setting, we recommend choosing a large value for each . Through this prior specification, we are able to marginalize out the nuisance parameter b_k and obtain . Consequently, we can write the collapsed model of (4) as To complete the model specification, we impose an independent Bernoulli prior on γ as , where ω is interpreted as the probability of a time point being a change point a priori. We further relax this assumption by allowing ω ∼ Be(a_ω, b_ω) to achieve a beta-Bernoulli prior. In practice, we suggest a constraint of a_ω + b_ω = 2 for a vague hyperprior of ω (Tadesse et al., 2005). In addition to that, we make the prior probability of γ equal to zero if two adjacent time points are jointly selected as change points. In other words, a segment should consist of at least two time points.

3.4 Estimating basic reproduction numbers via a stochastic SIR model

Given the multiple change points γ, we estimate the basic reproduction number R₀ for each segment via a stochastic version of the standard SIR model as shown in (1), which only needs the cumulative confirmed case numbers C. This is because recovery data only exist in few states in the U.S., which makes both model inference and predictions infeasible. This model considers both of the removed and actively infectious cases as missing data and mimic their relationship in spirit to some compartmental models in epidemiology. Specifically, we assume the number of new removed cases at time t, i.e. , is sampled from a Poisson distribution with mean γI_t−1, that is, , where γ should be pre-specified. Based on this simplification, we rewrite the discrete version of the first equation in (1) as, resulting in We further assume the new case number observed at time t, i.e. , is sampled from an NB model, as it automatically accounts for measurement errors and uncertainties associated with the counts. Following most epidemiological models, we assume this stochastic process is a Markov process, where the present state (at time t) depends only upon its previous state (at time t − 1). For the prior distribution of the dispersion parameter ϕ, we choose a gamma distribution, ϕ ∼ Ga(a_ϕ, b_ϕ). We recommend small values, such as a_ϕ = b_ϕ = 0.001, for a non-informative setting. This model, on average, mimics the epidemic dynamics and is more flexible than those deterministic epidemiological models. We assume β comes from a gamma distribution with hyperparameters that makes both mean and variance of the transformed variable ℛ₀ = β/γ equal to 1.

4.1 MCMC algorithms for detecting multiple change-point

Our primary interest lies in identifying the change-point locations via the vector γ based on the actively infectious cases I. According to Section 3.3, the full data likelihood and the priors of the change-point detection model are written as, where we use the superscript (γ) to denote the transmissibility parameters α_k, the design matrix X_k, and the number of time points n_k in segment k under a specific configuration of γ. At each MCMC iteration, we perform the following steps.

4.2 MCMC algorithms for estimating basic reproduction numbers

Once the change points are determined, we aim to estimate the basic reproduction numbers ℛ₀′s across different segments and quantify their uncertainties based on the cumulative confirmed cases C only. According to Section 3.4, the full data likelihood and the priors of the stochastic SIR model are written as, where β = (β₁, …, β_K) and ϕ = (ϕ₁, …, ϕ_K), i.e. the collections of transmission and dispersion rates of all segments. For the hyperparameters, we set a_β = 1 and b_β = 1/γ so that both of the expectation and variance of the basic reproduction number ℛ₀ = β_k/γ are equal to one. We specify a_ϕ = b_ϕ = 0.001 to obtain a weakly informative gamma prior for all dispersion parameters. With a pre-defined removal rate γ, we propose the following updates in each MCMC iterations.

4.3 Posterior inference

We explore posterior inference on the parameters of interest by postprocessing the MCMC samples after burn-in. We start by obtaining a point estimate of the change-point indicator γ by analyzing its MCMC samples {γ^(b), …, γ^(B)}, where b indexes the MCMC iteration after burn-in. One way is to choose the γ corresponding to the maximum-a-posteriori (MAP), The corresponding can be obtained by taking the cumulative sum of . An alternative estimate relies on the computation of posterior pairwise probability matrix (PPM), that is, the probabilities that time point t and t^′ are assigned into the same segment, . This estimate utilizes the information from all MCMC samples and is thus more robust. After obtaining this T -by-T co-clustering matrix denoted by P, a point estimate of z can be approximated by minimizing the sum of squared deviations of its association matrix from the PPM, that is, The corresponding can be obtained by taking the difference between consecutive entries in and set the first entry to one. To construct a “credible interval” for a change point, we utilize its local dependency structure from all MCMC samples of γ that belong to its neighbors. Due to the nature of the MCMC algorithm described in Section 4.1, if a time point t is selected as a change point, i.e. γ_t = 1, then its nearby time points must not be a change point. Thus the correlation between the MCMC sample vectors and tends to be negative when u is small. Thus, we define the credible interval of a change point as the two ends of all its nearby consecutive time points, of which MCMC samples of γ are significantly negatively correlated with that of the change point. This could be done via a one-sided Pearson correlation test with a pre-specified significant level, e.g. 0.05. Although quantifying uncertainties of change points is not formal, it performs very well in the simulation study and yields reasonable results in the real data analysis.

We constructed the confidence interval for a selected change point t by evaluating its MCMC samples . In particular, according to the swap step that updates γ in the MCMC algorithm, the change point was typically switching between the selected time point t and its close neighborhood in the MCMC iterations. This is due to the fact that the change in patterns between segments may not be abrupt, and therefore the time points closed to time t could also be accepted as a change point in an iteration. Hence, for each change point with γ_t = 1, we calculated the correlation between and Here, we set different values for Δt (e.g., Δt = −3, −2, …, 1, 2) such that the collection of t+Δt included a reasonable neighborhood for time t. Then, we calculated the correlation between and . Given the fact that only one change point would appear in the neighborhood region, vector tended to have a negative correlation with .

Once the change points are determined, an approximate Bayesian estimator of the transmission rate β_k for each segment k can be simply obtained by averaging over all of its MCMC samples, . In addition, a quantile estimation or credible interval can be obtained. Lastly, we summarize the basic reproduction number in each segment k as .

4.4 Prediction

Conditional on the multiple change-point locations, we can predict the cumulative or new confirmed cases at any future time T_f by Monte Carlo simulation based on the information in the last segment K only. Specifically, from time T + 1 to T_f, we sequentially generate Then, both short and long-term forecasts can be made by summarizing the (T_f − T)-by-B matrix of MCMC samples. For instance, the predictive number of cumulative and new confirmed cases at time T + 1, in average, are and , respectively.

5.1 The generative model

The three trajectories S, I, and R with length T = 120 were generated in the following ways. We first divided the T = 120 time points into K = 4 segments with the same length, that is, the true change points were t = 31, t = 61, and t = 91. In other words, the first 30 entries of z equal to 1, the following 30 entries equal to 2, the next 30 entries equal to 3, and the last 30 entries equal to 4. To mimic the disease transmissibility dynamics across different segments, we chose segment-varying transmission rates β_k while fixing the removal rate γ = 0.03. We considered four scenarios of the set (β₁, β₂, β₃, β₄), corresponding to 1) ℛ₀ = (3.0, 1.2, 2.0, 0.8); 2) ℛ₀ = (3.0, 2.3, 1.5, 0.8); 3) ℛ₀ = (3.0, 1.8, 0.8, 1.6); 4) ℛ₀ = (3.0, 2.0, 1.1, 0.5). Then based on the stochastic version of the standard SIR model, we sampled S_t and R_t from negative binomial (NB) distributions, and obtained I_t, sequentially from t = 1 to T through where N = 1, 000, 000, the initial I₀ = 100 and R₀ = 0, and the dispersion parameters ϕ_S = ϕ_R = 10. Note that the generative scheme was with an NB error structure, which was different from our model assumption based on a Poisson error structure. We repeated the above steps to generate 50 independent datasets for each setting of ℛ₀. Figure 2 displays the temporal patterns of the simulated infectious counts I for the four scenarios.

• Download figure
• Open in new tab

Figure 2:

Simulation study: The simulated actively infectious data I under the four scenarios. Each curve represents a replicated sequence of I under a scenario. The red dashed lines mark the true change-point locations. The blue curve under Scenario 4 was randomly chosen for evaluating the model fitting, of which results are shown in Figure 3.

5.2 Evaluation criteria

To evaluate the multiple change-point detection, we may rely on either the binary change point indicator vector γ or the time point allocation vector z. For the choice of γ, a change point is considered to be correctly identified if its location is within a local window of the true position (Killick and Eckley, 2014). The selection of the window size is ad hoc and may bias the evaluation.

In addition to that, since change points and non-change points are usually of very different sizes, most of the binary classification metrics are not suitable for model comparison here. Thus, we chose those metrics that quantify the agreement between the true and estimated allocation vectors, i.e. z and . The two classic performance metrics for the analysis of clustering results are the adjusted Rand index (ARI) and mutual information (MI), proposed by Hubert and Arabie (1985) and Steuer et al. (2002), respectively. ARI is the corrected-for-chance version of the Rand index (Rand, 1971), as a similarity measure between two sample allocation vectors. Let ; and be the number of pairs of time points that are a) in the same segment in both of the true and estimated partitions; b) in different segments in the true partition but in the same segment of the estimated one; c) in the same segment of the true partition but in different segments in the estimated one; and d) in different segments in both of the true and estimated partitions. Then, the ARI can be computed as The ARI usually yields values between 0 and 1, although it can yield negative values (Santos and Embrechts, 2009). The large the index, the more similar between z and , thus the more accurately the method detects the actual times at which change points occurred. An alternative metric choice is MI, which measures the information about one variable that is shared by the other (Steuer et al., 2002). Let be the number of time points shared between the k-th segment in the true z and the k^t-the segment in the estimated one . Then, the MI can be computed as where is the number of segments and ’s are the segment lengths for segment in . It yields non-negative values. The larger the MI, the more accurate the partition result.

To quantify how well a method estimates the dynamic transmissibility across different segments, we used the root mean square error (RMSE) that measures the deviation between the true and estimated values of ℛ₀ over all T time points: A smaller value of RMSE indicates a more accurate estimation of ℛ₀’s.

5.3 Results

As for the MCMC setting of change point detection, we set 40, 000 MCMC iterations and discarded the first half as burn-in. We adopted the weakly informative setting by setting a_ω = 0.1 and b_ω = 1.9 in the Beta-Bernoulli prior for the change-point indicator vector γ. We set with and as the covariance matrix in the prior distribution of b_k’s. Finally, we let take ten equally spaced values ranging from 0.0001 to 0.01 at the logarithmic scale (base 10) in the PSIS-LOO cross validation. In fitting the stochastic SIR model, we set 100, 000 MCMC iterations with the first half as burn-in. As suggested in Waqas et al. (2020), the value of removal rate γ could be estimated by for each simulated dataset. Then, we set a_β = 1 and b_β = 1/γ so that both of the expectation and variance of the basic reproduction number ℛ₀ = β_k/γ are equal to 1.

We first checked the performance of the BayesSMILES on a single simulated dataset, which was randomly selected from the 50 replicates in Scenario 4 (marked as the blue line in Figure 2). Figure 3(a) demonstrates the change-point detection result based on the Poisson segmented regression model. The red dashed and the blue solid lines represent the true and the estimated change point locations, respectively, while the gray ribbons represent the 95% confidence intervals for those identified change points. As we can see, the BayesSMILES successfully detected the three true change points in general, as each of the 95% confidence intervals covered the truth. The resulted values of ARI and MI were 0.93 and 1.28, respectively. Later on, the stochastic SIR model introduced in Section 3.4 was then fitted to quantify the disease transmissibility in each segment bounded by the identified change points. Figure 3(b) shows the posterior distributions of ’s for from their MCMC samples. The red dashed and blue solid line pinpoints the true and posterior mean of ’s, while the two black solid lines mark the boundary of their 95% credible intervals. Clearly, those true values were within their corresponding 95% credible intervals. The final RMSE for ℛ₀ estimation was 0.38 for this single simulated dataset.

• Download figure
• Open in new tab

Figure 3:

Simulation study: The model fitting results based on a randomly selected simulated dataset (see the blue curve under Scenario 4 in Figure 2). (a) The locations of change points (blue solid lines) estimated from the posterior pairwise probability matrix (PPM) and their credible intervals (gray ribbons). The red dashed lines mark the true change-point locations; (b) The posterior distributions of ℛ_0k’s for k = 1, 2, 3, 4 estimated from the segmented time-series data, given the three identified change points as shown in (a). The red dashed and blue solid lines are the true and estimated values of ℛ_0k’s, respectively. The two black solid lines are the lower and upper bounds of the 95% credible intervals.

To our best knowledge, there is no method like the BayesSMILES that can detect latent change points while characterizing the transmission dynamics through an SIR model. Thus, in setting up a comparison study, we therefore considered a two-stage approach that first identifies multiple change-point of time-series data based on a likelihood based framework, and then estimates the basic reproduction numbers between each pair of nearby change points, following the stochastic SIR model introduced in Section 3.4. The alternative change point model assumes time points within one segment follow a normal distribution with distinct mean and/or variance from its nearby segments (Hinkley, 1970; Jen and Gupta, 1987), and it uses likelihood ratio test (LRT) to detect multiple change points. An algorithm named binary segmentation (Edwards and Cavalli-Sforza, 1965; Sen and Srivastava, 1975) is commonly used to compute the test statistics for the LRT with high efficiency (Killick et al., 2012). In our case, to detect multiple change points using this alternative approach named likelihood ratio test with binary segmentation (LRT-BinSeg), we input the logarithmic scale of I into the function cpt.meanvar in the related R package changepoint (Killick and Eckley, 2014) for each of the simulated datasets. We specified to use the binary segmentation algorithm with up to 5 of change points to search for. Note that this restriction was not applicable to the alternative algorithms provided in the changepoint package. In practice, we found that the alternative algorithms tended to over-select the number change points.

Figure 4(a) and (b) exhibit the change-point detection performances for the four scenarios of ℛ₀. Our BayesSMILES performed much better than the likelihood ratio test with the binary segmentation algorithm (LRT-BinSeg) with respect to multiple change-point detection under both performance metrics, ARI and MI. For instance, the ARI by the BayesSMILES increased 39.29% to 122.16% over the LRT-BinSeg among the four scenarios, while the growth in MI could be up to 60.54%. Figure 4(c) compares the ability to capture the transmission dynamics in terms of RMSE, which depends on the change-point detection accuracy. As expected, Our BayesSMILES yielded smaller RMSE values across all scenarios since its identified change point locations were more accurate. In all, the simulation study demonstrated that the strengths of BayesSMILES.

• Download figure
• Open in new tab

Figure 4:

Simulation study: The violin plots of (a) adjusted Rand index, (b) mutual information, and (c) R0 root mean square error from 50 replicated datasets generated under the four scenarios. Red and blue violins correspond to the results obtained by the BayesSMILES and LRT-BinSeg.

6 Analysis of COVID-19 Data

In this section, we applied the BayesSMILES to the COVID-19 daily report data provided by JHU-CSSE (COVID-19 Data Repository). Several recent COVID-19 studies also based their analyses on this resource (Dong et al., 2020; Zhou and Ji, 2020; Toda, 2020). We first performed a preprocessing step to ensure the quality of the infectious data I for the model fitting. Due to the fact that the recovery cases did not exist for some U.S. states, we treated I and R as missing data and reconstructed the two sequences according to the process described in Section 3.1. The observed cumulative confirmed case numbers C were collected for each U.S. state starting from an early stage of the pandemic outbreak. In particular, we chose the starting time for each state as when there were least ten confirmed COVID-19 cases for that state. We also set the removal rate γ = 0.1 as suggested by Pedersen and Meneghini (2020) and Weitz et al. (2020). Due to the fact that different states could have different starting times, we further trimmed the sequences I and R for each state based on the latest starting time available. Finally, we set March 22, 2020, as the new starting time (t = 1) for all 50 states, and let July 19, 2020, be the last observed time point (t = 120). We used the same hyperparameter and algorithm settings as described in Section 5.3.

6.1 Detecting change points for U.S. states

We limit our analysis to four U.S. states with the highest cumulative confirmed cases as of July 19, 2020, to keep the paper in a reasonable length. They are New York, Texas, California, and Florida. The results for the 46 remaining states are available in https://shuangj00.github.io/BayesSMILES/ (see details in Section 8). Figure 5 displays the detected change points, as well as the estimated basic reproduction number ℛ₀’s cross segments, for the four states. The associated confidence interval to each identified change point is represented by a gray ribbon. In general, those change points detected by the BayesSMILES indeed captured the important COVID-19 events that might affect the transmission rates. For instance, some change points reflected the positive effects of the preventative strategies such as lockdown, while others explained the “bounce back” in confirmed cases after the reopening. Table 2 lists the change point locations and their potentially related events for the four states.

• Download figure
• Open in new tab

Figure 5:

Case study: The change point detection and basic reproduction number estimation for the states of New York, Texas, Florida, and California. The black circles are the actively infectious case numbers divided by the total population (in thousands) at logarithmic scale, i.e. log(It/N). The black solid lines pinpoint the change point locations, with the associated gray ribbons indicating the credible intervals. The red dashed lines describe the variation in the basic reproduction numbers ℛ₀ across segments. The posterior means and the 95% credible intervals for ℛ_0k’s are given by red numbers.

In New York, the first change point was estimated to be March 28. We estimated the posterior mean of the basic reproduction number decreased from 2.24 (between March and to March 27) to 1.63 (between March 28 and April 8). Notably, March 28 was the date when the Centers for Disease Control and Prevention (CDC) issued a 14-day domestic travel advisory for non-essential persons, which presumably alleviated the situation for the populated states such as New York. The second change point appeared around April 9, and the R₀ of the third segment dropped to 0.98 with a 95% credible interval of [0.76, 1.25]. This matched the exact day when the New York state posted its first drop of the ICU admissions since the COVID-19 outbreak began. The third change point was around April 27. Though there was no direct intervention issued in late April, we noticed that the mayor of New York City announced that all the major events had been canceled starting from April 20. This action could bring a positive effect in controlling the outbreak, and our estimation from the SIR model suggested a further decrease in the basic reproduction number down to 0.66 with a 95% credible interval of [0.54, 0.81]. We observed another change point around June 18, which was close to the Phase II reopening of New York state on June 22. During the Phase II reopening, restaurants were allowed for outdoor dining, stores opened for in-person retail, and more services resumed operational under strict limitations. Thus, we saw a little “bounced back” in ℛ₀ from 0.66 to 0.82. The last change point was on June 29. As expected, the basic reproduction number increased to 1.04 with a 95% credible interval of [0.84, 1.29] in the last segment. Although there was no public announcement around June 29 with a confidence interval from June 28 to July 4, we suspect that the increased social interaction during the Independence Day long weekend (between July 3 and July 5) could be responsible for the increase of the transmission dynamics.

In Texas, there were five change points detected. The first change point was estimated to be March 28, the same day as the first one for New York state. Due to the similar reason, the policy of mandatory 14-day quarantines for travelers entering Texas could bring a decrease in terms of the basic reproduction number (decreased from 2.97 to 2.07). The second change point was around April 9 with a further drop of ℛ₀ to 1.14 with a 95% credible interval [0.96, 1.35]. We found that the Texas Governor had extended the state’s disaster declaration for an additional 30 days on April 12. The extension aimed at protecting the health and safety of Texans by ensuring adequate capabilities to support communities. Organizations such as the State Operations Center and the Strategic National Stockpile would continuously supply the state government with resources needed to protect residents. May 25 was detected as the third change point, and it was the first time that ℛ₀ increased after the two drops. The estimated basic reproduction number was 1.29 with a 95% credible interval [1.02, 1.62]. This increase appeared around May 25 could be due to the Governor’s updated executive order issued on May 26 that allowed additional services and activities to open for phase II reopening. The next change point was around June 16, and ℛ₀ further increased to 1.72 with a 95% credible interval [1.40, 2.11]. According to the prediction reported by UT Austin’s COVID-19 Modeling Consortium at the end of May, there might be a significant increase in the number of cases and hospitalizations beginning mid-June (News from kxan). Here, the change point location, as well as the increased basic reproduction number, were consistent with the results of this report. The last change point was around June 28 with an estimated decrease in ℛ₀ to 1.42 with a 95% credible interval [1.17, 1.72]. Notably, the Texas Governor issued multiple executive orders around late June to early July to mitigate the disease spreading. For instance, the executive order on June 26 reemphasized the limited occupancy for all business establishments in Texas. According to an executive order on July 2, all Texans were required to wear a face-covering in public spaces in counties with 20 or more positive COVID-19 cases. On the same day, the Governor announced an update regarding the executive order on June 26 with additional measures to slow down the spreading of COVID-19.

In Florida, the first estimated change point was April 3. It was two days after the statewide stay-at-home order for Florida. We estimated that the basic reproduction number decreased from 2.70 to 1.28 after the change point. The second change point appeared around the middle of April. Starting from April 13, some counties such as Osceola county started the requirement of wearing a face covering while in public. It could explain the reason why we observed a slight decrease in ℛ₀, from 1.28 to 0.92 with a 95% credible interval [0.73, 1.15]. The next change point was located around May 13, and ℛ₀ in this new stage went above 1 again, with a posterior mean of 1.19 and a 95% credible interval [0.96, 1.50]. We noticed that Florida entered the phase I reopening on May 18, which could lead to the “bounced back” situation. The fourth change point was around June 7, two days after the phase II reopening in Florida. Changes in the phase II reopening included Universal Orlando opened the parks to the general public for the first time in months, and we observed that ℛ₀ increased again to 1.81 with a 95% credible interval [1.50, 2.19]. In the last segment (after June 27), our result revealed a slight drop in the basic reproduction number from 1.81 to 1.45. This change was potentially related to the consequence of requiring facial coverings in the four most populated cities in Florida: Tampa, Orlando, Miami, and Jacksonville. The face mask mandates went into effect for the four cities starting from June 19, 20, 25, and 29, respectively. Therefore, the drop in the transmissibility at the end of June may be explained by the effectiveness of wearing face masks as a non-pharmaceutical practice.

In California, we detected two change points. California was the first state to announce lockdown in the COVID-19 pandemic and its stay-at-home order became effective on March 19. Our results in change point detection could miss these early actions since the data we analyzed started from March 22. The first selected change point was on April 5, with the value of the basic reproduction number decreased dramatically when transitioning to the second segment (from 2.20 to 1.20). The second change point was on June 17, and we saw that ℛ₀ increased to 1.35 in the last segment with a 95% credible interval [1.12, 1.62]. According to California Governor, higher-risk businesses and venues (e.g. movie theaters, bars, gyms) were allowed to reopen with restrictions on June 12. Hence, increasing in the basic reproduction number could be the consequence of reopening. This observation was also observed in New York and Texas.

6.2 Clustering U.S. states based on their change-point locations

We applied BayesSMILES on all 50 U.S. states. Based on the results, we seek to derive an overall picture of the COVID-19 dynamics across states. We summarized the temporally detected change points of the 50 states into common patterns, and then we labeled each state by matching its specific change point pattern to the common patterns. In particular, for each state, we calculated the marginal posterior probability of inclusion (PPI) for all time points, where the PPI for a time point t was calculated based on the B of MCMC samples after burn-in: . Then we obtained the vector . Each entry in p^PPI is a value between 0 and, 1, representing the proportion of time t selected as a change point among all iterations. Next, we computed the overall pattern by averaging over the vector p^PPI across 50 states. We noticed that some time points were rarely or never selected as change points. That naturally hinted us to group the time points. To illustrate this, we trimmed the top 20% values of for each time t (Figure 6). The trimming step provided a clear pattern and highlighted the groups of dates that were commonly identified as change points. We observed three time spans as shown in Figure 6: March 27 - April 11, May 1 - May 10, and May 22 - July 3. For each state, we defined its cluster label based on the corresponding change point detection results. If this state had at least one change point (including the confidence interval) between March 27 - April 11, the first element in its cluster label is “Change”. Otherwise, the first element in the group label was set to “Stable”. We repeated the same process to determine the second and third elements of the class label for each state. In the end, each state was assigned to a cluster label from “Change-Change-Change”, “Change-Stable-Change”, or “Stable-Change-Change”.

• Download figure
• Open in new tab

Figure 6:

Case study: The averaged marginal posterior probability of inclusion (PPI) for each time point to be selected as a change point over all 50 U.S. states, after trimming the top 20% PPI values. The black dashed lines partition the whole time range into three segments: March 27 - April 11, May 1 - May 10, and May 22 - July 3.

The map in Figure 7 colors each of the 50 states based on its cluster label, where green, yellow, and pink correspond to temporal patterns “Change-Stable-Change”, “Change-Change-Change”, and “Stable-Change-Change”, respectively. Interestingly, three out of the four states we have analyzed, New York, Texas, and California, belonged to the same category, i.e. “Change-Stable-Change”. Several other states that are also in this category include Georgia, Arizona, North Carolina, and Louisiana. All of these states were in the top ten states with the most COVID-19 confirmed cases. We noticed that the phase I statewide reopening for all these states occurred in mid-May (May 15 for Georgia, May 13 for Arizona, May 8 for North Carolina, May 15 for Louisiana). Therefore, our model did not report any change points for these states between May 1 and May 10. The rest states in the 10 states with the most COVID-19 confirmed cases, including Florida, Illinois, and New Jersey, were labeled as “Change-Change-Change”, and all of them had a change point between May 1 and May 10. As discussed in Section 6.1, Florida had a change point around May 13 with a confidence interval [May 8, May 18]. According to the Executive Order 2020-32 issued by the Illinois governors, the state entered the phase II reopening starting on May 1 with a modified stay-at-home order. For New Jersey, the statewide state-at-home order was not lifted until June 9. However, our model suggested a change point around the end of April with a confidence interval [April 25, May 3] with a drop in R₀ from 1.11 to 0.68 (details available at https://shuangj00.github.io/BayesSMILES/). We noticed that on May 3 the New Jersey Governor announced a multi-state agreement to develop a regional supply chain for personal protective equipment, other medical equipment and testing. This joint-state protective measure allowed for efficient delivery and reliability of medical equipment for states, and therefore best utilized the life-saving resources in the face of the COVID-19 outbreak.

• Download figure
• Open in new tab

Figure 7:

Case study: The temporal patterns of the COVID-19 transmission dynamics based on change points across the 50 U.S. states. Green, yellow, and pink correspond to “Change-Stable-Change”, “Change-Change-Change”, and “Stable-Change-Change” patterns, respectively.

6.3 Predicting new confirmed cases for U.S. states

Short-term forecasting of the new daily confirmed COVID-19 cases at a future time T_f is important for designing public policies. Here, we used the estimation results from BayesSMILES to predict the new confirmed cases for the four states that we analyzed: New York, Texas, Florida, and California. As suggested in Section 6.1, the basic reproduction number ℛ₀ for one state could go above or below the threshold value of 1 at different disease stages. Hence, for short-term forecasting, the predictions given by the BayesSMILES were based on the last segment’s information, which ensured that the predictions fully reflected the most recent disease characteristics. In contrast, we evaluated an alternative approach that ignored the underlying disease dynamics and used all the data to fit a single SIR model for the prediction. We named this method FullDataSIR. Here, we illustrated the advantages of incorporating the change point information into prediction. In particular, we compared the prediction accuracy of the two methods: the BayesSMILES and FullDataSIR.

Figure 8 showed the true value of the new daily confirmed COVID-19 cases and the two types of predictions for the four states. For each state, the 7-day forecast was from July 20 to 26. First of all, we observed that the predicted mean by the FullDataSIR model tended to be larger than that from the BayesSMILES for the same day. This was because the data in the early stage (in our case, from late March to early April) corresponded to a much larger basic reproduction number, as compared to the data from the late stage, and the BayesSMILES would ignore the early stage observations in making the prediction. Second, the closer the predicted value to the true value on each day, the better the prediction performance. We thus quantified the prediction accuracy using the mean absolute percentage error (MAPE). The MAPE for the 7-day forecast was defined as where was the observed new confirmed cases at future time T_f, and was the corresponding prediction. The smaller the MAPE value, the more accurate the prediction.

• Download figure
• Open in new tab

Figure 8:

Case study: The 7-day forecast (between July 20 and 26) of the daily confirmed COVID-19 case numbers for the states of New York, Texas, Florida, and California. The red and blue circles and bars are the predictive means and 95% intervals by the BayesSMILES and FullDataSIR, respectively. The black thick lines indicate the observed truth.

We calculated the MAPE for the 7-day forecast of the four states analyzed in Section 6.1, and Table 1 summarized the results. For the first three states, the MAPEs from the BayesS-MILES were much smaller than the MAPEs from the FullDataSIR model, suggesting a better performance of BayesSMILES. For California, the conclusion was reversed since the FullDataSIR model showed a slightly lower MAPE. However, if we look at Figure 8(d), the real data (black horizontal lines) and the predictions from the BayesSMILES (red circles) both showed a down-going trend. In this case, the consistency in trend could be more important to evaluate two types of predictions, as compared to the slight difference in MAPEs. Lastly, for most of the days shown in Figure 8, the 95% credible interval by BayesSMILES tended to be much narrower than the interval given by the FullDataSIR model under the same day, especially for Texas, Florida, and California. Combined with a smaller value of MAPE, the BayesSMILES was able to provide a more accurate and precise prediction as compared to the alternative approach. In all, incorporating the change point location information helped improve the prediction accuracy for the 7-day forecast of new confirmed COVID-19 case numbers.