key: cord-0602965-zhumf1mt
authors: Giffin, Andrew; Gong, Wenlong; Majumder, Suman; Rappold, Ana G.; Reich, Brian J.; Yang, Shu
title: Estimating intervention effects on infectious disease control: the effect of community mobility reduction on Coronavirus spread
date: 2021-03-07
journal: nan
DOI: nan
sha: 50754bb8062d3556cfbb1d09ff70fb259e223ee7
doc_id: 602965
cord_uid: zhumf1mt

Understanding the effects of interventions, such as restrictions on community and large group gatherings, is critical to controlling the spread of COVID-19. Susceptible-Infectious-Recovered (SIR) models are traditionally used to forecast the infection rates but do not provide insights into the causal effects of interventions. We propose a spatiotemporal model that estimates the causal effect of changes in community mobility (intervention) on infection rates. Using an approximation to the SIR model and incorporating spatiotemporal dependence, the proposed model estimates a direct and indirect (spillover) effect of intervention. Under an interference and treatment ignorability assumption, this model is able to estimate causal intervention effects, and additionally allows for spatial interference between locations. Reductions in community mobility were measured by cell phone movement data. The results suggest that the reductions in mobility decrease Coronavirus cases 4 to 7 weeks after the intervention.

Integrating causal inference into this spatial framework presents some challenges. Existing SIR models are generally used for prediction and forecasting, but not for ascertaining causal effects of interventions. Causal methods generally require adjusting for confounders that are related to both intervention and response, which allows for identification of the causal intervention effect as opposed to the observed association between intervention and response. However, if interventions at one location affect the response at other locations -a phenomenon known as interference -then the standard causal framework which adjusts only for local confounding is rendered ineffective. With a virus that spreads across space, clearly some degree of interference is involved, and a causal analysis must take this into account. This is typically done by placing assumptions on the form of interference. Commonly used assumptions on interference (for a recent review, see Reich et al., 2020) include "partial interference" in which the population is divided into isolated clusters (Halloran and Struchiner, 1991; Sobel, 2006) , as well as "network interference" in which interference is allowed along a pre-specified network, and "spatial interference" in which interference is tempered by the distance between locations . This paper uses a parsimonious assumption on interference, which allows for interference from neighboring counties only. This can be thought of as a simple version of both the spatial and network interference assumptions. Intervention effects are then divided between direct (within county) interventions and indirect (between county) interventions.

The framework that we propose approximates the standard SIR model, in a way that allows for estimation of the causal effects of interventions. In particular, we use the fact that infection rates are likely smooth across space, to obtain an approximation to the rate of infection at across space and time. This approximation is a function of the intervention variables, covariates, and a spatiotemporal model, and allows us to fit the observed data and obtain intervention effect estimates.

The remainder of the paper proceeds as follows. The motivating data are described in Section 2. Section 3 introduces the statistical method and corresponding theoretical properties. The method is evaluated using a simulation study in Section 4 where we show that the approximate model can recover causal effects of data generated from the SIR model. Section 5 applies the proposed model to estimate the effect of decreases in mobility on Coronavirus cases, finding that decreases in mobility cause a decrease in local Coronavirus cases 4-7 week later. Section 6 concludes.

We estimate the effect of mobility on the number of observed Coronavirus cases. There are three primary types of data for analysis: intervention (mobility reduction), response (Coronavirus cases) and covariates. All variables are defined temporally by the week (t) and spatially by the county (j). Data from March 6, 2020 through October 8, 2020 for 3,137 counties or county equivalents are used in this study.

The response data Y j (t) are new, recorded cases of Coronavirus in a given county/week. These data are taken from the publicly available Johns Hopkins University Coronavirus Resource Center (Dong et al., 2020) , which aggregates Coronavirus case counts and provides a daily, cumulative count of cases in each county.

