Changes in Reproductive Rate of SARS-CoV-2 Due to Non-pharmaceutical Interventions in 1,417 U.S. Counties

Characterizing Groups of Similar Counties as Clusters

Since our hypothesis is that local conditions affect the spread of COVID-19, we differentiate among groups of counties using a data-driven approach known as clustering. Each group – or cluster – of counties is characterized by having similar demographic and socioeconomic qualities. For instance, cluster 1 consists of low-population, mostly rural counties with little public transit capacity and the lowest median household income. Cluster 2 and 3 are similar in size, having a mean population of 45,000 and 52,000 respectively, but cluster 2 includes higher-income, suburban areas where-as cluster 3 has lower income areas with a large land area. Cluster 4 consists of the densest metropolitan areas with a high proportion of 18- to 65-year-olds, high household income, a high public transit score, and small land area. Finally, cluster 5 has the highest mean population besides cluster 4, but is less densely populated, has poor public transit, and a lower household income on average. It is important to note that although we refer to these clusters with numbers 1–5, these are merely labels returned by the clustering algorithm without any meaning inherent to the ordering Figure 1. shows our clustering for all U.S. counties with this data available, including those without any incidence of COVID-19.

• Download figure
• Open in new tab

Figure 1:

Cluster labels based on demographic and socioeconomic conditions are used to aggregate data and specialize epidemiological models. Here, one can see how cluster 1 and 3 primarily cover rural areas, while clusters 5, 2, and 4 consist of increasingly urban counties.

Once we have identified the clusters, we infer the reproductive rate of SARS-CoV-2 over time. by jointly optimizing the parameters of a Bayesian mechanistic model on each cluster. This process is described in greater detail in Section 4. Figure 2 shows an example of this progression for counties ordered by their public transportation use.

• Download figure
• Open in new tab

Figure 2:

Relationship between public transit capacity and the time-dependent reproductive rate of SARS-CoV-2 for U.S. counties. Colors correspond to clusters as in Figure 1. In the cluster of counties that rely least on public transport, shown in yellow, we observe a substantial drop in reproductive rate already on 3/15/20, possibly indicating a strong effect of advisory NPIs that were in effect at that time.

We observe that although the basic reproductive rate of SARS-CoV-2 starts at a similar level for all clusters, the speed at which counties in each cluster respond to the disease, and reduce its R_t differs. The reproductive rate in metropolitan counties (cluster 4), which tends to have higher reliance on public transportation decreases over the entire period. This is especially apparent going from March 15 to March 25. On the other hand, cluster with little transportation use, such as cluster 5, was estimated to have quickly reduced transmission rates. This suggests that if people have a higher reliance on public infrastructure, more stringent interventions are necessary to reduce transmission. The staggered drops in reproductive rates are attributed to both varied efficacy of NPIs (to be discussed in the next section) and counties implementing NPIs at different times. Comparisons with more features are shown in Section 4.

Estimates of Initial and Current Reproductive Rates and Number of Infected

From Table 1, we observe that most counties exhibited an initial reproductive rate R₀ above 3. As of May 28^th, however, most have successfully reduced the reproductive rate to R_t ≈ 1 after implementing NPIs. It is estimated that 18.4 to 42.7% of the population (95% confidence interval) has been infected in New York, NY, which has the highest number of cases. Based on the initial reproductive rates between 2 and 4, herd immunity is reached only after 50 – 70% of the population has recovered,^8,9 suggesting that all U. S. counties are far from resilience. Consequently, easing restrictions is likely to result in subsequent waves of the epidemic. We state these findings for 15 heterogeneous counties in Table 1 and provide the same metrics for all 1,417 counties online at https://github.com/JieYingWu/npi-model, where we also provide code and data required for reproduction.

Learned Effects of NPIs

Using a cluster-specialized model, we quantify the effectiveness of NPIs, as shown below in Table 2 and Table 3.

We note that our model estimates different behavior across counties. While all counties have implemented a similar set of interventions, their estimated effects are substantially different in each respective cluster. For example, metropolitan counties (cluster 4) were estimated to have a strong response to stay-at-home orders, while rural and suburban areas (clusters 1, 2, 3, and 5) responded more to national-level interventions according to our model.

