#  A demonstration of Modified Treatment Policies to evaluate Evaluating
shifts in mobility and COVID-19 case rates in U.S. counties: A demonstration
of modified treatment policies for causal inference with continuous exposures

Joshua R. Nugent, A.B., M.Ed.1 Laura B. Balzer, PhD2

(1University of Massachusetts Amherst, School of Public Health and Public
Health Sciences, Department of Biostatistics and Epidemiology, 413-545-4603,
<EMAIL_ADDRESS>of Massachusetts Amherst, School of Public Health
and Public Health Sciences, Department of Biostatistics and Epidemiology,
413-545-4603<EMAIL_ADDRESS>
###### Abstract

Mixed evidence exists of

Background: Previous research has shown mixed evidence on the associations
between mobility data and COVID-19 case rates. Evaluating the potential impact
of decreased mobility on COVID-19 transmission is important, but complicated
by differences between places on factors influencing both behavior and health
outcomes.

Methods: We aimed to evaluate the county-level impact of reducing mobility on
new shifting the distribution of mobility on the growth in COVID-19 cases in
summer/fall 2020 in the United States, while demonstrating the use of modified
treatment policies (MTPs), a recent approachfor causal inference with
continuous exposures. Using MTPs, we examined case rates from June 1 -
November 14, 2020. To do so, we utilized a _modified treatment policy_ (MTP)
approach, which considers the impact of shifting the observed mobility
distribution on an exposure away from its observed value. The MTP approach
facilitates studying the effects of continuous exposures while minimizing
parametric modeling assumptions. Ten mobility indices were selected to capture
several aspects of behavior expected to influence and be influenced by
COVID-19 case rates: Facebook’s proportion of residents staying at home and
number of locations visited; PlaceIQ’s density of people at locations visited;
Descartes Lab’s median and maximum distances traveled; and Google’s indices of
retail/recreation business visits, workplace attendance, public transit use,
and time spent at home. The outcome was defined as the number of newly
reported new cases per 100,000 residents two weeks ahead . We considered ten
mobility indices capturing behaviors expected to influence COVID-19
transmissionof each mobility measure. Primary analyses used targeted minimum
loss-based estimation (TMLE) with Super Learner and considered a Super Learner
ensemble of machine learning algorithms, considering 20+ potential confounders
, including capturing counties’ recent case rates ; as well as social,
economic, health, and demographic variables. For comparison, we also
implemented unadjusted analyses.

Results: Analyses included 1,182 of the most populous U.S. counties, covering
90% of the population. For most weeks considered, unadjusted analyses
suggested strong associations between mobility indices and subsequent new
growth in case rates. However, after confounder adjustment, none of the
indices showed consistent associations under mobilityreduction. We
demonstrated how using MTP with TMLE and Super Learner facilitated the
definition of causal effects for continuous exposures without relying on
arbitrary discretizations or parametric modeling assumptions. after
hypothetical shifts to reduce mobility. Sensitivity analyses varying the
confounder set and time lag between exposure and outcomes yielded similar
results.

Conclusion: Application of MTP facilitated examination of associations between
reductions in mobility and subsequent COVID-19 cases in U.S. counties while
minimizing parametric assumptions. After adjusting for measured confounders,
common mobility metrics were generally not associated with new COVID-19 cases
in Summer and Fall of 2020. While identifiability concerns limit our ability
to make causal claims in this analysis, MTPs are a powerful and underutilized
tool for studying the effects of continuous exposuresin epidemiology and
public health..

##  1 Background

During a pandemic caused by a virus novel to the human immune system, such as
today’s SARS-nCoV-2 (COVID-19), policymakers have only non-pharmaceutical
interventions (NPIs) to rely upon while treatments and vaccines are developed,
tested, and distributed. Many NPIs are designed to reduce the number of close
contacts a person has in order to curtail community spread. However, NPIs are
only effective if they actually create the desired changes in behavior and if
the changes in behavior have a causal impact on disease transmission.

To aid in COVID-19 response, private companies and research groups
[huang_twitter_2020, pepe_covid-19_2020, couture_measuring_2020,
gao_mapping_2020] developed and maintain data sets aiming to quantify mobility
and social distancing in the U.S. and worldwide [google_covid19_nodate,
facebook_movement_nodate]. These mobility indices are generated from
anonymized mobile device tracking data aggregated at some geographic level,
such as ZIP code or county. Despite privacy concerns, these indices may be a
useful tool for researchers and public health officials
[buckee_aggregated_2020]. However, it is not clear which of these indices, if
any, are meaningful proxies for behaviors that are related to COVID-19 spread,
and whether any associations identified persist over time [grantz_use_2020].

Estimating the effect Our aim was to evaluate the potential causal impact of
mobility on new COVID-19 case ratesis complicated by , among other factors, ,
a goal complicated by the presence of confounders, varying epidemic arrival
time across different areas, and changing governmental policies. A brief
review of the literature on mobility and COVID-19 transmission is provided in
Appendix A; our work contributes to that corpus in four critical ways. First,
we applied a modified treatment policy (MTP) to define the causal effect of a
continuous exposure [haneuse_estimation_2013], while avoiding examined the
relationships between mobility and COVID-19 case rates over a long (four
month) timeframe. Second, we employed a _modified treatment policy_ (MTP)
approach, described next, to define relevant causal effects and avoid
arbitrary categorizations of that exposure or strong parametric assumptions to
summarize the causal “dose-response” curve; specifically, we considered
additive and multiplicative shifts of the exposuredistribution. To the best of
our knowledge, this is the first application of MTP in an infectious disease
setting, and, thus, this serves as a demonstration of the novel approach for
epidemiologists studying the effects of continuous exposures. Second, the
continuous exposure. Third, we used targeted minimum loss-based estimation
(TMLE) with Super Learner for statistical estimation and inference, allowing
for flexible adjustment for us to adjust for a rich set of county-level
confounders [laan_targeted_2011, laan_targeted_2018]. Third, we examined the
relationships between mobility and COVID-19 case rates over a longer (four
month) timeframe than previously considered in the literaturewhile relaxing
parametric assumptions used in much of the other modeling research. Finally,
we examined a wide set of mobility indices to determine which, if any, might
be associated with COVID-19 case rates, and hence be possible proxies for the
underlying behaviors that drive COVID-19 spread.

###  1.1 Why use MTP?

All ten of the mobility indices, our exposures of interest, were continuous.
As repeatedly demonstrated in All of our exposure data was continuous. Given
the extensive literature on causal inference for binary exposures, one
approach to examining mobility effects would be to discretize the exposure
data into “high” and / “low” levels and ask what would be the expected
difference in counterfactual COVID-19 case rates if all counties had “high”
versus “low” levels that level of mobility. However, selecting a binary cut
point for this categorization is arbitrary and risks losing information.
Additionally, this approach requires a strong assumption on data support:
there must be a non-zero probability of having a “high” and a “low” level of
mobility within all values of measured adjustment variables
[petersen_diagnosing_2012]. In our application, for example, counties with a
high proportion of essential workers might not be able to reduce their overall
mobility to the specified “low” level [bembom_practical_2007], violating the
positivity assumption .

that all counties must have nonzero probability of receiving the specified
exposure level within all values of measured adjustment variables
[petersen_diagnosing_2012]. As an alternative to categorizing the exposure, we
could target the parameters of a working marginal structural model (MSM),
providing a lower-dimensional summary of the true causal
“[robins_marginal_1998, neugebauer_nonparametric_2007] or directly attempt to
estimate a causal ‘dose-response” curve [robins_marginal_1998,
robins_marginal_2000]. However, even when considered as a “working” model
[neugebauer_nonparametric_2007], MSMs ultimately rely on strong parametric
assumptions to summarize the relationship between the continuous exposure and
’ curve to capture how a continuous measure of mobility impacts the expected
counterfactual outcome and, thus, to define the causal parameter of interest.
In our application, for example, it was unclear a priori whether the
relationship between mobility and counterfactual COVID-19 cases rates was
linear, quadratic, or had some other form.

Several exciting developments can help overcome the challenges associated with
categorizing the continuous exposure or correctly specifying the MSM
[westling_causal_2019, kennedy_nonparametric_2019,
rothenhausler_incremental_2020]; however, each has its own limitations. Causal
isotonic regression [westling_causal_2019], for example, requires a monotonic
relationship between the exposure and outcome, which may not be safe to assume
in our context. Alternatively, employing incremental odds-based propensity
score interventions [kennedy_nonparametric_2019] or incremental causal effects
[rothenhausler_incremental_2020] may be unintuitive and/or difficult to
translate into policy. Finally, dose-response curves estimated with kernel
smoothing [kennedy_non-parametric_2017] may not meet the asymptotic
convergence rates, required for accurate confidence interval constructioncase
rates [westling_causal_2019, kennedy_non-parametric_2017]. However, as before,
there may be challenges with data support from structural and practical
violations of the positivity assumption.

Here, we applied MTPs to nonparametrically define the causal effect of
interest and to demonstrate their use in epidemiology. In the MTPframework,
the research question is translated into a causal effect by considering the
Instead, we took non-parametric approach to specify a causal effect that
minimized parametric assumptions and naturally satisfied positivity
[laan_targeted_2011, laan_targeted_2018]. We defined our research question
using an MTP, which aims to quantify the impact on the counterfactual outcome
of an intervention that shifts the exposure for each observation to a new
level $A=a^{d}$, set by a deterministic function $d(a,w)$ of the observed
“natural” ‘natural’ exposure level $A=a$ and covariates $W=w$
[haneuse_estimation_2013]; for through some shift function $d(A,W)$
[haneuse_estimation_2013]. For example, we might consider an MTP where each
county’s mobility is reduced by the mobility shift amount depended on
population density: sparse counties would be assigned a mobility reduction of
10% . MTPs also improve the relevance of our causal questions by allowingwhile
dense counties would be assigned 20%. However, for simplicity and
demonstration, we ignored covariates when determining shift amounts, instead
examining the expected counterfactual outcome if the mobility distribution
across counties were shifted by an additive or multiplicative constant - say,
if the percentage of people staying in their home all day were increased by 5
points from the observed level in each county.

The MTP approach was attractive in our context because it allows, as per
Haneuse and Rotnitzky [haneuse_estimation_2013], “the set of feasible
treatments for each subject [to depend] on attributes that are not fully
captured by baseline covariates but which are reasonably captured by the
treatment actually received.” Thereby, we can avoid the There is no assumption
that, for example, the densely populated Manhattan (New York County ) could
would be shifted to the mobility level of a rural area, even if some observed
covariate values (say, demographics) were similar. MTPs are similar to
stochastic interventionsFull details of the assumptions and estimation
strategy are given in Section 2. MTP is similar to a _stochastic intervention_
, where the exposure $a^{d}$ is a random draw from a shifted distribution
$A^{d}$ of the observed exposure exposures $A$ [didelez_direct_2012,
munoz_population_2012, stock_nonparametric_1989], possibly depending on
covariate values. The MTP approach has been used in several biomedical
settings [diaz_causal_2020, kamel_relationship_2019, hubbard_time-
dependent_2013] and in at least one public policy context [rudolph_when_2021],
but to our knowledge this is the first time it has been applied to infectious
disease.

##  2 Methods

At the county levelcounty-level, we aimed to evaluate the effect of shifting
mobility on new COVID-19 cases per 100,000 residents two weeks later during
the period from June 1 - November 14, 2020. Data on confirmed COVID-19 cases
were obtained from the New York Times [the_new_york_times_coronavirus_2020]
and USAFacts [noauthor_coronavirus_nodate]. In line with other studies
[badr_association_2020, li_association_2020], the a two-week “leading” period
was selected based on the reported incubation time for COVID-19
[linton_incubation_2020, lauer_incubation_2020, li_early_2020] and to account
for delays in case reporting. In sensitivity analyses, we examined the impact
on one-week leading case rates.

We considered ten 10 mobility indices, generated from mobile phones and
selected to represent a range of behaviors that potentially impact and are
impacted by COVID-19 case rates. Specifically, Google’s “residential” index
[google_covid19_nodate] and Facebook’s “single tile user” and “tiles visited”
indices [facebook_movement_nodate] were expected to capture shelter-in-place
behavior behavior related to sheltering at place at home. Additional indices
from Google [google_covid19_nodate] captured changes in workplace attendance
(“workplaces”), non-essential travel (“retail and recreation”), and use of
public transportation (“transit stations”). Likewise, the “Device Exposure
Index” (DEX) and “Adjusted Device Exposure Index” (DEX-A)
[couture_measuring_2020] captured the density of human traffic at locations
visited. Finally, Descartes Labs’ “m50” and “m50 Index”
[descartes_labs_mobility_nodate] captured distance traveled outside the home
(full index definitions in Appendix B.1). Each of these constructs was
expected to capture distinct elements of social interactions, direct or
passive, that could contribute to COVID-19 transmission. As detailed below, we
studied the impact of shifts in the distribution of each index independently.

The ten indices were only reported in counties where sufficient observations
were available to provide accurate and anonymous data, typically more populous
counties. To address this and make the results more comparable across
different indices, we excluded counties with fewer than 40,000 residents. This
led to the exclusion of 1,948 of 3,139 counties (62%) in our main analysis,
many of them rural and sparsely populated. The 1,182 included counties
represent 90% of the U.S. population. A comparison of the counties
included/excluded is in Appendix B.2. Notably, for interpretation of results,
this inclusion/exclusion criteria changes our target population from all to
just the most populous U.S. counties. For the included counties, Figure E.3
summarizes the distributions of each of the ten mobility indices over the
study period. The medians of the distributions vary only mildly over our study
time-frame, and the overall distributions tend to be either symmetrical or
right-skewed, some of which is due to the definition of the index (for
example, the “DEX” cannot be negative) and some to the inherent variability of
the counties.

Violin plots of observed distributions of the ten mobility indices over time.
Descriptions of the indices are provided in Appendix B.1.

To account for measured differences between counties on factors potentially
influencing mobility and case rates (i.e., confounders), we collected a wide
range of covariates from public sources. These spanned several domains:
socioeconomic (ex: levels of poverty, unemployment, education, median
household income, crowded housing, population density, urbanization); health-
related health (ex: rates of smoking and obesity, pollution, extreme
temperatures); demographic (ex: distribution of age, race, sex); and political
(ex: 2016 presidential election result, mask mandate level). Additionally,
since the one-week-behind case rate could influence that week’s mobility
(e.g., people may be less mobile if the perceived danger is higher) and
subsequent transmission, previous week case rates were also included in the
confounder set $W$. The full list of variables collected and links to their
sources are provided in Appendix B.3.

###  2.1 Data structure, causal model, & causal parameter

To assess impact over time while avoiding holiday- and weekday-level case
reporting idiosyncrasies [noauthor_analysis_nodate], mobility and case data
were binned into weeks using a simple mean. We then repeated our analysis in
cross-sectional snapshots over each week during the period beginning June 1,
2020 and ending November 14, 2020 and separately for each mobility index. More
formally, let $t=\\{1,\ldots,24\\}$ denote the week and $m=\\{1,\dots,10\\}$
denote the mobility index. Then for time $t$ and for mobility index $m$, let
$W_{tm}$ denote Our observed data consisted of our outcome $Y$, the new case
rate; exposure $A$, the mobility level; and the set of potential confounders ,
$A_{tm}$ the observed mobility level, and $Y_{tm}$ the new COVID-19 cases per
100,000 residents for a given county. We $W$. For each week and mobility
index, we assumed that the following nonparametric structural equation model
(NPSEM) characterized the county-level data generating process for week $t$
and mobility index $m$ [pearl_causality_2000]:

$\displaystyle\DIFdel{W_{tm}}$ $\displaystyle\DIFdel{=f_{W_{tm}}(U_{W_{tm}})}$
$\displaystyle\DIFdel{A_{tm}}$
$\displaystyle\DIFdel{=f_{A_{tm}}(W_{tm},U_{A_{tm}})}$
$\displaystyle\DIFdel{Y_{tm}}$
$\displaystyle\DIFdel{=f_{Y_{tm}}(W_{tm},A_{tm},U_{Y_{tm}})}$
$\displaystyle\DIFadd{W}$ $\displaystyle\DIFadd{=f_{W}(U_{W})}$
$\displaystyle\DIFadd{A}$ $\displaystyle\DIFadd{=f_{A}(W,U_{A})}$
$\displaystyle\DIFadd{Y}$ $\displaystyle\DIFadd{=f_{Y}(W,A,U_{Y})}$

where $U_{W_{tm}}$, $U_{A_{tm}}$, and $U_{Y_{tm}}$ $U_{W}$, $U_{A}$, and
$U_{Y}$ are unobserved random variables and $f_{W_{tm}}$, $f_{A_{tm}}$, and
$f_{Y_{tm}}$ $f_{W},f_{A}$, and $f_{Y}$ are unspecified (nonparametric)
functionsthat may vary by mobility index and week. For each county
$i=\\{1,\dots,1182\\}$, we assumed its observed data
$O_{tmi}=(W_{tmi},A_{tmi},Y_{tmi})$ were generated by sampling from a
distribution $\mathbb{P}_{0,tm}$ compatible with this causal model . With
$T=24$ time points and $M=10$ mobility indices, we had $T\times M=240$
observed sets, each consisting of 1non-parametric) functions. This causal
model also encoded our assumption of independence between counties. (We
discuss the plausibility of this assumption below.) For each week and index,
182 observations (counties). As previously discussed, we conducted repeated
cross-sectional analyses, taking each week and mobility index in turn.
(Longitudinal approaches to examine the cumulative effects of mobility shifts
over time were of great interest, but beyond the scope of this demonstration
paper.) Importantly, we included the previous week’s new case levels in the
confounder set for the effect of current mobility on cases two weeks later.