The intervention data are publicly available measures of mobility -as measured by Google devices and provided by Google LLC (2020) -which have been shown to have strong associations with Coronavirus case data (Chang et al., 2021; Yilmazkuday, 2020) . We use an aggregate mobility measure that includes the categories: "retail and recreation", "grocery and pharmacy", "transit station", "workplace", and "residential" mobility. These measure both the number of visits to such locations as well as the length of stay in comparison to baseline. These data are given as percentage change from a baseline level which was taken over January 3, 2020-February 6, 2020. Intuitively, in the months since the pandemic began, "retail and recreation", "grocery and pharmacy", "transit station" and "workplace" mobility have decreased substantially from the baseline and "residential" mobility (i.e., amount of time spent at home) has increased. The specific metric that we will use as intervention A j (t) for a given county/week is the mean of percentage change from baseline of "retail and recreation", "grocery and pharmacy", "transit station", "workplace", and negative "residential" mobility.

Because of privacy concerns, county/days with too few users are not provided by Google. When a subset of the categories are provided for a given county/day, the mean of the available categories is taken. When none of the categories are provided for a given count/day, we impute the missing value from any available first and second degree neighbors (with values weighted by 1/distance, where weights are summed to one).

In addition to intervention and response we include a number of important covariates. From the 2016 American Community Survey we have poverty rate, population density, median household income, and total population, and from the U.S. Census we have median age and the number of foreign born residents (U.S. Census Bureau, 2016) . The Bureau of Labor Statistics provides the number of people employed in healthcare and social services, classified as North American Industry Classification System, Sector 62) (U.S. Bureau of Labor Statistics, 2019). The Environmental Systems Research Institute (ESRI) provides publicly available data on the number of hospital beds and the number of ICU beds (ESRI, 2020). Lastly, because Wu et al. (2020) shows that PM 2.5 (particulate matter smaller than 2.5 micrometers) may interact significantly with COVID-19, we include county level PM 2.5 estimates from 2016 (Alexeeff et al., 2015; Chudnovsky et al., 2012) . We also include non-static county-level daily meteorological data, as this could inform levels of mobility and spread of disease. We use daily mean, maximum and minimum temperature (degrees Celsius), relative humidity and mean dew point (degrees Celsius) from NOAA's Global Surface Summary of the Day using the "GSODR" R package (Sparks et al., 2017) . Then we obtain county-level data by predicting at the mean county latitude and longitude coordinates from the 2010 Census using thin plate spline regression from the "fields" R package (Nychka et al., 2014) .

Let Y j (t) be the number of new cases of Coronavirus reported during week t ∈ {1, ..., T } and county j ∈ {1, ..., J}, where T = 31 and J = 3,137. For county j, denote N j as the population size and X j (t) as a vector of covariates for week t. The direct intervention variable, A j (t), is the mobility percentage change from baseline for county j and week t. In addition to the direct intervention, we allow for an indirect (spillover) interventionà j (t). This variable captures the interventions received by neighboring counties. Including this terms allows for interference between neighbors, since the response at county j can now be impacted by the interventions in surrounding counties. The spatial configuration of the counties is summarized by their adjacency, with c jj = 0, c jk = 1 if counties j and k are adjacent and c jk = 0 otherwise.à j (t) is defined as the mean direct intervention over the m j = J k=1 c jk adjacent regions:à j (t) = J k=1 c jk A j (t)/m j . Similarly, the mean of the adjacent covariate values isX j (t) = J k=1 c jk X j (t)/m j .

To motivate the proposed statistical framework, we begin by describing the discrete-time, spatial SIR model proposed in Paeng and Lee (2017) . Let S j (t), I j (t) and R j (t) be the number of susceptible, infected and recovered individuals in county j at time t. These three states evolve over time but always satisfy the constraint that N j = S j (t) + I j (t) + R j (t).

For each county j the state evolution is similar to the classical SIR model (Kermack and McKendrick, 1927) :

where γ j > 0 is the recovery rate and λ j (t) = β

is the rate new infections in county j and week t. Here W jk is proportional to the contact rate between individuals in region j with individuals in region k.

We expand on this by allowing the infection rate β j (t) > 0 to vary with region and time:

This allows us to model spatiotemporal variation in β j (t) as a function of the covariates X j (t) and direct and indirect intervention variables A j (t) andà j (t). In the absence of other information on connectivity, we simply set W jj = (1 − φ) and W jk = φ · c jk /m j for φ ∈ [0, 1] so that large φ leads to strong spatial dependence and vice versa. A major difficulty in using model (1) is that we do not observe the latent states S j (t), I j (t) or R j (t). We learn about these latent processes via the reported number of new cases Y j (t). Further complicating the statistical model is that there may be under-reporting and a lag time between an increase in the true and reported infection rate due to the latency of the disease. Because of these issues, we link Y j (t) and λ j (t) in the SIR model by the over-dispersed Poisson distribution

where p ∈ (0, 1) accounts for under-reporting, g j (t) iid ∼ Normal(0, τ 2 ) captures over-dispersion, and l ∈ {0, 1, 2, ...} is the reporting lag.

The SIR model in Section 3.2 is an elegant way to model and forecast the spread of a disease through a population. However, the solution to the difference equations is complex, which hinders the ability to estimate model parameters in a statistical manner. Previously methods have been proposed to mimic the mechanistic model to provide realistic forecasts, e.g., Buckingham-Jeffery et al. (2018) . Here, we propose an approximation to the SIR model to allow for estimation of intervention effects on local infection rates.

The key insight is that the rate of new infections in (2) can be re-expressed as

Both θ j (t) and v j (t) in (4) are latent spatial and temporal processes due to the unobserved S j (t) and I j (t), which we approximate by spatiotemporal models. We model θ j (t) as a separable spatiotemporal continuous autoregressive model (Stein, 2005) . Denote θ(t) = {θ 1 (t), ..., θ J (t)} and θ = {θ(1) , ..., θ(T ) } . Then the spatiotemporal conditional autoregressive model (STCAR) for θ is multivariate normal with mean zero and covariance

where M is diagonal with diagonal elements m j and the (j, k) element of C is 1 if sites j and k are adjacent and zero otherwise.

, so that ρ s determines the strength of spatial dependence and σ is the scale parameter. We denote this model as θ(t) ∼ CAR(σ, ρ s ). Similarly, temporal dependence is governed by Ω(ρ t ), which has the same form as Σ(ρ s ) except for temporal adjacency structure with time t having t − 1 and t + 1 considered neighbors. We refer to this as the STCAR(σ, ρ s , ρ t ) model.