One caveat is that the effects of interventions that were implemented close together are difficult to disentangle. Many local governments implemented formal NPIs in immediate response to the federal guidelines, leading a quick succession of NPIs coming into effect. Disentanglement of the effects of NPIs is explored more in Section 4.

Another limitation of our model is that the federal level interventions may have come before some counties have seen any cases. This introduces ambiguity between R₀ estimates and the weight attributed to federal NPIs, i.e., federal guidelines and travel ban. Similar ambiguity is less prominent for other interventions, since they are not generally implemented before a county has seen cases. Additionally, since the epidemiological model used here is mechanistic, it only allows for changes of the transmission rate at the time of implementation of interventions. Other events, such as high-profile cases and cancellations of prominent festivals, likely contributed to increasing awareness of the disease and may also have effects on individual behavior, and therefore the R_t. These effects cannot be attributed to a specific date and may thus affect the weights of interventions that come into effect at around the same time.

Method Overview

Figure 3 shows an overview of the proposed method. We train a Bayesian mechanistic model (BMM) model in a cluster-specialized manner and compare with the baseline model trained on all counties. We show the stability of our model in two ways, first by withholding random days during training and second by withholding given counties. We compare the predictions based on incomplete data with those based on complete data.

• Download figure
• Open in new tab

Figure 3:

The fitting process for our model. We compare the performance of the baseline model (a), which is fit to all eligible counties, with cluster-specialized models BHI 1–5. (b) BMM 1-5. (c) Two validation methods for our model. Validation A tests the stability of the model by withholding 3 random days and comparing the outputs. Validation B withholds a set of counties for comparison rather than days.

County-level Clustering

Quantifying changes in COVID-19’s reproductive rate is complicated when considering differences at the county – rather than the national – level. Parametric epidemiological models are optimized to describe fatality counts, the volume and reliability of which decreases outside of highly populated regions. To make up for this scarcity, we leverage a balanced clustering of U.S. counties to aggregate data from similar counties in the same state, treating them as a single “super-county,” if they have identical NPI implementation dates. This has the advantage of considering counties that would otherwise be excluded without assuming that the spread of the disease in those counties follows the same trend of more advanced regions in the same state or country. Figure 4 shows how for a given state-in this case Texas, and a given cluster, the counties with 1–49 cumulative deaths from COVID-19 are treated as a single entity. In addition to this data aggregation strategy, we fit a cluster-specialized model to each set of counties, quantifying the possibly disparate effects of the NPIs in each type of county, as detailed below.

• Download figure
• Open in new tab

Figure 4:

(a) The total confirmed deaths caused by COVID-19 for counties in Texas, as of May 28^th. (b) Cluster labels for each Texas county, based on demographics, education, density, and other factors. (c) Texas counties in cluster 1 having 1–49 cumulative deaths as of May 28^th, 2020 and the same NPI implementation dates.. To enable robust epidemiological models, these counties are treated as a single “super-county.”

To generate this clustering, we partition 3,059 U.S. counties into five clusters based on variables which directly affect disease spread, using a Gaussian mixture model.¹⁰ Table 4 summarizes these variables, which include demographic, economic, and public transit capacities that we have gathered, processed for machine readability, and released in a publicly available dataset.¹¹ Sources include the United States Census Bureau, the United States Department of Agriculture Economic Research Service and the Center for Neighborhood Technology. A full list of sources can be found on the corresponding website.¹¹ To incorporate potential exposure, we consider county population density, housing density, and land area. Additionally, we consider portions of the population for age- and gender-based demographic categories, due to COVID-19’s disparate effects on these groups.^12–16

Our clustering is also based on socioeconomic variables, which may indicate behavioral traits relevant to the spread of COVID-19. For instance, workers with tertiary education are more likely to hold office-type jobs which can be done from home.¹⁷ At the same time, many secondary-education jobs have been deemed essential, resulting in a high contact rate, which in turn increases the likelihood of infection. Thus, our clustering considers college education, poverty, unemployment, and median household income for each county as a reflection of the overall job composition in the local area. Finally, we include a population-weighted transit score due to the likelihood of transmission in the enclosed, possibly crowded space that public transport entails.

Figure 5 and Figure 6 show the plotted over median income and density for counties and super-counties in the different clusters. Super-counties are visualized as a single point using their population-weighted average for that feature. The plots show the distribution of the cluster over the features and its correlation with how R_t changes.