For ease of notation for the remainder of the paper, we will drop the
subscripts denoting time $t$ and index$m$ when discussing the causal model
$W=f_{W}(U_{W});A=f_{A}(W,U_{A});Y=f_{Y}(W,A,U_{Y})$, the observed data
$O_{i}=(W_{i},A_{i},Y_{i})$ for county $i$, and the observed data we assumed
the county-level observed data $O=(W,A,Y)$ were generated by sampling from a
distribution $\mathbb{P}_{0}$ . We emphasize that this was to simplify
notation; we are not introducing new assumptions that the data generating
processes were common across time or exposure. Each causal and statistical
analysis was specific to each week-mobility combinationcompatible with this
causal model.

Following the MTP framework in [diaz_nonparametric_2021], we generated
counterfactual case rates $Y^{d}=f_{Y}(W,A^{d},U_{Y})$ by intervening on then
modified $Y=f_{Y}(W,A^{d},U_{Y})$ in the causal model to shift the observed
exposure by some user-specified function: $A^{d}=d(\cdot)$. While this shift
function $d$ can depend on covariates, we focused on simpler interventions to
shift the observed exposure $A$ by an additive constant $c$ or multiplicative
constant $k$: $d(A,c)=c+A$ generate counterfactual case rates under the
shifted exposure distributions. In our application, shifts did not depend on
the covariates $W$, so simplified the treatment rule to $d(A;M)=MA$ and
$d(A,k)=k\times A$$d(A;C)=C+A$ for shifts by a multiplicative constant $M$ or
additive constant $C$, respectively. By this definition, setting $c=0$ or
$k=1$ recovers the observed exposure $A$ and the observed outcome $Y$.