The second term in (4), v t (t), sums over the m j regions adjacent to region j, and is a function of the local differences, I k (t) − I j (t). Assuming the number of infected individuals varies smoothly across space, these local differences should be small. Therefore, one approximation we consider is simply setting v j (t) = 0 for all j and t. For cases where these terms cannot be removed, we note that because they are functions of local differences they should be less spatially correlated than the spatial process itself, and thus a second approximation is log{v

. A crucial point is that neither our model for θ t (s) nor v t (s) attempt to retain the mechanistic properties of I j (t) and S j (t). Therefore, this model would not provide reliable forecasts about future disease spread. However, forecasting is not our objective. Our aim to is provide a computationally feasible method to use observed data to estimate effects of intervention variables on local infection rates. The following section verifies that the approximation indeed has a causal interpretation under standard assumptions from causal inference, and proposes an estimation procedure for these effects.

For the final component, modeling the spatiotemporal infection rate β j (t), we regress log{β j (t)} on covariates X j (t) andX j (t) as well as mobility variables A j (t) andà j (t) as

where δ 1 and δ 2 quantify the direct and indirect (spillover) causal effects of mobility on infection rate. The covariate vectors X j (t) andX j (t) can include both the original covariates or covariates derived as propensity scores (Section 3.4).

Since it is impossible to identify the reporting rate p in (3) and intercept terms α 0 and µ v we fit the final model

The final term in (7) combines the overdispersion term g j (t) and the nugget term v j (t). When the dimensions of the covariates are large, we propose dimension reduction for log{β j (t)} using the generalized propensity score (Imbens, 2000) in Section 3.4. To complete the Bayesian model, we specify noninformative independent Normal(0, 10 2 ) priors for α 0 , δ 1 , δ 2 , µ v and the elements of α 1 and α 2 , and noninformative priors for covariance parameters σ 2 , τ 2 , σ 2 v ∼ InvGamma(0.1, 0.1) and ρ s , ρ t ∼ Uniform(0, 1).

In this section, we provide a causal interpretation of δ 1 and δ 2 under the potential outcomes framework (Rubin, 1974) . To ease discussion, let Y j (t) be the total number of cases, rather than reported number (i.e., p = 1). We use the potential outcomes framework to define the causal effect of mobility on infection rate. We use overbar to denote all history. Let a j (t) = (a j (1), . . . , a j (t)) be the trajectory of mobility level at region j through time t. Let a(t) = (ā 1 (t), . . . ,ā J (t)) be the trajectory of mobility level for all regions through time t. Define Yā Assume the potential outcome model for the SIR model with θ (5) and (6), and

Consider two regimesā(t) andā (t), where all components are the same except for a j (t) = a j (t) + 1, and v j (t) = 0. Then under model (9),

whereX(t) is the full history of covariates and intervention through time t, andS(t − 1) andĪ(t − 1) are the full history of susceptible and infected,

Then exp(δ 1 ) entails the risk ratio of new cases by increasing mobility level by 1 unit locally (direct effect). If v j (t) is nonzero but small in comparison to βā Similarly, if we consider two regimesā(t) andã (t) where all components are the same except forā j (t) =ā j (t) + 1, and v j (t) = 0, then

encodes the risk ratio of new cases by increasing mobility level by 1 unit non-locally (indirect effect). For the effect from neighbor k, consider two regimesā(t) andā (t), where all components are the same, except for a k (t) and a k (t) + 1, and v j (t) = 0. Then

which encodes the indirect effect of area k on area j.

To link the potential variables and the observed variables and identify the direct and indirect effects, we invoke the following assumptions.

Assumption 1 (Interference) Given the past information through t, the intervention A(t) affects the potential infection rate at time t and region j only through A j (t) andà j (t).

To facilitate recovery of causal estimates, we define a generalized propensity score for the direct/indirect effects

In practice, we often assume that the direct and indirect interventions are independent, and can write this score as the two separate components

denoted by e jt {κ;X(t),S(t − 1),Ī(t − 1)} andẽ jt {κ;X(t),S(t − 1),Ī(t − 1)}, respectively.

Lastly we require a positivity assumption that ensures all regions can take any plausible mobility level. This is formulated based on the propensity score:

Assumption 4 (Positivity) For allX(t),S(t−1),Ī(t−1) with f {X(t),S(t−1),Ī(t−1)} > 0,ē jt {κ 1 , κ 2 ;X(t),S(t − 1),Ī(t − 1)} > 0, ∀t, κ 1 , κ 2 where f (·) denotes the generic probability function; i.e., it is a probability density function for a continuous variable and a probability mass function for a discrete variable.

Under Assumptions 2 and 3, the induced model from (9) for Y j (t) is Y j (t) ∼ Poisson{λ j (t)}. One can fit the induced model with the observed data to estimate the causal parameters. BecauseX,Ī, andS are high-dimensional, fitting the above model directly may become a daunting task. We consider the generalized propensity score approach for dimension reduction. Under Assumption 3, and following Rosenbaum and Rubin (1983) or Imbens' generalized propensity score approach, we can show that

In practice, following Giffin et al. (2020) we model the two generalized propensity scores using linear regressions of A j (t) andà j (t) onto covariates and past interventions and responses. Then, their distributions for a given j and t (i.e., their generalized propensity scores values) can be completely summarized with a single sufficient statistic: the estimated conditional mean (i.e., the fitted value from the regression). For this reason, henceforth e j (t) andẽ j (t) (without the index for κ) will refer to these sufficient conditional means, for which the result (Balancing) applies. This shows that the generalized propensity score serves as a balancing score; i.e., at the same level of the propensity score, the distribution of the history of confounding variables are the same across different intervention levels.

Allowing X andX to include propensity score components, the induced model for Y j (t) given the intervention and generalized propensity score is

The derivation is provided in the Appendix. Therefore, we can fit the model

To model the propensity of A t given the past history, we require some structural assumption to reduce the dimension of the historical variables. We posit a model for A j (t) adjusting for its direct causes A j (t − 1), A j (t − 2), X j (t), X j (t − 1) and Y j (t − 1) to get e j (t).à j (t) can be modeled similarly with the neighbor-averaged variablesà j (t − 2),X j (t),X j (t − 1) andỸ j (t − 1). In practice, simply taking the average over neighbors of e j (t) gives similar estimates forẽ j (t).

In this section we conduct a simulation study to evaluate the statistical properties of the causal effects that stem from the approximate SIR models. Data are generated from the full SIR model in Section 3.2 and analyzed using the approximated spatial model in Section 3.3. The objectives are to determine when the approximations give unbiased point estimation and valid interval estimation, and to compare different approximations using these criteria.

We generate data on a 15 × 15 square grid of J = 225 regions for T = 30 time steps. Rook neighbors are considered adjacent. Each region has population N j = 100, 000 and initial states are generated as I j (1) = 100 exp{U j }, R j (1) = 0 and S j (1) = N j − I j (1), where (U 1 , ..., U J ) ∼ CAR(1, ρ s ). A confounding variable X j (t) is generated from the STCAR(1, ρ s , ρ t ) distribution and the intervention variable is generated as A j (t) = ρ x X j (t)+ 1 − ρ 2 x E j (t) where E j (t) is also generated from the STCAR(1, ρ s , ρ t ) distribution. The latent states for times t ∈ {2, ..., T } are generated from the full mechanistic model given by (1) and (2) with recovery rate γ = 0.1 and infection rate log{β j (t)} = α 0 + X j (t)α 1 +X j (t)α 2 + A j (t)δ 1 +Ã j (t)δ 2 for α 0 = −3, α 1 = 0.5, α 2 = 0.3, δ 1 = 0.5 and δ 2 = 0.2. The data are then generated as Y j (t) ∼ Poisson{pλ j (t − l)} where the reporting rate is p = 0.5 and the lag is l = 2.

We compare six simulation scenarios defined by the spatial (ρ s ) and temporal (ρ t ) dependence of the intervention variable A j (t), the strength of confounding (ρ x ) and the strength of spatial dependence (φ) in the SIR model: 1. Base model: ρ s = 0.9, ρ t = 0.5, ρ x = 0.5 and φ = 0.4 2. Strong spatial dependence in A j (t): ρ s = 0.99, ρ t = 0.5, ρ x = 0.5 and φ = 0.4 3. Strong temporal dependence in A j (t): ρ s = 0.3, ρ t = 0.9, ρ x = 0.5 and φ = 0.4 4. Strong spatiotemporal dependence A j (t): ρ s = 0.9, ρ t = 0.9, ρ x = 0.5 and φ = 0.4 5. Strong confounding: ρ s = 0.9, ρ t = 0.5, ρ x = 0.9 and φ = 0.4 6. Weak SIR spatial dependence: ρ s = 0.9, ρ t = 0.5, ρ x = 0.5 and φ = 0.2 For each scenario we generate 100 datasets.

The propensity score e j (t) is the fitted value from a least squares regression of A j (t) onto A j (t − 1), A j (t − 2), X j (t), X j (t − 1) and Y j (t − 1). We then model log{β j (t)} as a linear combination of A j (t),Ã j (t), X j (t),X j (t), e j (t),ẽ j (t), e j (t) 2 ,ẽ j (t) 2 and e j (t)ẽ j (t). For each simulated dataset we fit the full model in (7) with the propensity scores included as covariates ("Full"), and compare this model to the full model but with v j (t) = 0 ("No nugget"), the full model but without the propensity scores e j (t) =ẽ j (t) = 0 ("No PS") and the full model but without spatial dependence ρ s = 0 ("Non-spatial"). For the non-spatial model we also set v j (t) = 0 because the MCMC algorithm did not converge with these terms included. The models are fit using responses from time points t ∈ {5, ..., n t } to accommodate lagged predictors. We fit the model using MCMC with 10,000 iterations and the first 2,000 samples are discarded as burn-in to approximate the posterior median and 95% interval of the direct (δ 1 ) and indirect (δ 2 ) effects. Table 1 . Spatiotemporal model simulation study results. The six data-generating scenarios are given in Section 4.1. The two recommended models, the full model ("Full") and the model without a nugget term ("No nugget"), are given in bold; and are compared to the model without the generalized propensity score ("No PS") and the model without spatial dependence ("Non-spatial") in terms of bias and empirical coverage of 90% prediction intervals, separately for the direct (δ 1 ) and indirect (δ 2 ) effects. Bias is multiplied by 100 and standard errors are given in subscripts. 

The results are given in Table 1 . Broadly speaking the Full model and the No Nugget model perform well. Giving similar precision and coverage results under all scenarios, they generally outperform the other models in terms of coverage and bias. Moreover, their coverage rates remain reasonably close to their 90% target, although they appear to have a small, positive bias. The Non-spatial model performs worse than the Full and No Nugget models in both precision and coverage, but it remains competitive. The No PS model performs the worst, with extremely poor performance for the strong temporal/spatiotemporal dependence scenarios. The four models seem to have the most difficulty estimating effects when there is strong temporal/spatiotemporal dependence in A j (t) or strong confounding. All models appear to estimate the direct effects better than the indirect effects.

We implement similar models that were used in the simulation study on the real data: A Full model, a No Nugget model, a No PS model, and a Non-spatial model. The models are fit using responses from time points t ∈ {8, ..., 31} (April 24, 2020 through October 8, 2020) to accommodate lagged predictors. The direct propensity score is estimated with least squares regressions of A j (t) onto the previous intervention A j (t−1), the local covariates X j (t) and X j (t − 1), the number of weeks since the first case in county j (enters both linearly and quadratically), time t (enters both linearly and quadratically), a time and intervention interaction t · A j (t − 1), the baseline level of mobility A j (1), and log{Y j (t − 1)} − log(N j ). The propensity score e j (t) is set to the fitted values of the resulting model. The propensity of the spillover intervention is estimated as the local average of its neighboring direct scores: e j (t) = k∼j e k (t)/m j . The direct and indirect propensity scores are added to X and X, respectively, as in (8), along with l-lagged interventions and covariates. Each model is implemented with lags of l = 0, . . . , 7. The No Nugget, No PS, and Non-spatial models were run for 100,000 MCMC iterations, with a burn-in tuning period of 20,000 iterations. The Full model required more iterations to converge, so these values were increased to 200,000 and 40,000, respectively. Table 2 gives the estimated posteriors for the direct and spillover effects for the four different models across lags l = 0, . . . , 7, and Figure 1 shows these results graphically for the Full model. Rather than present the estimated values for δ 1 and δ 2 , the estimates shown are 100 · {exp(50 · δ i ) − 1}, i = 1, 2, and correspond to the expected percentage increase in cases from an increase in mobility 50% above baseline. Because it is not clear that the observed data can identify the correct lag, we provide model results for all lags, and speculate that lags of 3-5 weeks are most appropriate to see effects of changes in mobility. The Centers for Disease Control and Prevention guidelines suggest that COVID symptoms appear 2-14 days after virus exposure (Centers for Disease Control and Prevention, 2020) . Given this, 3-5 weeks would allow the virus to spread through roughly several generations of peoplewhich would be enough to see changes in county-level case counts.

The key trend evident in these results is that the direct effect appears to be positive and significantly different from zero for higher lags. The Full model gives positive significant Table 2 . Posterior estimates for direct and indirect treatment effects for lags 0-7 weeks for all models. Estimates correspond to the expected percentage increase in cases from an increase in mobility 50% above baseline (100 · (exp(50 · δ) − 1)). 95% confidence intervals are shown in parentheses, and estimates significant at the 95% level are denoted with * .

Lag ( -7.5, 66.5) No Nugget 0 -0.8 (-11.6, 11.7) -21.7 (-41.1, 4.1) 1 -1.5 (-11.8, 9.8) -28.5 (-46.4, -4 .6)* 2 3.6 (-7. 24.2 (6.8, 44.5) * 6.1 (-11.1, 26.5) estimates on δ 1 for lags 4-7 weeks; the No Nugget for lags 5-7 weeks; the No PS model for lags 2-7 weeks; and the Non-spatial model for lags 5-7. The size of this effect differs between models/lags, but the Full model for lags 5-7 predicts that a mobility level 50 percent above baseline will increase Coronavirus cases by between 11 and 41 percent. This trend is demonstrated in the grouping of posteriors above zero in the top sub-figure of Figure 1 . This provides evidence that a decrease in local mobility should decrease local Coronavirus cases roughly 4-7 weeks after the change in mobility. The indirect effect estimates are more variable with wider credible intervals, as illustrated in Figure 1 . For the Full model, a lag of 7 gives a positive but non-significant estimate for the indirect effect, although the indirect effect posteriors for other lags seem centered at zero, with the exception of lag 1 which is significant and negative. For lag 1 the No Nugget model also shows this counter-intuitive result. The No PS model shows a number of positive and significant indirect effects, and the Non-spatial model paradoxically shows several significant indirect effects -some positive and some negative. It is not entirely clear why there is more variability in the indirect effects or why we see these counter-intuitive results. Some of these could simply be spurious significance resulting from estimating many lags. Another possibility is that our measure of connectivity does not precisely measure the levels of contact between different counties. This could explain the wide intervals, as the connectivity is overestimated between some counties and under-estimated between others. In any event, the indirect effect story appears less clear cut than that of the direct effect. Table 3 gives parameter estimates from the No Nugget/Lag 5 model. In addition to the intervention effects and accounting for the propensity scores, the covariate coefficients tell an interesting story about their association with Coronavirus case counts. For example, expected Coronavirus case counts are decreasing in median age of population, population density, county population, and relative humidity. Alternatively, Coronavirus case counts tend to increase with the poverty rate, PM 2.5 , number of hospital beds, number of ICU beds, and percentage healthcare-employed residents.

This analysis uses a spatiotemporal model to provide insights into how mobility affects the spread of Coronavirus. The analysis uses a propensity score framework to obtain causal estimates of this effect, while allowing for potential interference between nearby counties. We show via simulation study that the causal estimates from our computationally-efficient approximate model have good statistical properties for data generated from a mechanistic process. In particular, the data suggest that decreases in mobility appear to cause a statistically significant decrease on the number of local Coronavirus cases after roughly 5 weeks. This is significant, as such a long lag might be overlooked in a naive analysis of the relationship between mobility and Coronavirus.

There are several noteworthy limitations of this analysis. The intervention data are Estimates correspond to the expected percentage increase in cases from an increase in mobility 50% above baseline: 100 · {exp(50 · δ i ) − 1}, i = 1, 2. The direct effect estimates for lags 3-5 are 3.5, 15.7, and 24.8, with the latter two significant at the α = 0.05 level. The indirect estimates for lags 3-5 are -9.6, -7.6, and 0.3, though none are significant. Taken from a small minority of smartphone users, they provide a limited view into the distribution of mobility among residents for a given county and week. Moreover, the five relevant measures of mobility provided by Google (which we average over) may have different intervention effects. We also do not have data on other interventions to reduce the spread of Coronavirus, such as wearing masks, statemandated closures of different types of businesses, social distancing, washing hands, etc. These are all plausibly correlated with mobility, so it is possible that the intervention effect of mobility incorporates the causal effect of these measures too.

The measure of connectivity that we use -that of adjacency between counties -is also somewhat crude. For example, while the first viral hotspot in the United States is thought to have been Seattle, it quickly travelled to New York by air. A more comprehensive measure of connectivity and travel between regions may be helpful to incorporate to the model.

The response data also has limitations. Because the availability of Coronavirus tests in the United States has varied over time and space, the number of new positive tests results in a particular county/week might reflect both the prevalence of Coronavirus, but also the availability of tests. One solution to this would be to use COVID-19-related hospitalizations or deaths as a response. However, these are necessarily less common -providing fewer response data -and they would add to the already substantial lag between the moment of infection and positive test.

Finally, as we make the "no unmeasured confounders" assumption, it is possible that there are key confounders that have not been included. Because of this, we include a number of possible confounders, but more can always be done in this regard.

As the pandemic continues, many of these data issues are being addressed: surveys are now being done to ascertain levels of mask-usage and other social distancing practices. Publicly available data on the availability of Coronavirus tests by location are emerging. Further, as the number of cases continues to grow, using hospitalizations or deaths as a response becomes more feasible.

Consequences of kriging and land use regression for PM2.5 predictions in epidemiologic analyses: insights into spatial variability using high-resolution satellite data

A multi-city epidemic model

A fully adaptive numerical approximation for a twodimensional epidemic model with nonlinear cross-diffusion

Joint spatio-temporal analysis of multiple response types using the hierarchical generalized transformation model with application to coronavirus disease 2019 and social distancing

Gaussian process approximations for fast inference from infectious disease data

Modelling the spatial-temporal evolution of the 2009 A/H1N1 influenza pandemic in Chile

Asymptotic behaviour of reaction-diffusion systems in population and epidemic models

Symptoms of Coronavirus

Mobility network models of covid-19 explain inequities and inform reopening

Numerical modelling of an SIR epidemic model with diffusion

Prediction of daily fine particulate matter concentrations using aerosol optical depth retrievals from the Geostationary Operational Environmental Satellite (GOES)

Impact assessment of non-pharmaceutical interventions against coronavirus disease 2019 and influenza in hong kong: an observational study

Quantifying the effect of quarantine control in Covid-19 infectious spread using machine learning

Inferring COVID-19 spreading rates and potential change points for case number forecasts

An interactive web-based dashboard to track covid-19 in real time

Definitive Healthcare: USA Hospital Beds

Generalized propensity score approach to causal inference with spatial interference

Google COVID-19 Community Mobility Reports

Study designs for dependent happenings

A diffusive SI model with Allee effect and application to FIV

Modeling the spread of influenza among cities

The role of the propensity score in estimating dose-response functions

A contribution to the mathematical theory of epidemics

The TVBG-SEIR spline model for analysis of COVID-19 spread, and a Tool for prediction scenarios

A spatial-temporal transmission model and early intervention policies of 2009 A/H1N1 influenza in South Korea

The effect of public health interventions on the spread of influenza among cities

The role of residence times in two-patch dengue transmission dynamics and optimal strategies

Early transmission dynamics in wuhan, china, of novel coronavirus -infected pneumonia

Coronavirus Disease 2019 (COVID-19) in Italy

COVID-19 Time-series Prediction by Joint Dictionary Learning and Online NMF

Machine Learning the Phenomenology of COVID-19 From Early Infection Dynamics

Bayesian Inference of COVID-19 Spreading Rates in South Africa

Sir model with directed spatial diffusion

fields: Tools for spatial data

Continuous and discrete sir-models with spatial distributions

The effect of control strategies to reduce social mixing on outcomes of the covid-19 epidemic in wuhan, china: a modelling study

Covid-19 epidemic analysis using machine learning and deep learning algorithms

A review of spatial causal inference methods for environmental and epidemiological applications

A two-phase epidemic driven by diffusion

Spatial dynamics of airborne infectious diseases

The central role of the propensity score in observational studies for causal effects

Estimating causal effects of treatments in randomized and nonrandomized studies

A structured epidemic model incorporating geographic mobility among regions

Simulating the effect of quarantine on the spread of the 1918-19 flu in central Canada

What do randomized studies of housing mobility demonstrate? causal inference in the face of interference

GSODR: Global Summary Daily Weather Data in R

Space-time covariance functions

Health Care and Social Assistance: NAICS 62

American Community Survey Public Use Microdata Samples

Cross diffusion-induced pattern in an SI model

Exposure to air pollution and COVID-19 mortality in the United States: A nationwide cross-sectional study

Stay-at-home works to fight against covid-19: international evidence from google mobility data

The authors thank Sukanya Bhattacharyya, Can Cui, Matt Miller, Parker Trostle, Laura Wendelberger, Stephen Xu and Zun Yin (North Carolina State University), Yawen Guan (University of Nebraska) and Pulong Ma (Duke University) for help compiling the data. This work was partially supported by NIH grant R01ES031651-01.

Disclaimer: The views expressed in this manuscript are those of the individual authors and do not necessarily reflect the views and policies of the U.S. Environmental Protection Agency. Mention of trade names or commercial products does not constitute endorsement or recommendation for use.Note that the conditional expectation on the last line is a function ofē(t). The second equality follows from the causal consistency assumption, the third from the ignorability assumption on the intervention process, the fourth from the form of model (9), and the fifth from the balancing property (Balancing) of the generalized propensity score. The final equality assumes that the generalized propensity scores components have been included in X andX.