• Download figure
• Open in new tab

Figure 5:

Scatterplot and density distribution plot for counties and super-counties comparing R_t over time to median household income. Colors indicate which cluster the county or super-county belongs in, as indicated by Figure 4.

• Download figure
• Open in new tab

Figure 6:

Scatterplot and density distribution plot for counties and super-counties comparing R_t over time to population density. Colors indicate which cluster the county or super-county belongs in, as indicated by Figure 4.

Notably, we exclude ethnic demographics from the variables considered during clustering. because we assume that no direct relationship exists between race and incidence of COVID-19. Rather, such a relationship would be mediated by socioeconomic factors that are already included for similarity assessment.

Modelling the Effects of NPIs

Validation

We propose three validation schemes to validate our approach. First, we show that our model is stable, so a small change in the input produces a small change in the output. Second, we separate our counties and super-counties into train and test sets. We fit models for each cluster as well as for all counties and super-counties on either, both the train and test set or with only the train set. We evaluate using the parameters of the “correct” cluster model to predict fatalities for the held-out regions. Comparing these predictions to those made when the regions are included during training confirms that the model is not over-reliant on each data point. Third, we evaluate how well our model discerns the effects of individual NPIs and discuss biases it may have to attribute more weight to certain NPIs.

Validating the Stability of the Model

To show stability in our model, we train the model, withholding three non-consecutive days chosen at random, and compare the predictions for these days to the baseline model fit on full data. Only days with non-zero death counts are selected as potential leave-out candidates. We set the threshold for data selection at counties with 50 cumulative deaths by May 28^th, which results in 211 counties. (For model validation, we do not use super-counties.) We fit the model over 300 total iterations, with 150 of those as warmup iterations. From the uncertainty distributions obtained as results from the models, we use the expected values of the deaths of the held-out days to compare the two models and report our observations in Figure 7a. The clustering of points around the optimal fit visually confirms that the expected values of the baseline model fit those from the model with held-out days. In Figure 7b, we compare the average county-wise expected value of deaths of the withheld days and observe a more concentrated distribution around the line of optimal fit.

Using the Pearson-correlation coefficient to quantify the similarity between the expected values, we find r-values greater than 0.99 for the held-out days and the average held-out days which suggest that the model is robust against small input perturbations.

• Download figure
• Open in new tab

Figure 7:

The predicted values of the validation and the baseline model. (a) Predictions for withheld days, as shown by the optimal fit, are remarkably close to the true values. (b) The three held-out days per county are averaged and compared.

Validating the Advantages of Clustering

We show the advantage of cluster-specialized models for quantifying the effects of NPIs by comparing their performance with models trained on different clusters or on the national level. While a clustering of U.S. counties may or may not be interesting, its value related to COVID-19 comes from its ability to identify epidemiologically meaningful differences among various regions in the country. When fitting each model, we withhold a validation region from each cluster and use the remaining counties to fit Equation(1), for both cluster-specialized and baseline models. We then use the learned weights α_i to initialize a fixed-α BHI model for each withheld region.

With this scheme, we show the advantage of clustering U.S. counties in three ways. First, we show that within the same cluster, we obtain comparable predictions for a county whether or not it is included in training. Second, we observe that α values learned from a different cluster produce substantially different predictions. Finally, we apply the α values learned at the national level to similar effect. This demonstrates how cluster-specialized models can reveal trends in the spread of COVID-19 that are not apparent under the assumption that NPIs have universal effect at the national or state level.

Figure 8 shows the first and second validation strategy for cluster 3, which consists of less-populous counties with large land area. For such regions, which necessarily have fewer cases to support local models, it is vital to gain accurate insight from the entire training pool and verify these insights through the aforementioned validation process. Thus, in this illustrative example, we use data up to May 18^th, at which point the continuing spread of cases was not immediately clear. Indeed, Figure 8a shows the predictions obtained for seven counties in Arizona (treated as a single super-county). One can readily observe the model takes into account random outliers while following a general upward trend. Figure 8b exhibits less awareness of outliers but still follows the same trend, despite using α values obtained from training on cluster 3 with these same counties withheld. On the other hand, Figure 8c shows the predictions made when using the wrong α values, from cluster 4, which consists of densely populated city centers. As can be seen, this results in a markedly different prediction for the trend of the disease than actually occurred up to May 28^th in Figure 8b. This is a critical difference, which could result in very different decisions for implementing or removing NPIs.