For each index, we selected the scale and the shift amount We selected shift
amounts to reduce mobility based on the index ’s each index definition and
range of observed values(Figure E.3; Appendix B.1). For example, Facebook’s
“single tile user” is a proportion of users who stayed in a single location
(presumably, their residence) in a day and is thus bounded by 0 and 1. Based
on its observed distribution (pooling over all follow-up weeks), we considered
a relative increase relative increases of 5percentage points in primary
analyses and 7-8 percentage pointsin sensitivity analyses; all shifts were ,
7, and 8 percentage points, truncated at the upper bound of 1. For some
indices, such as the DEX-A, which measures the density of mobile devices
visiting commercial locations, larger multiplicative shifts were reasonable;
an example of a multiplicative reduction shift of 0.75 to bring each county to
75% of its observed mobility is shown in the left panel of Figure 1. As a
third example, consider Google’s “retail and recreation” index, defined as the
percentage of visits to non-essential businesses relative to a pre-COVID
baseline of 100%. Here, we examined additive shifts of -4, -5, and -6
percentage points, bounded by the index limit of -100, shown in the right
panel of Figure 1. A list of shifts examined for each variable is given in
Appendix B.1.

Figure 1: Examples of MTP interventions: a multiplicative shift of 0.75 in the
Example distribution of DEX-A on August 26, 2020 (left) and an additive shift
of -5 percentage points in the distribution of Google’s “retail and
recreation” index on on August 26, 2020 (right)shifts exposures. Observed
values are shown in black, and shifted values in redoverlay.

For each week and Given the selected shift $A^{d}$ for each mobility index, we
specified our causal parameter $\phi$ as the difference in the then specified
the causal parameter as the expected counterfactual COVID-19 case rate under
the shifted mobility and the expected COVID-19 case rate under the observed
mobilityif the mobility index had been shifted to $A^{d}$ versus if mobility
index remained at its observed level $A$:

$\DIFdelbegin\DIFdel{\phi^{\Delta}}\DIFaddend\DIFaddbegin\DIFadd{\psi}\DIFaddend=\mathbb{E}\left[Y\DIFaddbegin\DIFadd{(A}\DIFaddend^{d}\DIFaddbegin\DIFadd{)}\DIFaddend\right]-\DIFdelbegin\DIFdel{\mathbb{E}}\DIFdel{Y}\DIFaddend\DIFaddbegin\DIFadd{\mathbb{E}[Y(A)],}\DIFaddend$
(1)

where the expectation is over the (most populous) U. S. counties of our target
population. with $d(A;M)=MA$ for multiplicative shifts or $d(A;C)=C+A$ for
additive shifts. Setting $M=1$ or $C=0$ in the treatment rule $d$ recovers the
observed exposure $A$ and outcome $Y$.

###  2.2 Identification & the statistical estimand

For the causal parameter (Eq. 1, involving written in terms of a summary
measure of the distribution of counterfactuals) to be identified in terms of
the observed data distribution, several assumptions would be were required,
discussed in more detail in Appendix C [munoz_population_2012,
haneuse_estimation_2013, diaz_nonparametric_2021]:

  * •

No unmeasured confounding: There are no unmeasured common causes of the
mobility index and subsequent case rates, which we can formalize as
$Y^{d}\perp\\!\\!\\!\perp A|W$$Y(A^{d})\perp\\!\\!\\!\perp A|W$. In our
context, this assumption would be violated if, for example, an unmeasured
variable influenced both the county-level mobility and the case rate. One
possible example could be issuing of and compliance with NPIs; an NPI might
reduce mobility and also influence non-mobility-related behavior that is
associated with future case rates. We attempted to improve the plausibility of
this assumption by considering a very large confounder set including known
correlates for our outcome and exposure. Nevertheless, we cannot guarantee
that this assumption holds, and therefore limit our interpretations to
statistical associations rather than causal effects.

  * •

Positivity: If $(a,w)$ is within the support of $A,W$, then $(d(a,c),w)$ for
additive shifts or $(d(a,k),w)$ for multiplicative shifts $(d(a,w),w)$ must
also be within the support of $A,W$. In practice, this means that for any
given time-slice and set of adjustment covariates, there is a positive
probability of finding a county with the same covariate values and a mobility
level matching the shifted value. We attempted to improve the plausibility of
this assumption by first screening for instrumental variables (Appendix B.3)
and then considering small shifts, which avoided extreme density ratios
(described below and in Appendix B.1), while recognizing that this approach
makes our causal effect data-adaptive.

The additional assumptions of independence between

  * •

Independence of counties, consistency, and time-orderingwere already made when
specifying the causal model in Eq. LABEL:Eq:SCM and discussed further in
Appendix C. .

If these identifiability assumptions held, we could specify and focus our
estimation efforts on a statistical estimand that equals the wished-for causal
effect. In the (likely) case that they are not satisfied, we could still
specify and focus our estimation efforts on a statistical estimand that is as
close as possible to the causal parameter. Factoring the joint distribution of
the observed data $\mathbb{P}_{0}$ into
$\mathbb{P}_{0}(O)=\mathbb{P}_{0}(Y|A,W)\mathbb{P}_{0}(A|W)\mathbb{P}_{0}(W)$,
Haneuse and Rotnitzky showed $P_{0}$ into
$P_{0}(O)=P_{0}(Y|A,W)P_{0}(A|W)P_{0}(W)$, it was shown in
[haneuse_estimation_2013] that the statistical estimand corresponding to
expected counterfactual outcome under shift treatment rule $d$,
$\mathbb{E}[Y^{d}]$$\mathbb{E}[Y(A^{d})]$, is given by

$\DIFdelbegin\DIFaddend\DIFaddbegin\begin{aligned}
\psi_{0}(A^{d})=\int\mathbb{E}(Y\mid
A=a^{d},W)\text{d}F_{A,W}(a,w),\end{aligned}\DIFaddend$ (2)

with $\text{d}F_{A,W}(a,w)$ as the joint density of received exposures $A$ and
covariate levels $W$ being integrated over, equivalent to the extended
g-computation formula of Robins et al. [robins_effects_2004]. Here and
throughout, the subscript $0$, indexing the statistical parameters, is used to
emphasize these are functions of the observed data distribution
$\mathbb{P}_{0}$. We refer to $\psi_{0}(A^{d})$ as the “shift parameter.”
Under no shift, the expected outcome was trivially identified as
$\psi_{0}(A)=\mathbb{E}(Y)$. Therefore, Extending this, our statistical
estimand of interest, corresponding to the expected difference in COVID-19
case rates under shifted and observed
mobility($\phi^{\Delta}=\mathbb{E}\left[Y^{d}\right]-\mathbb{E}\left[Y\right]$,
given in Eq. 1), was , is given by $\psi_{\Delta}$:

$\DIFdel{\psi^{\Delta}_{0}=\psi_{0}(A^{d})-\psi_{0}(A)=\int\mathbb{E}(Y\mid
A=a^{d},W=w)\text{d}F_{A,W}(a,w)-\mathbb{E}(Y).}$
$\DIFadd{\psi_{\Delta}=\psi_{0}(A^{d})-\psi_{0}(A)=\int\mathbb{E}(Y\mid
A=a^{d},W)\text{d}F_{A,W}(a,w)-\mathbb{E}_{P_{0}}(Y),}$

All shifts where $Y$ denotes the new case rate and $P_{0}$ the observed data
distribution. In words, our statistical estimand $\psi_{\Delta}$ is the
difference in expectation of new cases per 100,000-residents under a specified
mobility shift $A^{d}$, after adjusting for measured confounders, and the
observed expectation of new cases per 100,000-residents
$\mathbb{E}_{P_{0}}(Y)$. All interventions we examined were reductions in
mobility; hence, under this definition, we expected $\psi_{0}^{\Delta}$ expect
$\psi_{\Delta}$ to be negative if mobility reductions were are associated with
lower case rates. We reiterate that each of the causal and statistical
parameters were specific to the week $t$ and mobility index $m$, but
additional subscripting is suppressed for notational convenience.

###  2.3 Estimation and inference

##  3 Estimation

To estimate the expected outcome under the observed exposure
$\psi_{0}(A)=\mathbb{E}(Y)$$\psi_{0}(A)$, we used the empirical mean outcome.
For estimating the shift parameter $\psi_{0}(A^{d})$, given in Eq. 2, shifted
parameter $\psi_{0}(A^{d})$, we used TMLE [laan_targeted_2011,
laan_targeted_2018]. TMLE is a general framework for constructing
asymptotically linear estimators with an optimal bias-variance tradeoff for a
target parameter. ; in our case, $\psi_{0}({A^{d}})$. Typically, TMLE uses the
factorization of the observed data distribution $P_{0}(O)$ into an outcome
regression $\bar{Q}(A,W)=\mathbb{E}(Y\mid A,W)$
$\bar{Q}(A,W)=\mathbb{E}_{P_{0}}(Y\mid A,W)$ and an intervention mechanism
$g(A\mid W)=P(A\mid W)$. Then, , following which an initial estimate of the
outcome regression $\bar{Q}$ is updated by a fluctuation that is a function of
the intervention mechanism. The update step tilts the preliminary estimates
toward a solution of the efficient influence function estimating equation,
endowing it with desirable asymptotic properties such as double-robustness,
described in Appendix D. Further, it allows practitioners practictioners to
utilize machine learning algorithms in the estimation of $\bar{Q}$ and $g$
while still maintaining valid statistical inferenceunder regularity and
convergence conditions [laan_targeted_2011].

TMLE is a general framework; the exact algorithm method may vary depending on
the data structure and question of interest. In this application, we used the
implementation of Díaz [diaz_nonparametric_2021]and refer the reader to their
paper for the technical details. Here, we . We summarize briefly the
algorithmfor the shift parameter $\psi_{0}(A^{d})$, given in Eq. 2; we again
emphasize that we repeated these analyses for each of the 240 week-mobility
combinations. ; further technical details are in [diaz_nonparametric_2021].

First we defined First we define the density ratio

$\displaystyle r(a,w)$ $\displaystyle=\frac{g^{d}(a|w)}{g(a|w)}$ (3)

where $g^{d}(a|w)$ is the conditional density of $a$ under the shifted
distribution and $g(a|w)$ the conditional density as observed. Unfortunately,
estimating the conditional density $g$ directly can be difficult in high
dimensions. Few data-adaptive algorithms exist [diaz_nonparametric_2021], and
for those that do [benkeser_highly_2016], computation time can be demanding
and data sparsity can make estimates unstable. Luckily, an equivalent estimate
for $r$ can be found by recasting it as a classification problem
[diaz_nonparametric_2021]. In this formulation, a data set consisting of two
copies of each observation is created, one assigned the observed exposure and
one the shifted exposure. Further, each is given an indicator $\Lambda$ for
whether it received the natural or shifted value of treatment. Using Bayes’
rule and the fact that $P(\Lambda=1)=P(\Lambda=0)=.5$,

$\displaystyle r(a,w)=\frac{g^{d}(a|w)}{g(a|w)}$
$\displaystyle=\frac{P(a,w\mid\Lambda=1)}{P(a,w\mid\Lambda=0)}=\frac{P(\Lambda=1\mid
a,w)}{P(\Lambda=0\mid a,w)}$ (4) $\displaystyle=\frac{P(\Lambda=1\mid
a,w)}{1-P(\Lambda=1\mid a,w)}.$

In this formulation, the density ratio estimate can be calculated as a
function of $P(\Lambda=1\mid a,w)$, which can be estimated using any binary
classification method and which we will write as $\tilde{\lambda}$. This
innovation by Díaz [diaz_nonparametric_2021] drastically improved computation
time compared to estimating $g$ directly (as was is done in several algorithms
for stochastic interventions [hejazi_txshift_2020, hejazi_tmle3shift_2021]).

To implement estimation of $\tilde{\lambda}$ and the conditional outcome
regression $\bar{Q}$ flexibly, we used the ensemble machine learning algorithm
Super Learner[laan_super_2007], allowing us to avoid the parametric
assumptions inherent in relying on a single regression model and, instead, to
harness the strengths of multiple approaches, including nonparametric
algorithmsminimize parametric assumptions [laan_super_2007]. Super Learner
creates a convex combination of candidate algorithms chosen through cross-
validation to minimize a user-specified loss function. In our case, the Super
Learner ensemble included generalized linear models, generalized additive
models, gradient boosted decision trees, random forest, and a simple mean.

Given estimates of the nuisance functions $\bar{Q}$ and $\tilde{\lambda}$,
Appendix D outlines With nuisance function estimates in hand, step-by-step
estimation of the shift parameter $\psi_{0}(A^{d})$ and then the association
of interest $\psi_{0}^{\Delta}$. The Appendix also discusses asymptotic
properties of TMLE, under specific regularity and convergence conditions
[laan_targeted_2011], and how those asymptotic properties can be used to
generate Wald-style 95% confidence intervals from the estimated influence
curve of $\hat{\psi}^{\Delta}$and inference for $\psi_{0}(A^{d})$ is outlined
in Appendix D. All statistical analyses were analysis was performed in R
version 4.0.3 [r_core_team_r_2021][r_core_team_r_2020], in particular the
`ltmp_tmle` function in the `lmtp` package [williams_lmtp_2020]. Code and data
for reproduction are available at . is available at the authors’ GitHub
repository.

##  4 Results

Figure 3 shows estimates of $\hat{\psi}^{\Delta}$The following figures show
estimates of $\hat{\psi}_{\Delta}$, the difference between the COVID-19 case
rate that would be expected under the shifted distribution, after adjusting
for measured confounders, and the observed average case rate, with associated
95% confidence intervals (CIs), for each of the mobility indices and weekly
time slices from June 1, 2020 to November 14, 2020. For comparison, we also
show unadjusted estimates, which is defined analogously to
$\hat{\psi}^{\Delta}$$\hat{\psi}_{\Delta}$, simply without any covariates $W$.

Figure 2: Difference in mean predicted case rates two weeks ahead with shifted
distribution vs. observed case rates for weekly time slices from June 1 -
November 14, 2020. Unadjusted analyses show some correlation between mobility
and case rates (negative gray values), but after adjustment (colored values),
there is no consistent association across the set of weeks examined. The +5%
increase corresponds to a reduction in mobility (more people staying at home).

As an example, consider the first panel of Figure 3, showing
$\hat{\psi}^{\Delta}$ Figure 2, showing $\hat{\psi}_{\Delta}$ for “Facebook
single tile users,” which measures the proportion of mobile devices in a
county that stayed in a single location (presumably, a home) each day,
averaged over the week. In the first time-slice (June 1, 2020), the unadjusted
results, in gray, showed a small association between reduced mobility (5%
increase in the proportion of people staying at home) and new reduced case
rates: -3.2 (95% CI: -4.4, -1.9) cases per 100,000 residents two weeks later.
After adjusting for confounders (in red), the association shrunk change in
case rates shrinks to -0.7 (-2.2, 0.7). In the week of September 28, the
results were are similar but more extreme: a 5% increase in the proportion of
people staying at home in the unadjusted analysis suggested suggests a change
of -16.0 (-18.6, -13.3) in new cases per 100,000 residents two weeks later,
while that association was is attenuated to -4.5 (-7.0, -2.0) after confounder
adjustment. Across all of the time slices for that mobility index, we saw see
the same pattern. : Unadjusted estimates of $\hat{\psi}^{\Delta}$ varied, but
were $\hat{\psi}_{\Delta}$ vary, but are negative, as we would expect; lower
mobility appeared appears to be associated with lower new case rates. After
adjustment, however, the estimates shrunk shrink towards zero and many of the
CIs overlapped overlap with zero.