• Download figure
• Open in new tab

Figure 8:

The fatalities predicted by our model for seven counties in Arizona (treated as a single super-county), with fatalities data from Apache, Cochise, La Paz, Mohave, Navajo, Yavapai, and Yuma counties, using data up to May 18^th. (a) When the super-county is included in training, the cluster-specialized BHI model incorporates outlier events in deaths from COVID-19 while still describing the overall trend. (b) When the super-county is withheld from training but uses α values from the same cluster, the model still captures the same overall trend of increasing deaths. (c) On the other hand, using values from a different cluster makes it appear as though the curve has flattened, even though (d) shows the continued upward trend based on data up to May 28^th.

Finally, we compare the performance of our cluster-specialized model with the baseline model, which uses all eligible counties for training. In many cases, this difference is slight, and both the baseline and cluster-specialized models exhibit similar trends for a given region. However, for certain areas the difference can be profound. Take, for example, the District of Columbia, which has become a significant hot-spot for COVID-19. The baseline model, shown in Figure 9a, has only recently flattened, whereas the cluster-specialized model in Figure 9b reflects the reality that efforts to combat the disease have had much greater effect, significantly reducing the number of deaths. This is because when forced to accommodate counties across the U.S., the baseline model emphasizes federal guidelines with a mean α value of 0.804 (see Table 3). On the other hand, our cluster-specialized model found that stay-at-home orders were much more effective for counties in cluster 4, with a mean α value of 0.963. This difference illustrates the advantage of specializing a BHI model based county-level characteristics; it allows the model to include a greater number of counties within separate clusters while not being forced to overgeneralize and, in doing so, compromise some certainty.

• Download figure
• Open in new tab

Figure 9:

Predicted fatalities in Washington D.C. up to May 28^th, 2020. (a) The fatalities predicted by the baseline model, which observes all eligible counties during training. (b) The fatalities predicted by a cluster-specialized model, showing a very different trend in the course of the spread of COVID-19. Although this difference is more pronounced for the District of Columbia than for most regions, it underscores the potential for misleading predictions when the United States’ heterogeneity is not taken into account.

Disentanglement

One drawback of thes model is that it cannot disentangle implementations that came into effect at the same time. For example, states often closed public schools at the same time as federal guidelines were issued so it is difficult to discern the individual effect on reducing R_t. Additionally, in Table 5, we observe that banning gatherings of 50 or more people often occurs at the same time as banning gatherings 500 or more, and restaurants and entertainment venues are often closed together. This suggests that these pairs of interventions may not or be poorly disentangled.

To further investigate the model’s ability to disentangle intervention weights, we create simulated trajectories of counties’ deaths and cases counts based on their R₀ and the dates on which the interventions came into effect. Using all counties that have more than 50 cumulative deaths on May 28^th without super-counties, we seed each county with 200 cases in each of the first 6 days. To simulate county-specific trajectories, we construct a sets of generated timeseries. In the first set, we assign intervention weights α_i to be randomly generated from a Gamma distribution, the same distribution as our prior on the Bayesian mechanistic model adjusted to be in the range of our learned weights. We then calculate what the R_t on each day must have been based on the R₀ and the interventions in place. Once we have the seeded infection and the R_t trajectory for each county, we can calculate daily infections and thus expected fatalities. Using the simulated trajectories, we fit the model. Table 6 compares the weights used for generation with the weights that the model learned.

We observe that the effects of individual NPIs are not well disentangled in general. The model tends to attribute more weight to few NPIs rather than spread out the weight evenly. Specifically, the model tends to put more weight on stay-at-home orders. This may be because interventions I2-I8 are often implemented close together (see Table 5) and it is difficult to attribute effect to any single one of them on a national scale. While we can conclude that the trajectories the model predicts are reliable, due to their match to measured death, and therefore the overall decrease in R_t is reliable, attributing decreases to any individual NPI is challenging whenever the difference in implementation date is small. This observation seems to be in line with the similarly large confidence intervals reported in previous work on varied models and regions.^2–7