Figure 3: Difference in mean predicted case rates two weeks ahead with shifted
distribution vs. observed case rates for weekly time slices from June 1 -
November 14, 2020. Unadjusted analyses show some correlation between mobility
and case rates (gray values with 95% CIs below the null line of 0), but after
adjustment (colored values), no indices show consistent associations across
the weeks examined. For most indices, the majority of confidence intervals
include 0 after adjustment. Recall that the two stay-at-home measures,
Facebook’s ‘single tile users’ and Google’s ‘residential’, were assigned
positive shifts, which corresponded to a reduction in mobility.

Broadening the scope to all indices in Figure 3, some patterns emergedemerge.
First, almost all unadjusted associations were estimates are larger than the
adjusted ones, although point estimates were though point estimates are
usually in the same direction after adjustment. Second, CIs increased increase
in October and November, suggesting higher outcome variance, more treatment
effect variability, or higher density ratios (perhaps from near practical
positivity violations). Third and most notable, after adjustment for our
confounder set, associations were generally small; , CIs overlapped with zero,
and point estimates showed almost no consistent trends of positive or negative
association.

For example, only six of the indices showed even short-term consistent
associations with the outcome after adjustment, perhaps due to seasonal
effects: Facebook’s ‘single tile users’ metric (Aug, Sept, Oct), PlaceIQ’s
‘DEX’ and ‘DEX-A’ measures of device density (Jun, Aug), Decartes Labs ‘m50
index’ of distance traveled from the home (Aug, Sept, Oct), Google’s workplace
index (Sept, Oct), and Google’s stay-at-home residential index (Sept, Oct).
However, these same indices had uncertain associations in other time periods,
with CIs overlapping zero, and in several instances had point estimates (and
sometimes and CIs) unexpectedly above zero.

Surprisingly, three indices showed positive unadjusted point estimates in June
and July (Facebook ‘tiles visited’ and Google’s ‘retail & recreation’ and
‘residential’ indices), with Facebook’s ‘tiles visited’ index showing an
unusually large swing between summer and fall. After adjustment, the estimates
were strongly attenuated. These data points may be displaying most clearly the
confounding effect of the covariates in our adjustment set. More discussion of
possible explanations for the mixed results are presented in Section 5. We
performed several sensitivity analyses that did not meaningfully change our
results (Appendix E).

##  5 Discussion

To the best of our knowledge, we are the first to apply MTPs to study the
effects of shifting the distribution of a continuous exposure in an infectious
disease setting; specificallya doubly-robust MTP/TMLE approach to provide
insights into a pressing question in infectious disease research.
Specifically, we aimed to answerthe question of whether : Did reductions in
mobility, as captured by the mobile phone data, led to reductions lead to
reduction in subsequent COVID-19 case rates across U.S. counties in the summer
and fall of 2020. 2020? We also hope this paper will serve as a demonstration
of how the MTP approach MTP can be used by other researchers to
nonparametrically define causal effects with continuous exposures; in our
case, we examined the impact of additive and multiplicative shifts of the
observed exposure distribution. More commonly, continuous exposures are
handled by discretizing into binary ‘high’ / ‘low’ categories (which may be
defined arbitrarily, or violate positivity assumptions) , or, alternatively,
the causal “dose-response” curve is summarized with an MSM, relying on strong
parametric assumptions (i.e., a generalized linear model to relate the
expected counterfactual outcome to the continuous exposure). Additionally, by
defining the intervention as a function of the observed exposure value, MTPs
can help us to define causal effects that are better supported by the data.

The MTP approach is also simple to implement in software. The recent
innovation by [diaz_nonparametric_2021] to recast estimation of density ratios
into a binary classification problem overcomes many of the prior computational
challenges for estimation of statistical parameterscorresponding to MTPs and
stochastic interventions. This innovation is integrated into `ltmp` R package,
which also provides doubly-robust and flexible estimation through TMLE with
Super Learner, making it a powerful off-the-shelf option for researchers
seeking to implement MTP. Altogether, we believe the MTP framework is an
attractive minimize (likely unjustifiable) parametric assumptions when
defining causal effects and how TMLE can be applied to estimate the
corresponding statistical parameters. We again note that the causal inference
literature for binary exposures is vast and well-established, especially as
compared to the corresponding literature for continuous exposures. We believe
MTP is a powerful and under-utilized approach to studying the defining effects
for continuous exposures.

Our application of MTP to investigate the relationships findings on how the
associations between mobility and new COVID-19 case rates suggested that
associations varied over time . Although COVID spread vary over time and with
confounder adjustment are in line with the results found in the earlier
research outlined in Appendix A. Though our unadjusted analyses seem to
suggest that decreased mobility was associated with reductions in COVID-19
case rates for the majority of timepoints examined, adjusted analyses were
strongly attenuated with most 95% confidence intervals crossing the null
value. This provides evidence that unadjusted associations results suggest
that the associations between commonly studied indices and COVID-19 case rates
may be highly confounded by county characteristics and current case rates. For
example, our adjustment set included variables for county-level political
lean, which may be a proxy for adherence to NPIs like mask-wearing mandates.
It could appear that counties with lower mobility have fewer new cases, but if
low-mobility counties are also more consistently adhering to mask-wearing due
to underlying political dynamics, the association with mobility would be
attenuated vanish after adjustment. Including one-week lagged case rates in
the adjustment set could also help control for (as we argue is necessary in
the infectious disease context) also demonstrates important confounding by
epidemic arrival time and intensity. Our results contextualize and help
explain findings in the earlier research, outlined in Appendix A, which showed
mixed associations depending on the timeframe, mobility index, geographic
area, and potential confounder set considered.

There are several limitations to our applicationanalysis, some of which may
explain the inconsistent patterns of association seen we see over time and
suggest avenues for future investigation. First, while we chose a timeframe
when COVID testing was widely available, differential ascertainment of case
rates across counties could bias our estimates. Similarly, while we examined
ten distinct mobility indices, we recognize that each is a proxy for
underlying mobility patterns in a given county. While we controlled for
factors potentially related to use of mobile phones and certain apps (e.g.
Facebook), such as age, we cannot rule out the potential for measurement error
in the exposure and therefore the potential for bias in our results. Second,
while county-level mobility and covariate data were are easily accessible,
there might be important sub-county level effects that were are masked when
aggregating to the county-levelcounty level. Performing this analysis at the
zip code or census tract level might reveal stronger or more consistent
associations. Indeed, a minimally adjusted an analysis of zip-code level
mobility data and case rates in five cities [glaeser_how_2020], while finding
a strong association, also reported substantial heterogeneity across space and
time. On the other hand, the policy relevance of such results is complicated
considering that most public heath measures aiming to limit close contacts
(and, indirectly, mobility) are implemented at the county and state level.
Third, our assumption of independence between counties is unlikely to hold,
given state-level policy enactment and spillover effects between counties.
Indeed, the results of spatial and network dependence Network effects can be
unpredictable [lee_network_2021]; , so we hesitate to speculate about the
directions in which they might have biased our point estimates or standard
errors, but they are likely in play.

Fourth, despite our large adjustment set, it is possible that we excluded
important confounders, or that we have unmeasured confounding, as is common
with observational data. For example, we had no data on the strength of a
county’s public health infrastructure or compliance with NPIs, which could
presumably impact both the mobility and case rates.

Fifth, while MTPs relax the positivity assumptions of binary, static
interventions (requiring all levels of the exposure to be supported within all
possible values of adjustment variables [petersen_diagnosing_2012]), MTPs do
require support within the shift value. While reasonable in theory, this is
only practically supported when the specified shifts $A^{d}$ are not “too”
drastic and the adjustment set Related to the dimension of the adjustment set
$W$not “too” large . Here, we examined the largest shifts possible where the
mobility-specific density ratio $r(a,w)$ was consistently less than 10, which
we took as a reasonable cutoff for practical positivity violations. This meant
that examining extreme shifts was not possible, because during the study
period many counties already had significantly reduced their mobility compared
to pre-pandemic baselines. However, sensitivity analyses examining slightly
larger shifts, still somewhat supported by the data, yielded similar results
(Appendix F.2). , practical positivity violations may have also limited our
ability to detect associations. After adjusting for demographics, political
lean, and economic variables, the variability in exposure and outcome across
counties may have been small, limiting the amount of shift we could apply
without generating extreme density ratios and the accompanying large standard
errors, since the density ratio is a key component of the efficiency bound. If
mobility index data were available for the sparsely populated counties that we
excluded, it might be possible to expand our target population and investigate
a wider range of shifts intervention shift values that might detect stronger
or more consistent associations. Indeed, the MTP approach will generally be
more informative when a large range of shifts are supported by the data.

Sixth, Fifth, there are still many unknowns about the dynamics of COVID
transmission over space and timeas well as their interplay with human
behavior. Including one week lagged case rates in our adjustment set may have
helped address some of those unknowns, but is likely not able to account for
all of the complex dynamics of COVID transmission. For example, the return to
schools in Fall fall 2020 may partially explain how several indices show
different behavior in the summer and fall: the adjusted DEX and DEX-A point
estimates change sign and the adjusted m50 estimates shift from minimal
associations to stronger ones. Additionally, while MTPs facilitated our
examination of shifting the mobility distribution, we did not explicitly
discuss the policies, such as non-essential business closures, that could
generate such shifts. Instead in this demonstration paper, we sought to
understand if these mobility indices were consistently associated with
subsequent COVID-19 case rates; they were not. Finally, our choice of a cross-
sectional time-series analysis in weekly slices may have hidden complex
patterns over time. A longitudinal analysis that examines the full scope of
mobility over time and its impact on cumulative cases per capita could provide
a more complete view of the full impact (or lack thereof) of these mobility
indices mobility on case transmission. Doubly-robust longitudinal MTP methods
are available in the `ltmp` R package and offer a promising avenue to answer
those questions.

##  6 Conclusion

Altogether, our applied analysis suggested that county-level mobility indices
are not generally associated with reductions in subsequent COVID-19 case rates
, after adjusting for a wide range of measured confounders. Our analysis,
founded in a causal inference framework [petersen_causal_2014,
petersen_introduction_nodate, wong_evaluating_2021], also highlighted some of
the limitations of prior research in this area, including measurement error,
dependence between observational units, unmeasured confounding, limited data
support, and the many unknowns about COVID-19.

Nonetheless, using Despite these limitations, our analysis provides evidence
that mobility indices may not be a simple proxy measure for future case rates
with airborne infectious diseases. Using MTP allowed us to consider shifts in
mobility supported by our observed data, and TMLE allowed for the leveraging
of machine learning to minimize parametric assumptions while still providing
the basis for valid statistical inference. MTPs can be For continuous
exposures, MTP can be a useful tool for studying the impacts of continuous
exposures in epidemiology and beyondimportant research and policy questions.

##  Appendix A Brief review of the literature on mobility & COVID-19
transmission

Recent papers investigating the associations between mobility, COVID-19
spread, and NPIs can be roughly categorized into one of three bins: work that
uses mobility as an input into a more complex model to forecast COVID-19 cases
and deaths (see, e.g., [chang_mobility_2020, chernozhukov_causal_2021,
karaivanov_face_2020]); work that assesses the effectiveness of, or compliance
with, NPIs, using mobility as an outcome variable (see, e.g.,
[jay_neighbourhood_2020, gupta_tracking_2020, lai_unsupervised_2020, yabe_non-
compulsory_2020, allcott_what_2020, wellenius_impacts_2021]); and finally work
that treats mobility as an exposure variable and new COVID-19 case rates or
deaths as the outcome. This work falls clearly in the last category.

Most papers focusing on the early-pandemic (Spring 2020) timeframe found
strong correlations between mobility and COVID-19 spread [li_association_2020,
badr_association_2020, xiong_mobile_2020, persson_monitoring_2021,
wellenius_impacts_2021]. However, these early studies are potentially
confounded by the spatial dynamics of where the pandemic arrived first and
unmeasured variables such as underlying changes in behavior including mask-
wearing or the propensity to socialize outdoors rather than in closed indoor
environments (see, e.g., the heterogeneity identified by [allcott_what_2020]).
In addition, they may be subject to bias from undercounted cases due to lack
of testing capacity in many parts of the U.S. [noauthor_covid_nodate]. Given
these challenges, it is perhaps unsurprising that the associations between
mobility and COVID-19 cases were strongly attenuated in studies extended
through mid-Summer of 2020 [badr_limitations_2020, gatalo_associations_2020].

Our findings on how the associations between mobility and COVID spread vary
over time and with confounder adjustment are in line with other results. In an
analysis using data from the 25 hardest-hit U.S. counties through mid-April
2020, Badr et al. [badr_association_2020] fit an unadjusted statistical model
for future case growth rates based on a mobility measure using number of trips
made per day compared to a pre-pandemic baseline. While strong associations
were found between changes in their mobility index and county-level case
growth 9-12 days later, a later updated analysis [badr_limitations_2020],
using data through September 16 and expanding the dataset to all U.S.
counties, showed that the association disappeared after April. Similarly,
Gatalo et al. [gatalo_associations_2020], using a slightly different
methodology, found strong correlations between mobility and case growth from
March 27 to April 20, but only a weak correlation from then until July 22.
Another analysis [xiong_mobile_2020], using data through June 9, 2020, treated
the inflow of traffic into a county as the exposure, adjusted for a small
number of socio-demographic covariates, and found a strong relationship
between inflow and new cases 7 days later. Finally, an examination
[li_association_2020] of the correlations between Google’s six mobility
indices (retail & recreation, transit stations, workplaces, residential,
parks, grocery & pharmacy) and 11-day ahead case growth rates at the county
level, stratifying by urban-rural classification, suggested strong
associations in the February to April timeframe, especially in more dense
counties, but those correlations were weaker when the timeframe was extended
to June.

Our work adds to other evidence of confounding in COVID-19 research; for
example, McLaren [mclaren_racial_2020] found that some correlations between
racial/ethnic groups and COVID-19 mortality disappeared once confounding
factors were taken into account, and Wong and Balzer [wong_evaluating_2021]
found that associations between mask mandates and future case rates were
attenuated after confounder adjustment. We suspect this confounding may
explain why our results contradict those of Wellenius et al.
[wellenius_impacts_2021], who found that mobility levels (using many of the
same county-level metrics) in Spring 2020 were strongly associated with COVID
case growth rates two weeks later. Their work differed from this analysis in
that the researchers did not include any socioeconomic covariates in their
parametric model (suggesting a different causal parameter and estimation
approach) nor extend their analysis beyond March and early April.

##  Appendix B Data details

###  B.1 Mobility index sources, definitions, and shifts considered

Below are the source and the description of each of the ten mobility indices
examined. For each, we chose the reported intervention shift value by: (i)
examining the observed exposure distribution, given in Figure E.3 of the main
text; (ii) selecting a set of shifts that seemed to have reasonable support;
and (iii) of those, selecting the largest shift that yielded an estimated
density ratio $\hat{r}(a,w)=\frac{\hat{g}^{d}(a|w)}{\hat{g}(a|w)}<10$. Recall
$g^{d}(a|w)$ is the conditional density of $a$ under the shifted distribution
and $g(a|w)$ the conditional density as observed.

Facebook Movement Range

  * •

Single tile users: Percentage of people who stay in a single ‘Bing tileChange
in movement: Percentage change in the number of ‘Bing tiles’ (roughly 600 m by
600 m) over an entire day. We examined multiplicative shifts of 1.05, 1.07,
and 1.08, and used 1.05 in the primary analysis.

  * •

Tiles visited: Percentage change in the number of Bing tiles a device pinged
in a given day, compared to the baseline average for that weekday in February
2020, excluding Presidents’ Day, aggregated at the county level. We examined
additive shifts of -3 and -4, and used -4 in the primary analysis.

PlaceIQ -4.

  * •

Device exposure index (DEX): Counts the number of distinct devices that
entered the same commercial spaces that day as the target device, then
averages that number across all devices in the county. We examined
multiplicative shifts of 0.75, 0.7, and 0.65, and used 0.75 in the primary
analysis.

  * •

Adjusted DEX (DEX-A): Scaled version of the DEX that adjusts for the
possibility of devices sheltering in place, generating no pings, and hence not
having their zero-DEX values be counted in the overall average. We examined
multiplicative shifts of 0.75, 0.7, and 0.65, and used 0.75 in the primary
analysis.

Descartes Labs

  * •

m50: Median maximum distance traveled from starting location by all tracked
devices in a countyStay put: Percentage of people who stay in a single Bing
tile over an entire day. We examined multiplicative shifts of 0.6, 0.7, and
0.8, and used 0.8 in the primary analysis.

  * •

m50 index: Relative value of m50 relative to median baseline value for that
day of the week between February 2, 2020, and March 7, 2020. We examined
multiplicative shifts of 0.9 and 0.85, and used 0.9 in the primary analysis.
1.05, 1.07, and 1.08.

Google COVID-19 Community Mobility Reports

For Google indices, all of the baseline values are defined as the median value
for that day of the week from January 3, 2020, to February 6, 2020.

  * •

Retail & recreation: Percentage point change in number of visits and length of
stay from baseline for restaurants, malls, movie theaters, and retail stores
(excluding grocery stores and pharmacies). We examined additive shifts of -4,
-5, and -6, and used -5 in the primary analysis. -6.

  * •

Workplace: Percentage point Transit stations: Percentage change in number of
visits and length of stay from baseline for mass transit stations (buses,
subways, etc). We examined additive shifts of -5, -10, and -15.

  * •

Workplaces: Percentage change in number of visits and length of stay from
baseline for an individual’s place of work. We examined additive shifts of -2,
-3, and -4, and used -3 in the primary analysis. -4.

  * •

Residential: Percentage point change from baseline of duration time spent at
one’s residence each day. Since the baseline values are high, the variation in
this index was is smaller than the others. We examined additive shifts of +1,
+2, and +3, and used +2 in the primary analysis. 3.

PlaceIQ DEX

  *   * •

Transit station: Percentage point change in number of visits and length of
stay from baseline for mass transit stations (buses, subways, etc)Device
exposure index (DEX): number of distinct devices that entered the same
commercial spaces that day as the target device, then averages that number
across all devices in the county. We examined multiplicative shifts of 0.75,
0.7, and 0.65.

  * •

Adjusted DEX: Scaled version of the DEX that adjusts for the possibility of
devices sheltering in place, generating no pings, and hence not having their
zero-DEX values be counted in the overall average. We examined additive shifts
of -5, -10, and -15, multiplicative shifts of 0.75, 0.7, and 0.65.

Descartes Labs m50

  * •

m50: Median maximum distance traveled from starting location by all tracked
devices in a county. We examined multiplicative shifts of 0.6, 0.7, and 0.8.

  * •

m50 index: Relative value of m50 relative to median baseline value for that
day of the week between February 2, 2020, and March 7, 2020. We examined
multiplicative shifts of 0.9 and used -5 in the primary analysis0.85.

###  B.2 Comparison of counties included/excluded

All mobility indices require sufficient sample size within a county to report
reliable data on a given day. In low-population counties, this can make the
mobility data sparse or nonexistent over the timeframe in question. Hence, to
improve comparability across different indices and address positivity
violations, we restricted our analysis to counties with 40,000 or more
residents. The table figure below compares the covariate distribution across
the included/excluded counties. While a vast number of counties were excluded,
they represent only 10% of the U.S. population, and , with the exception of
measures of population and 2016 presidential election results, their covariate
means were similar to the counties included.

| All counties | Included | Excluded
---|---|---|---
N | 3,158 | 1,182 | 1,926
Mean population | 105,362 | 250,886 | 16,053
Mean population density (/sq. mi) | 27,167 | 62,366 | 5,608
% Adults without health insurance | 14.3% | 13% | 15%
% Under 18 | 22.3% | 22.6% | 22.2%
% Households below poverty line | 15.2% | 13.8% | 16.1%
% Hispanic | 9.3% | 10.8% | 8.4%
% Women | 49.9% | 50.6% | 49.5%
% Black | 8.9% | 9.9% | 8.5%
% Clinton vote in 2016 election | 31.7% | 38.6% | 27.8%
% Drive alone to work | 79.2% | 80.5% | 78.9%

###  B.3 Adjustment variables and their sources

We gathered a rich set of social, economic, physical, and demographic
variables $W$ from the U.S. Census American Community Survey, U.S. Bureau of
Labor Statistics, U.S. Centers for Disease Control, U.S. Department of
Agriculture, University of Washington Population Health Institute and
Department of Political Science, National Oceanic and Atmospheric
Administration, and MIT Election Data and Science Lab.

Given plausible concerns about data support and desire to avoid controlling
for instrumental variables, we reduced These included the covariates below,
which we screened into the adjustment set $W$ to eight variables for each
index, based on which had the strongest univariate associations after
examining their correlations with the exposure (mobility index) and outcome
(leading new cases per 100,000 people), pooled over all weeks. In sensitivity
analyses (Appendix E), we also tested strategies where (i) covariates were
chosen week-by-week instead of being held constant across all time slices and
(ii) a smaller covariate set was chosen to make possible examinations of
larger shift values without practical positivity violations. These sensitivity
analyses produced similar results to our main analysis. .

The primary analyses included the following covariates in the adjustment set:

  * •

U.S. Census American Community Survey

    * –

Percentage of households below the poverty line

  * •

University of Washington Population Health Institute, via Robert Wood Johnson
Foundation County Health Rankings

    * –

Percentage of people that identify as Black

    * –

Percentage of people that identify as Hispanic

    * –

Percentage of people that drive alone to work

    * –

Percentage of people that are women (note: values under 49% are often an
indicator of a large prison being in the county)

    * –

Percentage of people that are under 18 years of age

    * –

Percentage of adults without health insurance

  * •

MIT Election Data and Science Lab

    * –

Percentage of vote for Hillary Clinton in the 2016 Presidential election

The following covariates were also examined, but not included in the primary
analysiswere removed from the final adjustment set after the screening
process.

  * •

U.S. Census American Community Survey

    * –

Population density

    * –

Natural logarithm of population density

    * –

Average household size

    * –

Percentage of people living below the poverty line

    * –

Percentage of households with income below the poverty line

    * –

Percentage of people in the county who live outside of a metropolitan area

  * •

University of Washington Population Health Institute, via Robert Wood Johnson
Foundation County Health Rankings

    * –

Percentage of adults who smoke

    * –

Percentage of people that identify as Asian

    * –

Percentage of people that are 65 years of age or older

    * –

Percentage of people with more than a high school education

    * –

Index of income inequality

    * –

Natural logarithm of median household income

    * –

Level of air particulate matter (pollution)

  * •

State-level social distancing policies in response to the 2019 novel
coronavirus in the U.S., University of Washington Department of Political
Science

    * –

Mask mandate level (statewide data only)

  * •

Centers for Disease Control and Prevention

    * –

Percentage of people who live in crowded housing

    * –

Percentage of people categorized as obese

    * –

Percentage of people living with diabetes

  * •

National Oceanic and Atmospheric Administration

    * –

Absolute daily high temperature deviation from 70 degrees Fahrenheit Farenheit

  * •

U.S. Bureau of Labor Statistics

    * –

Unemployment rate

  * •

U.S. Department of Agriculture

    * –

Presence of a level 5 meatpacking plant (largest plants)

    * –

Presence of a level 4 meatpacking plant (second largest plants)

    * –

Presence of a level 4 or level 5 meatpacking plant

##  Appendix C Details of identification assumptions

The positivity assumption in this approach is a weaker positivity assumption
than for other causal effects, which require nonzero probability of all
observations receiving a specified level of treatment (say, “low” or “high”
mobility) within values of adjustment variables. While reasonable in theory,
this is viable in practice only when the specified shifts in $A$ are not too
drastic and the adjustment set $W$ not too large. Our main analysis restricted
to shifts that were small enough to make this assumption tenable, while
recognizing that this approach makes our causal effect data-adaptive.

Given plausible concerns about data support and desire to avoid controlling
for instrumental variables, we reduced the adjustment set $W$ to eight
variables based on which had the strongest univariate associations with the
outcome and exposure, pooled over all time slices. In supplementary analyses,
we also tested strategies where (i) covariates were chosen week-by-week
instead of being held constant across all time slices and (ii) a smaller
covariate set was chosen to make possible examinations of larger shift values
without practical positivity violations.

Beyond no unmeasured confounding and positivity, there are three additional
assumptions that are implied by our adoption of the NPSEM. (1) Independence of
counties. This independence assumption is likely violated since there were
state-level policies (e.g., stay-at-home orders) that plausibly created
dependence between counties in a state. The independence assumption also
implies no interfence: The exposure level $a$ for a given county does not
affect the outcomes of the other counties. Given that infectious disease
spread does not respect town, county, state, or country boundaries, this
assumption may be unrealistic. (2) Consistency: If $A=a$ for any county, then
$Y(a)=Y$, and hence the full observed set of outcomes when $A^{d}=A$ is simply
$Y(A^{d})=Y$. This means that the counterfactual outcome for a county with its
observed mobility level is the observed outcome. This assumption is implied by
our definition of counterfactuals as derived quantities through interventions
on the NPSEM. (3) Time-ordering: The confounders $W$ precede the mobility
exposure $A$, which also precedes the case rate outcome $Y$.

##  Appendix D Estimation and inference

Recall from section 3 that we used Super Learner to obtain have nuisance
function estimates for the conditional mean outcome $\bar{Q}$ and the density
ratio $\tilde{\lambda}$ . Given these estimates, we first describe how to
estimate the shift parameter using Super Learner; our goal is to estimate
$\psi_{0}(A^{d})$and then the associational parameter $\psi_{0}^{\Delta}$. The
following steps have been modified for a cross-sectional analysis (i.e., a
single week $t$ and a single index $m$) . The steps below are modified for the
single time point case from the longitudinal version of
[diaz_nonparametric_2021]. In practice, we note that our bounded continuous
outcome will be scaled to $[0,1]$ (see [gruber_targeted_2010]) before
estimation of the conditional mean outcome $\bar{Q}$ and before the steps
below.

  * •

Define and calculate the density ratio under the specified shift
$\hat{r}(a_{i},w_{i})$ for each county weights
$\omega_{i}=\hat{r}(a_{i},w_{i})$ for $i=1,...,n$, calculated from
$\tilde{\lambda}$ estimated via Super Learner .

  * •

Generate initial conditional expectations for each county, denoted
$\bar{Q}(A,W)$.

  * •

Fit the following After scaling our bounded continuous outcome $Y$ to
$\tilde{Y}\in[0,1]$ [gruber_targeted_2010], fit the logistic regression on the
observed (scaled) $Y$ values, using the logit of the initial estimates as an
offset $\tilde{Y}$ values:

$\text{logit}(\DIFdelbegin\DIFdel{{}^{*}(}\DIFaddend\DIFaddbegin\tilde{Y}\DIFadd{|}\DIFaddend
A,W)\DIFdelbegin\DIFdel{)}\DIFaddend=\text{logit}(\bar{Q}(A,W))+\epsilon$

with weights $\hat{r}(a_{i},w_{i})$$\omega_{i}$, calculating the estimated
intercept $\hat{\epsilon}$ , where $\bar{Q}^{*}(A,W)$ is now a targeted
estimate of the conditional mean outcome.

  * •

Use the resulting $\hat{\epsilon}$ value to generate targeted predictions
under a shift, $\bar{Q}^{*}(A^{d},W)$: Update the $\bar{Q}$ estimates as

$\DIFdelbegin\DIFaddend\DIFaddbegin\begin{aligned}
\text{logit}(\tilde{Y}|A^{d},W)&=\text{logit}(\bar{Q}(A^{d},W))+\hat{\epsilon}\\\
\hat{\tilde{Y}}^{d}&=\text{logit}^{-1}\left[\text{logit}(\bar{Q}(A^{d},W))+\hat{\epsilon}\right]\end{aligned}\DIFaddend$

With After unscaling, with the updated estimates $\bar{Q}^{*}(A^{d},W)$ in
hand(after unscaling)$\hat{Y}^{d}$ in hand, define $\hat{\psi}_{tmle}$ as:

$\hat{\psi}_{tmle}=\frac{1}{n}\sum_{i=1}^{n}\DIFdelbegin\DIFdel{{}^{*}(A}\DIFaddend\DIFaddbegin\hat{Y}\DIFaddend^{d}\DIFdelbegin\DIFdel{,W)}\DIFaddend\DIFaddbegin\DIFadd{{}_{i}}\DIFaddend.$

This The TMLE estimator is “doubly-robust” ‘doubly robust’ in that it is a
consistent estimator of $\psi_{0}(A^{d})$ if estimation of either consistent
if either the estimators of $\bar{Q}$ or $\tilde{\lambda}$ is consistentare
consistently estimated and converge at appropriate rates. Further, if both are
consistentmeet those conditions, along with certain regularity conditionsand
multiplicative $\sqrt{n}$ convergence rates [laan_targeted_2011], the TMLE
converges to a normal distribution with mean 0 and variance equal to the
nonparametric non-parametric efficiency bound, as shown in Theorem 3 of
[diaz_nonparametric_2021]. While convergence at the required $n^{-1/4}$ rate
for each nuisance parameter is not guaranteed, the inclusion of a wide variety
of algorithms in our Super Learner ensemble maximizes our chance of achieving
these rates [laan_super_2007]. The efficiency bound for this TMLE depends on
the variance of the outcome conditional on $A$ and $W$, the amount of
treatment effect variability, and the intensity of shift applied, as measured
by the density ratio $r(a,w)$ [diaz_nonparametric_2021].

The inclusion of a wide variety of algorithms in our Super Learner ensemble
improves our chances of meeting the required regularity conditions and, thus,
obtaining an asymptotically linear estimator, which is normally distributed in
the large data limit. Additionally, since recent work
[balzer_demystifying_2021] has highlighted the danger that tree-based
algorithms may not meet required convergence rates or the Donsker condition,
we reran the analyses excluding random forest and gradient boosting from our
Super Learner. The exclusion of these algorithms did not notably change our
results. Future work using cross-validation (a.k.a. cross-fitting) as part of
the estimation process would obviate the Donsker condition and allow for more
aggressive machine learning algorithms to be included.

Using the empirical mean $\bar{Y}=\frac{1}{n}\sum_{i=1}^{n}Y_{i}$ to estimate
$\psi_{0}(A)=\mathbb{E}(Y)$, we obtain an estimate of the statistical
association parameter $\psi_{0}^{\Delta}=\psi_{0}(A^{d})-\psi_{0}(A)$ with to
estimate $\mathbb{E}_{P_{0}}(Y)$, we hence have our desired estimate

$\hat{\psi}\DIFdelbegin\DIFdel{{}^{\Delta}}\DIFaddend\DIFaddbegin\DIFadd{{}_{\Delta}}\DIFaddend=\hat{\psi}_{tmle}-\DIFdelbegin\DIFaddend\DIFaddbegin\DIFadd{\frac{1}{n}\sum_{i=1}^{n}Y_{i}}\DIFaddend.$

Under the above conditions, TMLE is an asymptotically linear estimator of the
shift parameter $\psi_{0}(A^{d})$ with influence curve
$IC^{d}(O)=r(A,W)[Y-\mathbb{E}\left[Y\mid A,W\right]]+\mathbb{E}\left[Y\mid
A^{d},W\right]-\mathbb{E}\left[Y(A^{d})\right]$[diaz_nonparametric_2021].
Furthermore, From this result, variance estimation for $\hat{\psi}$ can be
calculated using the estimated influence curves from $\hat{\psi}_{tmle}$ and
the empirical mean$\bar{Y}$ is a linear estimator of the population mean
$\psi_{0}(A)$ with influence curve $IC(O)=(Y-\bar{Y})$. By applying the Delta
Method, we obtain the influence curve for the TMLE of the associational
parameter $\psi_{0}^{\Delta}$:

$\displaystyle\DIFdel{IC^{\Delta}(O)}$
$\displaystyle\DIFdel{=IC^{d}(O)-IC(O)}$
$\displaystyle\DIFdel{=r(A,W)}[\DIFdel{Y-\mathbb{E}\left[Y\mid
A,W\right]}]\DIFdel{+\mathbb{E}\left[Y\mid
A^{d},W\right]-\mathbb{E}\left[Y(A^{d})\right]-(Y-\bar{Y}).}$
$\DIFadd{Var(\hat{\psi}_{\Delta})=\frac{\frac{1}{n}\sum_{i=1}^{n}(IC(O_{i})-(Y_{i}-\bar{Y}))^{2}}{n},}$

Replacing all of the terms above with the observed data $Y_{i}$ and the
estimates from $\bar{Q}$ and $\tilde{\lambda}$ to calculate
$\hat{IC}^{\Delta}(O_{i})$, the variance of $\hat{\psi}^{\Delta}$ is well-
approximated by the sample variance of the influence curve :

$\DIFdel{Var(\hat{\psi}^{\Delta})=\frac{1}{n}\sum_{i=1}^{n}\hat{IC}^{\Delta}(O_{i})^{2},}$

allowing with the influence curve defined as follows:
$IC(O)=r(A,W)[Y-\mathbb{E}\left[Y\mid A,W\right]]+\mathbb{E}\left[Y\mid
A^{d},W\right]-\mathbb{E}\left[Y(A^{d})\right]$[diaz_nonparametric_2021]. This
allowed for the calculation of the 95% Wald-style confidence interval
$\hat{\psi}^{\Delta}\pm
1.96\sqrt{Var(\hat{\psi}^{\Delta})}$$\hat{\psi}_{\Delta}\pm
1.96\sqrt{Var(\hat{\psi}_{\Delta})}$.

##  Appendix E Sensitivity analyses

###  E.1 Using one-week leading case rates as outcome

###  E.2 Using current-week cases in the adjustment set instead of one-week
lagged cases

###  E.3 Observed exposure distributions over time

##  Appendix F Supplementary analyses

###  F.1 Using a time-varying confounder set

We tested three of the indices allowing for a different confounder set $W$,
screened in each week time-slice by univariate correlations with the exposure
and outcome. As with the main analysis, the adjusted results (shown below for
several shift value possibilities) show inconsistent patterns over time and
CIs frequently overlap with zero.

###  F.2 Using a smaller confounder set with larger shifts

We tested three of the indices with a smaller, time-varying confounder set $W$
(four rather than eight), screened in each week time-slice by univariate
correlations with the exposure and outcome, with larger shift values. As with
the main analysis, the adjusted results (shown below for several shift value
possibilities) show inconsistent patterns over time and CIs frequently overlap
with zero.