evidence from China (Lau et al., 2020) and Italy (Sjödin et al., 2020 ) that a strict lockdown and practice of social distancing could have a substantial effect on reducing disease transmission, but social distancing has economic, psychological and societal costs (Acemoglu et al., 2020; Atalan, 2020; Grover et al., 2020; Sheridan et al., 2020; Venkatesh and Edirappuli, 2020) . How much social distancing is needed to achieve the desired public health effect? We conducted a pilot study in March to investigate the effect of social distancing during the first week of President Trump's 15 Days to Slow the Spread campaign (March 16-22, 2020) on the influenza-like illness (ILI) percentage two and three weeks later. We tested the causal null hypothesis and found some weak evidence (p-value = 0.08) that better social distancing had an effect on ILI percentage three weeks later. In
Supplementary Material A, we described in detail our pilot study. A protocol of the design and analysis was posted on arXiv (https://arxiv.org/abs/2004.02944) before outcome data were available and analyzed.
In addition to testing a causal null hypothesis, our main interest in this article is to explore the "dose-response relationship" between the degree of social distancing and potential public health outcomes under various degrees of social distancing. Infectious disease experts seemed to express sentiments that the effect of social distancing on public health outcomes might be small or even negligible under a small degree of social distancing, but much more substantial under a large degree of social distancing. In an interview with the British Broadcasting Corporation (BBC Radio 4, 2020), director of the National Institute of Allergy and Infectious Diseases (NIAID), Dr. Anthony S. Fauci said:
"We never got things down to baseline where so many countries in Europe and the UK and other countries did -they closed down to the tune of about 97 percent lockdown.
In the United States, even in the most strict lockdown, only about 50 percent of the country was locked down. That allowed the perpetuation of the outbreak that we never did get under very good control". Perhaps Dr. Fauci was proposing a hypothesis that the treatment dose, i.e., level of social distancing, played a very important role, and the causal effect of social distancing as a public health strategy combating coronavirus transmission is likely to be very different depending on the extent to which it is practiced. As statisticians, we would like to formalize and test the hypothesis concerning a dose-response relationship between social distancing and public health outcomes.
Starting late April and early May, many states in the U.S. started phased reopening. States and local governments differed in their reopening timelines; people in different states and counties also differed in their social mobility during the process: some ventured out; some continued to stay at home as much as possible. Figure 1 plots the 7-day rolling average of percentage change in total distance traveled (data compiled and made available by Unacast ™ ) of all counties in the U.S., from mid-March to late May. It is evident that as many counties started to ease social distancing measures, we saw less reduction in distance traveled; in fact, in many counties, distance traveled started to return to and even supersede the pre-coronavirus level. 1: County-averaged 7-day rolling average (black solid line), middle 50% (dark shade), and middle 90% (light shade) of percentage change in total distance traveled, e.g., −0.35 corresponds to 35% reduction in total distance traveled compared to the pre-coronavirus period. Average percentage change during the first week of President Trump's 15 Days to Slow the Spread campaign is in red and during the initial reopening phases is in green.
In this article, we leverage the county-level social distancing data since phased reopening in the U.S. to study the dose-response relationship between the level of social distancing and its treatment effect on public health outcomes. One potential model is the following dose-response kink model:
Y n (z) = Y n (z * ), ∀z ≤ τ, and Y n (z) − Y n (τ ) = β(z − τ ), ∀z > τ, ∀n = 1, 2, · · · , N, for some τ and β,
where z is a measure of change in social mobility (such that smaller z corresponds to more reduction in total distance traveled and better social distancing), z * a baseline social distancing level, and Y n (z) the potential outcome under dose z. Model (1) states that the potential health-related outcome (e.g., death toll, test positivity rate, etc) remains the same as the potential outcome under the baseline social distancing level when the treatment dose z is less than a certain threshold τ , and then increases at a rate proportional to how much z exceeds the threshold. Model (1) succinctly captures two key features policy makers may be most interested in: τ the minimum dose that "activates" the treatment effect, and β how fast the potential outcome changes as the dose changes after exceeding the threshold. Model (1) may remind readers of the "broken line regression" models in regression analysis; see Zhang and Singer (2010, Chapter 4) . The key difference here is that Model (1) and other dose-response relationships in this article are about the potential outcomes, not the observed outcomes in a regression analysis. There are two closely related public health questions of interest. First, is the dose-response kink model sufficient in explaining the observed data? Second, if it is deemed sufficient, how to make inference for the threshold τ and the proportional effect β?
We have three goals in this article. First, we propose a simple, model-free randomization-based procedure that tests if a structured dose-response relationship, e.g., the dose-response kink model, fits the data in a static setting with a time-independent, continuous or many-leveled treatment dose. To be specific, an empirical researcher posits a structured dose-response relationship that she finds scientifically meaningful, parsimonious, and flexible enough to describe data at hand; our developed procedure can then be applied to test if such a postulated dose-response relationship is sufficient to describe the causal relationship. If the hypothesis is rejected, empirical researchers are then advised to re-examine the scientific theory underpinning the postulated model; otherwise, the model seems a good starting point for data analysis. In this way, our method can be deemed as a model-free "diagnostic test" for a dose-response relationship, and more broadly a test of the underlying scientific theory. In our application, the treatment and outcome are both longitudinal;
our second goal is to generalize the proposed design and testing procedure to the longitudinal setting. We define a notion of cumulative dose for a time-varying treatment dose trajectory, and discuss how to embed observational longitudinal data into an approximate randomized controlled trial in order to permute two treatment trajectories. Finally, we closely examine our assumptions in the context of an infectious disease transmission mechanism and apply the developed design and testing procedure to characterize the dose-response relationship between social distancing and COVID-19 related death toll and case numbers during the reopening phases in the U.S. using county-level data we compiled from sources including Unacast ™ , the United States Census Bureau, and the County Health Rankings and Roadmaps Program (Remington et al., 2015) .
2 Investigating the dose-response relationship via nonbipartite matching 2.1 Observational data with a continuous treatment dose in a static setting
Suppose there are N = 2I units, indexed by n = 1, 2, · · · , N . Each unit is associated with a vector of observed covariates X n , an observed treatment dose assignment Z obs n , and an observed outcome Y obs n . The vector of observed covariates X n are collected before the treatment assignment and not affected by the treatment. Let Z n denote the treatment dose assignment of unit n, Z the set of all possible treatment doses, z ∈ Z a realization of Z n , and |Z| the cardinality of Z. For a binary treatment, |Z| = 2; for a continuous treatment dose, |Z| is an infinite number. In most applications, Z is an ordered set (either partially ordered or totally ordered) with a (partial or total) order defined in light of the application.
Let Y n (z) denote the potential outcome that unit n exhibits under the dose assignment z assuming no interference among units (Rubin, 1980 (Rubin, , 1986 . Each unit n is associated with a possibly infinite array of potential outcomes {Y n (z), z ∈ Z}. We will assume consistency so that Y obs n = Y n (Z obs n ). A causal estimand is necessarily a contrast between potential outcomes. Each unit n is associated with a collection of unit-level causal effects {f n (z, z ) = Y n (z) − Y n (z ), ∀z, z ∈ Z}. Table 1 summarizes all information regarding these N units, where we let Z = {0, 1, 2, · · · } be a countable set for ease of exposition. Table 1 is referred to as a science table in the literature (Rubin, 2005) . In a causal inference problem, the fundamental estimands of interest are the arrays of potential outcomes in Table 1 ; the task of uncovering the arrays of potential outcomes is challenging because one and only one of the potentially infinite array of potential outcomes for each unit is actually observed (Holland, 1986) .
One unique feature of problems with a continuous treatment dose assignment is that the unitlevel causal effect is an infinite set of comparisons between any two potential outcomes Y n (z) and Y n (z ), unlike with a binary treatment where the unit-level causal effect unambiguously refers to a comparison between Y n (1) and Y n (0). Let z * ∈ Z denote an arbitrary reference dose. Observe 
Covariates
, z ∈ Z} is sufficient in summarizing all pairwise comparisons of potential outcomes. With a binary treatment, a "summary causal effect" (Rubin, 2005) is defined as a comparison between Y n (1) and Y n (0) over the same collection of units, e.g., the mean unit-level difference for females.
With a continuous treatment dose, we first summarize the causal effects with a "unit-level doseresponse relationship" for each unit n. For example, one simple unit-level dose-response relationship states that Y n (z) − Y n (z * ) = τ 0 , ∀z ∈ Z; in words, for unit n, the causal effect when comparing treatment dose z to the reference dose z * is equal to a constant τ 0 regardless of the dose z ∈ Z.
We may then summarize such unit-level dose-response relationships for a collection of units. For example, one such summary may state that a structured dose-response relationship f (z; z * , θ) holds for all counties in the U.S.; this summary can be represented by the following null hypothesis:
, for all counties in the U.S. indexed by n, for some θ.
Another summary may state that a constant dose-response relationship holds for counties in the Sun Belt, and a dose-response kink model holds for the rest of the U.S., i.e.,
for all counties in the Sun Belt, for some τ 0 , and
for all counties outside the Sun Belt, for some τ and β.
We first develop a simple, randomization-based testing procedure to assess hypotheses of the form H 1 0 and H 2 0 . The work most relevant to our development is Ding et al. (2016) , who studied testing the existence of treatment effect variation in a randomized controlled trial with a time-independent binary treatment.
In a randomization-based inferential procedure, the potential outcomes (i.e., the infinite collec-tion {Y n (z), ∀z ∈ Z, ∀n} in Table 1 are held fixed and the only probability distribution that enters statistical inference is the randomization distribution that describes the treatment dose assignment.
The key step here is to properly embed the observational data into an approximately randomized experiment (Rosenbaum, 2002 (Rosenbaum, , 2010 Bind and Rubin, 2019) , as we are ready to describe.
2.2 Embedding observational data with a time-independent, continuous treatment into an as-if randomized experiment via nonbipartite matching
In a randomized controlled experiment, physical randomization creates "the reasoned basis" for drawing causal inference (Fisher, 1935) . In the absence of physical randomization as with retrospective observational data, one strategy is to use statistical matching to embed observational data into a hypothetical randomized controlled trial (Rosenbaum, 2002 (Rosenbaum, , 2010 Rubin, 2007; Ho et al., 2007; Stuart, 2010; Bind and Rubin, 2019) by matching subjects with the same (or at least very similar) estimated propensity score or observed covariates and forging two groups that are well-balanced in observed covariates.
One straightforward design to handle observational data with a continuous treatment is to dichotomize the continuous treatment based on some prespecified threshold and create a binary treatment out of the dichotomization scheme. For instance, let Z denote a measure of social distancing; one can define counties with the social distancing measure above the median as the "abovemedian," or "treated" group, and the others as the "below-median," or "control" (or "comparison") group. One may then pair counties in the "above-median" group to those in the "below-median" group via a standard bipartite matching algorithm (for instance, via the R package optmatch by Hansen, 2007) , and test the null hypothesis that social distancing has no effect on the outcome. Such a strategy is often seen in empirical research, probably because of its simplicity; however, dichotomizing the continuous treatment inevitably censors the rich information contained in the original, continuous dose and prevents researchers from studying the dose-response relationship.
To address this limitation, Lu et al. (2001 Lu et al. ( , 2011 proposed optimal nonbipartite matching. In a nonbipartite matching, units with similar observed covariates but different treatment doses are paired. Suppose there are N = 2I units, e.g., counties in the U.S in our application. In the design stage, distances {δ ij , i = 1, · · · , N, j = 1, · · · , N } are calculated between each pair of units and a N × N distance matrix is constructed (Lu et al., 2001 (Lu et al., , 2011 Baiocchi et al., 2010) . Some commonly used distances δ ij include the Mahalanobis distance between observed covariates X i and X j and the rank-based robust Mahalanobis distance. Researchers may further modify the distance to incorporate specific design aspects of the study. For instance, in a study that involves effect modification, researchers are advised to match exactly or near-exactly on the effect modifier (Rosenbaum, 2005; Zubizarreta et al., 2014) , e.g., the geographic location of the county, and such an aspect of design can be pursued by adding a large penalty to δ ij if county i and j are not from the same geographic region.
An optimal nonbipartite matching algorithm then divides these N = 2I units into I nonoverlapping pairs of two units such that the total within-matched-pair distance is minimized.
Nonbipartite matching allows more flexible pairing compared to bipartite matching based on a dichotomization scheme, and preserves the continuous nature of the treatment, which is essential for investigating a dose-response relationship.
Suppose that we have formed I matched pairs of 2 units so that index ij, i = 1, · · · , I, j = 1, 2, uniquely identifies a unit. We follow Rosenbaum (1989) and Heng et al. (2019) and define the following potential outcomes after nonbipartite matching.
Definition 1 (Potential Outcomes After Nonbipartite Matching). Let Z obs i1 ∨Z obs i2 = max(Z obs i1 , Z obs i2 ) and Z obs i1 ∧ Z obs i2 = min(Z obs i1 , Z obs i2 ) denote the maximum and minimum of two observed treatment doses in each matched pair i. We define the following two potential outcomes for each unit ij:
where we abuse the notation and use subscripts T and C to denote the potential outcomes under the maximum and minimum of two observed doses within each matched pair, respectively.
Write F = {X ij , Y T ij , Y Cij , i = 1, · · · , I, j = 1, 2}, where Y T ij and Y Cij are defined in Definition 1, Z obs ∨ = (Z obs 11 ∨ Z obs 12 , · · · , Z obs I1 ∨ Z obs I2 ), and Z obs ∧ = (Z obs 11 ∧ Z obs 12 , · · · , Z obs I1 ∧ Z obs I2 ). As always in randomization inference (Rosenbaum, 2002 (Rosenbaum, , 2010 Ding et al., 2016) , we condition on observed covariates, potential outcomes, and observed dose assignments, i.e., we do not model X or the potential outcomes, and rely on the treatment assignment mechanism to draw causal conclusions.
The law that describes the treatment dose assignment in each matched pair i is
and π i2 = 1 − π i1 . In an ideal randomized experiment, experimenters use physical randomization (e.g., coin flips) to ensure π i1 = π i2 = 1/2: for matched pair i with two treatment doses Z obs i1 and Z obs i2 , a fair coin is flipped; if the coin lands heads, the first unit is assigned Z obs i1 and the second unit Z obs i2 , and vice versa if the coin lands tails. The design stage of an observational study aims to approximate this ideal (yet unattainable) hypothetical experiment by matching units with similar covariates X so that π i1 ≈ π i2 after matching. In this way, nonbipartite matching embeds observational data with a continuous treatment dose into a randomized experiment; this induced randomization scheme will serve as our "reasoned basis" for inferring any causal effect including a dose-response relationship. As is always true with retrospective observational studies, a careful design may alleviate, but most likely never eliminate bias due to the residual imbalance in X or unmeasured confounding variables. The departure from randomization, i.e., π i1 = π i2 , is investigated via a sensitivity analysis (Rosenbaum, 1989 (Rosenbaum, , 2002 (Rosenbaum, , 2010 .
3 Randomization-based inference for a dose-response relationship 3.1 Randomization inference for fixed τ = τ 0 and β = β 0 in the dose-response kink model Endowed with the randomization scheme induced by nonbipartite matching, we now turn to statistical inference. We first consider testing the dose-response kink model for a fixed τ = τ 0 and β = β 0 for all units, i.e.,
Under H τ 0 ,β 0 0,kink , the entire dose-response relationship for subject ij is known up to Y ij (z * ). Fortunately, we do observe one point on the dose-response curve, namely Y ij (Z obs ij ); hence, the entire dose-response curve for the subject ij is fixed, and both potential outcomes
can then be imputed for each unit ij. In matched pair i, for the unit with
Analogously, for the unit with
(3) Table 2 illustrates the imputation scheme by imputing the missing potential outcome for each subject under the null hypothesis H τ 0 ,β 0 0,kink with τ 0 = 0.3 and β 0 = 1. 
For each Y obs ij ∈ Y obs max , define the following transformed outcome:
Let Y obs max = { Y obs ij , i = 1, · · · , I} denote the collection of transformed outcomes, and F tr max (·) its CDF. The null hypothesis H τ 0 ,β 0 0,kink can then be tested by comparing the following Kolmogorov-
evaluated at the observed data to a reference distribution generated based on the imputed science table (e.g., Table 2 ) and enumerating all 2 I possible randomizations: within each matched pair i,
In principle, any test statistic can be combined with this randomization scheme to yield a valid test. We motivate the test statistic (4) in Supplementary Material B.
We illustrate the procedure using the following simple example. We generate I = 200 matched (1) with τ = 1 and β = 0.5. We test the null hypothesis H τ 0 ,β 0 0,kink with τ 0 = 1 and β 0 = 0.5 using the test statistic (4). The left panel of Figure 2 plots the empirical distribution F min (y) and F tr max (y), and t KS (1, 0.5) = 0.065 for the observed data. Instead of enumerating all 2 I = 2 200 possible treatment dose assignments, we draw with replacement 100, 000 samples from all 2 200 possible configurations. The right panel of Figure 2 plots the reference distribution based on these 100, 000 samples. Such a "sampling with replacement" strategy is referred to as a "modified randomization test" in the literature (Dwass, 1957; Pagano and Tritchler, 1983 ) and known to still preserve the level of the test. In this way, a p-value equal to 0.803 is obtained in this simulated dataset and the null hypothesis H τ 0 ,β 0 0,kink with τ 0 = 1 and β 0 = 0.5 is not rejected. The p-value is exact as the procedure does not resort to any asymptotic theory and works in small samples.
Observe that
i.e., H K 0 is a composite hypothesis that is equal to the union of H τ 0 ,β 0 0,kink over all τ = τ 0 and β = β 0 . In other words, the activation dose τ and the slope β are nuisance parameters to be taken into account. One strategy testing H K 0 is to take the supremum p-value over the entire range of (τ, β); another commonly used strategy (Berger and Boos, 1994) is to first construct a confidence Being able to reject H K 0 suggests evidence against the postulated dose-response relationship; otherwise, the model is deemed sufficient to characterize the dose-response relationship for the data at hand. We illustrate the procedure using the following example. We generate I = 200 Figure 3 plots the p-values in log scale against τ 0 and β 0 . The maximum p-value is obtained at τ 0 = 3.8 and β 0 = 0.4 and equal to 0.004. The null hypothesis H K 0 , i.e., the dose-response relationship follows a kink model, can be rejected at level 0.05 for this simulated dataset.
Our discussion above suggests a general model-free, randomization-based framework to test any structured dose-response relationship. Here, we say a dose-response relationship is "structured" if it is characterized by a few structural parameters. Consider the following structured dose-response relationship model:
where z * ∈ Z is a reference dose, and f (· ; z * , θ) is a univariate function that satisfies f (z * ; z * , θ) = 0 and is parametrized by a p-dimensional vector of structural parameters θ ∈ R p . Algorithm 1 summarizes a general procedure testing H dose-response 0 at level α. In Supplementary Material D, we brief discuss and illustrate how to sequentially test a few dose-response relationships ordered in their model complexity.
4 Extension to longitudinal studies with a time-varying treatment 4.1 Treatment dose trajectory and potential outcome trajectory In our application, the treatment dose evolves over time and the public-health-related outcomes, e.g., county-level COVID-19 related death toll, may depend on the treatment dose trajectory.
Formally, let t 0 denote a baseline period and t 1 , t 2 , · · · , t i , · · · , T subsequent treatment periods. : Y ij (z) − Y ij (z * ) f (z; z * , θ) = 0, ∀i, j, for some θ 1. Construct CI θ , a level-γ confidence set for the structural parameter θ;
2. For each θ 0 ∈ CI θ , do the following steps:
(a) Compute the test statistic t obs . For each unit ij with Z ij = Z obs i1 ∨ Z obs i2 , i.e., the unit with maximum dose in each matched pair i, define the following transformed outcome
Let F tr max (·) denote the empirical CDF of { Y obs ij , i = 1, · · · , I} and F min (·) the empirical CDF of the collection of units ij with
(c) Generate a reference distribution. Sample with replacement MC = 100, 000 dose assignment configurations from the 2 I possible configurations. For each sampled dose assignment configuration Z k , calculate t (k) KS (θ 0 ) according to Step (a). Let F θ 0 denote the distribution of { t (k) KS (θ 0 ), k = 1, 2, · · · , MC}; (d) Compute the p-value p θ 0 by comparing t obs to the reference distribution F θ 0 , i.e.,
3. Let p max = sup θ∈CI θ p θ and reject the null hypothesis H dose-response 0 at level α if p max + γ ≤ α.
Fix t 0 ≤ t i ≤ t j and let Z be the set of all possible treatment doses at each time point. Let
denote the random treatment dose trajectory (Robins, 1986; Bojinov and Shephard, 2019) from t i to t j , z t i :t j one realization of Z t i :t j , and Z obs n,t i :t j = (Z obs n,t i , Z obs n,t i +1 , · · · , Z obs n,t j ) the observed treatment dose trajectory of unit n from t i to t j . In our application, t 0 denotes the start of the phased reopening and Z t i :t j the trajectory of daily percentage change in total distance traveled from t i to t j . We are interested in the effect of a sustained period of treatment on some future outcome. We assume that the treatment dose at time t temporally precedes the outcome at time t. Fix a time t and let Y n,t (z t 0 :t ) = Y n,t (z t 0 , z t 1 , · · · , z t ) denote the potential outcome of unit n at time t under the treatment dose trajectory Z n;t 0 :t = z t 0 :t . We assume consistency so that Y obs n,t = Y n,t (Z obs n;t 0 :t ). Finally, we let Y n;t i :t j (z t 0 :t j ) denote unit n's potential outcome trajectory from time t i to t j under the treatment dose trajectory z t 0 :t j .
One unique feature of longitudinal data is that the observed outcome trajectory up to time t−1 may confound the treatment dose at time t; this is particularly true in our application: if the COVID-19 related case and death numbers were high during the last week in a county, then residents may be more wary of the disease and reduce social mobility this week. Following the literature on longitudinal studies, we let L n,t denote the time-dependent covariate process of unit n up to but not including time t; L n,t contains both time-independent covariates X n and time-dependent covariates like the observed outcomes {Y obs n,t 0 , Y obs n,t 1 , · · · , Y obs n,t−1 }. We further assume the sequential randomization assumption (SRA) (Robins, 1998) , which states that conditional on the treatment history up to time t − 1 and covariate process up to time t, the treatment dose assignment at time t is independent of the potential outcome trajectories, i.e., Y n;t 0 :T (z t 0 :T ) |= Z n,t | Z n;t 0 :t−1 = z n;t 0 :t−1 , L n,t , ∀z t 0 :T .
This assumption holds if residents' adopting the social distancing measures at time t depends on
(1) their history of adopting social distancing measures, (2) One general recipe for drawing causal inference from longitudinal data is to model the marginal distribution of the counterfactual outcomes Y n,t (z t 0 :t ), or the marginal joint distribution of Y n;t i :t j (z t 0 :t ), as a function of the treatment trajectory and baseline covariates; see Robins (1986 Robins ( , 1994 Robins et al. (1999 Robins et al. ( , 2000 for seminal works. For example, one simplest model may state that N units are i.i.d. samples from a superpopulation such that the counterfactual mean of the outcome at time t depends on the treatment dose trajectory and the time-independent covariates X through a known functional form g(·), i.e., E[Y t (Z t 0 :t ) | X] = g(Z t 0 :t , X; β), and the interest lies in efficient estimation of the structural parameters β.
In the infectious disease context, modeling the potential outcomes is a daunting task and our interest here lies in testing a structural dose response relationship in a less model-dependent way.
To proceed, we generalize the notion of "dose" from the static to longitudinal setting. Consider the following weighted difference between two treatment dose trajectories z t i :t j and z t i :t j :
where W is a shorthand for the weight function W(t ) = {w(t ) | 0 ≤ w(t ) ≤ 1 and t i ≤t ≤t j w(t ) = 1}. Let z * t i :t j denote a reference trajectory, e.g., z * t i :t j = (−1, −1, · · · , −1) corresponding to 100% reduction in total distance traveled from t i to t j . For each treatment dose trajectory z t i :t j , we define its "cumulative dose" as the weighted difference between z t i :t j and z * t i :t j .
Definition 2 (Cumulative Dose). Let z t i :t j be a realization of the treatment dose trajectory Z t i :t j .
Its cumulative dose with respect to the reference trajectory z * t i :t j and the weight function W is
where · W is defined in (6).
Remark 1. The cumulative dose of a treatment dose trajectory is defined with respect to a reference trajectory and a weight function. The choices of the reference trajectory and weight function should be guided by expert knowledge so that the cumulative dose reflects some scientifically meaningful aspect of the treatment dose trajectory.
A collection of treatment dose trajectories is said to be "W-equivalent" if they have the same cumulative dose with respect to the same weight function and reference trajectory.
Definition 3 (W-Equivalence). Two treatment dose trajectories z t i :t j and z t i :t j are said to be W-equivalent w.r.t.to the reference trajectory z * t i :t j , written as z t
. Treatment dose trajectories that are equivalent to z t i :t j form an equivalence class and is denoted as
Equipped with Definition 2 and 3, we are ready to state a major assumption that facilitates extending a dose-response relationship to longitudinal settings.
Assumption 1 (Potential outcomes under W-equivalence). Let [z t 0 :t ] W be an equivalence class as defined in Definition 3 with respect to · W and a reference trajectory z * t 0 :t . Then unit-level potential outcomes at time t, Y n,t (·), satisfies:
Example 1. In a longitudinal study of the effect of zidovudine (AZT) on mortality, Z t 0 :t represents the AZT dose trajectory from t 0 to t. Let Y n,30 = 1 if unit n dies at time t = 30 and 0 otherwise. Assumption 1 applied to Y n,30 then states that patient n's 30-day mortality status depends on the AZT trajectory from t 0 to t only through some "cumulative dose" captured by CD(z t 0 :t ; z * t 0 :t , W). For instance, Robins et al. (2000) assumes the causal effect of the time-varying zidovudine (AZT) treatment on mortality is only through its aggregate dose, i.e., z * t 0 :t = (0, · · · , 0) and CD(z t 0 :t ; z * t 0 :t , W) = t 0 ≤t ≤t z t .
Remark 2. Though Assumption 1 and its variants are often assumed in the literature on longitudinal studies in order to reduce the number of potential outcomes (Robins et al., 2000, Section 7) , its validity needs to be evaluated on a case-by-case basis. We will evaluate Assumption 1 in the infectious disease dynamics context before invoking it in our application.
We now extend the dose-response relationship from a static setting to a longitudinal setting.
Definition 4 (Unit-Level Dose-Response Relationship in Longitudinal Studies). Let CD(z t 0 :t ; z * t 0 :t , W) be a cumulative dose defined in Definition 2 and f n (·; θ n ) a univariate dose-response model parametrized by θ n such that f n (0; θ n ) = 0. Suppose that Assumption 1 holds. A unit-level dose-response relationship for unit n states that
Remark 3. Observe that when z t 0 :t = z * t 0 :t , the LHS of (8) evaluates to 0 and the RHS evaluates to f n CD(z * t 0 :t ; z * t 0 :t , W); θ n = f n {0; θ n } = 0.
Remark 4. Let z t 0 :t and z t 0 :t be two treatment dose trajectories such that z t 0 :t = z t 0 :t but CD(z t 0 :t ; z * t 0 :t , W) = CD(z t 0 :t ; z * t 0 :t , W). For the dose-response relationship (8) to be well-defined, we necessarily have Y n,t (z t 0 :t ) = Y n,t (z t 0 :t ), which is guaranteed by Assumption 1.
Remark 5. Similar to the static setting considered in Section 2, the dose-response relationship (8) can be thought of as a parsimonious summary of unit-level causal effects from a sustained period of treatment.
Let i = 1, 2, · · · , I be I pairs of two units matched on the covariate process L i1,t = L i2,t but Z obs i1;t 0 :t = Z obs i2;t 0 :t . Units i1 and i2 are each associated with the following two potential outcomes at time t: Y ij,t (Z obs i1;t 0 :t ) and Y ij,t (Z obs i2;t 0 :t ), i = 1, · · · , I, j = 1, 2, in parallel with Definition 1 in the static setting. Write F t = {L ij,t , Y ij,t (Z obs i1;t 0 :t ), Y ij,t (Z obs i2;t 0 :t ), i = 1, · · · , I, j = 1, 2}. Let ij denote the unit with the minimum cumulative dose in matched pair i and ij the other unit, and write Z obs ∧;t 0 :t = {Z obs 1j ;t 0 :t , · · · , Z obs Ij ;t 0 :t }, and Z obs ∨;t 0 :t = {Z obs 1j ;t 0 :t , · · · , Z obs Ij ;t 0 :t }. By iteratively applying the sequential randomization assumption, it is shown in the Supplementary Material E that π i1 =P (Z i1;t 0 :t = Z obs i1;t 0 :t , Z i2;t 0 :t = Z obs i2;t 0 :t | F t , Z obs ∧;t 0 :t , Z obs ∨;t 0 :t ) =P (Z i1;t 0 :t = Z obs i2;t 0 :t , Z i2;t 0 :t = Z obs i1;t 0 :t | F t , Z obs ∧;t 0 :t , Z obs ∨;t 0 :t ) = π i2 = 1/2.
Remark 6. In the static setting, it suffices to match on observed covariates to embed data into an approximate experiment; in the longitudinal setting, one needs to match on the covariate process L t including the time-independent covariates and observed outcomes during the treatment period.
Remark 7. Our framework is different from the balance risk set matching of Li et al. (2001) .
According to Li et al. (2001) et al. (2000) . The entire treatment dose trajectory is the unit to be permuted and our design reflects this aspect.
Consider testing the following dose-response relationship in a longitudinal study:
=f CD(z t 0 :t ; z * t 0 :t , W); θ , ∀i = 1, · · · , I, j = 1, 2, for some θ,
where z * t 0 :t is a reference trajectory, CD(z t 0 :t ; z * t 0 :t , W) a cumulative dose, and f (·; θ) a dose-response relationship of scientific interest. Within each matched pair are two observed treatment dose trajectories Z obs i1;t 0 :t and Z obs i2;t 0 :t . We observe the potential outcome that i1 exhibits under Z obs i1;t 0 :t , i.e., Y i1,t (Z obs i1;t 0 :t ) = Y obs i1,t ; moreover, we are able to impute Y i1,t (·) evaluated at Z obs i2;t 0 :t , under H L 0 and Assumption 1:
Similarly, we have Y i2,t (Z obs i2;t 0 :t ) = Y obs i2,t and can impute: Table 3 summarizes the observed and imputed information. The problem has now been reduced to the static setting, except that instead of permuting the two scalar treatment doses, we now permute two treatment dose trajectories. Randomization-based testing procedure like the one discussed in Section 3 in a static setting can be readily applied to testing (1) θ = θ 0 for a fixed θ 0 and (2) the validity of a postulated dose-response relationship H L 0 . 
One unique aspect of our application is that there is typically a time lag between social distancing and its effect on public-health-related outcomes. We formalize this in Assumption 2.
Assumption 2 (Time Lag). The treatment trajectory is said to have a " -lagged effect" on unit n's potential outcomes at time t if Y n,t (z t 0 , z t 1 , · · · , z t− , z t− +1 , · · · , z t ) = Y n,t (z t 0 , z t 1 , · · · , z t− , z t− +1 , · · · , z t ), for all z t 0 , · · · , z t− , z t− +1 , · · · , z t , and z t− +1 , · · · , z t .
In words, Assumption 2 says that the outcome of interest at time t depends on the entire treatment dose trajectory only through up to time t − . Assumption 2 holds in particular when Y n,t measures the number of people succumbing to the COVID-19 at time t. Researchers estimated that the time lag between contracting COVID-19 and exhibiting symptoms (i.e., the so-called incubation period) had a median of 5.1 days and could be as long as 11.5 days (Lauer et al., 2020) , and the time lag between the onset of the COVID-19 symptoms and death ranged from 2 to 8 weeks (Testa et al., 2020; World Health Organization, 2020) . Therefore, it may be reasonable to believe that the number of COVID-19 related deaths at time t does not depend on social distancing practices days immediately preceding t for some properly chosen . Assumption 2 may be further combined with Assumption 1 to state that unit n's potential outcomes at time t depend on the entire treatment dose trajectory z t 0 :t only via some cumulative dose from time t 0 to t − by defining the cumulative dose with respect to some lag-incorporating weight function W lag that assigns 0 weights to treatment doses immediately preceding time t.
Remark 8. Suppose that the time lag assumption holds for potential outcomes Y n,t (·), · · · , Y n,t+ −1 (·), and let g : R l → R be a function that maps these potential outcomes to an aggregate outcome g{Y n,t (·), · · · , Y n,t+ −1 (·)}. One immediate consequence of Assumption 2 is that g{Y n,t (·), · · · , Y n,t+ −1 (·)} depends on the entire treatment dose trajectory Z t 0 :t+ −1 only via Z t 0 :t−1 ; moreover, we may invoke Assumption 1 and further state that g{Y n,t (·), · · · , Y n,t+ −1 (·)} depends on the entire treatment dose trajectory Z t 0 :t+ −1 only via some cumulative dose of Z t 0 :t−1 . Doseresponse relationships, statistical matching, and testing procedures described in Section 4.3 and 4.4 then hold by replacing Y n,t (·) with the aggregate outcome g{Y n,t (·), · · · , Y n,t+ −1 (·)} where appropriate. Details are provided in Supplementary Material F.
5 Social distancing and COVID-19 during phased reopening: study design
Epidemiologists and infectious disease experts routinely leverage mathematical models, e.g., the SIR (Susceptible-Infected-Recovered) model and its numerous variants, to understand disease transmission dynamics (Brauer et al., 2012; Daley and Gani, 2001) . One simple variant is the following SIRD (Susceptible-Infected-Recovered-Dead) model (Canabarro et al., 2020) :
where N represents the total number of people, S(t), I(t), R(t), and D(t) cumulative number of people susceptible, infected, recovered, and succumbed to the disease at time t, respectively, and κ, α, and γ the rate at which a susceptible person gets infected (transmission rate), an infected person recovers, and an infected person dies, respectively.
Specialize Model (12) to a county n. A fixed set of parameters (κ n , α n , γ n ) would then correspond to one trajectory of the number of cumulative deaths D n (t). In reality, the level of daily social distancing is likely to affect the transmission rate κ n , e.g., a larger reduction in total distance traveled at time t may correspond to a smaller κ n at time t, so that κ n is time-varying and denoted as κ n (t). In this way, different social distancing trajectories z n;t 0 :t correspond to possibly different time-varying transmission rate functions κ n (t), which then induce different outcome trajectories.
The causal question we have been asking, i.e., the effect of the social distancing on public-health-related outcomes, can then be understood as how the daily social distancing trajectory affects the public-health-related outcomes via altering aspects of the disease transmission dynamics.
Disease transmission dynamics and the interplay between social distancing and various aspects of the dynamics are much more complicated in reality. The primary goal of this case study is to grasp a high-level understanding of the causal relationship between social distancing level and publichealth-related outcomes via testing structured dose-response relationships. Understanding a highlevel causal effect, i.e., a dose-response relationship in this context, can be a useful first step toward understanding the causal mechanism. As elaborated in previous sections, our randomization-based testing procedure does not require modeling the often highly complicated disease transmission mechanism or the potential public-health-related outcomes.
In the disease transmission context, Assumption 1 states that social distancing trajectories with the same, properly defined cumulative dose would correspond to similar transmission mechanism, e.g., transmission rate function κ(t) in Model (12), and similar transmission mechanisms from t 0 to T would then correspond to roughly the same public-health-related outcome at time T , e.g., the COVID-19 related deaths per 100, 000 people at time T . In this section, we evaluate this assumption by assessing the relationship between the disease transmission mechanism and potential outcomes under the canonical epidemiological model (12).
The left panel of Figure 4 plots three sampled κ(t) curves from each of the four W-equivalence classes with respect to a weight function W(t) proportional to the density of an exponential distribution of rate 5 (hence more weight on the initial transmission rate) and a reference trajectory that is 0 everywhere, using the xsample function ( of γ and α are made in reference to Canabarro et al. (2020) and calibrated so that the number of deaths per 100, 000 people after 60 days is somewhere between 20 and 60. To put these numbers into context, the total number of deaths due to COVID-19 in the U.S. from May 1st to July 1st is approximately 70, 000, corresponding to 21.2 per 100, 000 people. Figure 4 seems to suggest that the transmission rate curves κ(t) with the same cumulative dose from t = 0 to t = 60 lead to relatively similar death toll trajectory D(t) and hence similar COVID-19 related deaths at time t = 60, i.e., D(60) − D(59). In this regard, Assumption 1 seems at least a reasonable approximation in the case study. 5.3 Data: time frame, granularity, cumulative dose, outcome, and covariate history
The first state in the U.S. that reopened was Georgia at April 24th. We hence consider data from April 27th, the first Monday following April 24th, to August 2nd, the first Sunday in August in the primary analysis. We choose a Monday (April 27th) as the baseline period and a Sunday (August 2rd) as the endpoint because social distancing and public-health-related outcomes data exhibited consistent weekly periodicity (Unnikrishnan, 2020) .
We analyze the data at a county-level granularity. Literature on causal inference typically makes the "no interference" assumption at a certain data granularity, e.g., at the dorm level (Sacerdote, 2001; Li et al., 2019) or the school level (Hong and Raudenbush, 2006) , where interference among units in the same dorm or school is allowed, but among units from different dorms or schools is not. In the case study, our county-level analysis allows interference among individual residents in each county, but not interference among different counties.
We use the county-level percentage change in the total distance traveled compiled by Unacast ™ as the continuous, time-varying treatment dose. We consider a two-month treatment period from April 27th (Monday) to June 28th (Sunday). The primary outcome of interest is the cumulative COVID-19 related death toll per 100, 000 people from June 29th (Monday) to August 2nd (Sunday), a total of five weeks. The county-level COVID-19 case number and death toll are both obtained from the New York Times COVID-19 data repository (The New York Times, 2020).
As discussed in Section 4.4, we match counties similar on covariates, including time-independent covariates and time-dependent covariate processes, in order to embed data into an approximate randomized experiment. Specifically, we matched on the following time-independent baseline covariates: female (%), black (%), Hispanic (%), above 65 (%), smoking (%), driving alone to work (%), flu vaccination (%), some college (%), number of membership associations per 10, 000 people, rural (0/1), poverty rate (%), population, and population density (residents per square mile).
These county-level covariates data were derived from the census data collected by the United States Census Bureau and the County Health Rankings and Roadmaps Program (Remington et al., 2015) .
Moreover, we matched on the number of new COVID-19 cases and new COVID-19 related deaths per 100, 000 people every week from April 20th -26th to June 23th -29th.
A total of 1, 211 matched pairs of two counties were formed using optimal nonbipartite matching (Lu et al., 2001 (Lu et al., , 2011 . We matched exactly on the covariate "rural (0/1)" for later subgroup analysis and balanced all other 32 covariates. We added a mild penalty on the cumulative dose so that two counties within the same matched pair had a tangible difference in their cumulative doses, and added 20% sinks to eliminate 20% of counties for whom no good match can be found (Baiocchi et al., 2010; Lu et al., 2011) . Following the advice in Rubin (2007) , the design was conducted without any access to the outcome data in order to assure the objectivity of the design.
Within each matched pair, the county with smaller cumulative dose is referred to as the "better social distancing" county, and the other "worse social distancing" county. Appendix A.1 shows where the 1, 211 better social distancing counties and the other 1, 211 worse social distancing counties are located in the U.S., and Figure 5 plots the average daily percentage change in total distance traveled and the average daily COVID-19 related death toll per 100, 000 people during the treatment period (April 27th to June 28th) in two groups. It is evident that two groups differ in their extent of social distancing, but are very similar in their daily COVID-19 related death toll per 100, 000 people during the treatment period. Finally, Appendix A.2 summarizes the balance of all 33 covariates in two groups after matching. All variables have standardized differences less than 0.15 and are considered sufficiently balanced (Rosenbaum, 2002 average daily COVID-19 related death toll per 100,000 people (solid lines) in 1, 211 better social distancing counties (dark green) and 1, 211 worse social distancing counties (dark red). We saw a sharp contrast in the level of social distancing but little difference in the COVID-19 related death during this treatment period.
6 Social distancing and COVID-19 during phased reopening: outcome analysis 6.1 Primary analysis: causal null hypothesis regarding the death toll Fix t 0 = April 27th and T = June 28th. Let Z t 0 :T = z t 0 :T denote treatment dose trajectory, z * t 0 :T a reference trajectory equal to the −0.50 for all t 0 ≤ t ≤ T (corresponding to 50% reduction in total distance traveled from April 27th to June 28th), W lag a weight function that assigns equal weight to all t such that t 0 ≤ t ≤ T and 0 otherwise, and a cumulative dose CD(z t 0 :T ; z * t 0 :T , W lag ) defined with respect to z * t 0 :T and W lag . Finally, let Y t (·) denote the potential COVID-19 related deaths at time t, = 35 so that T + corresponds to August 2nd, and invoke Assumption 2 on Y t (·) for all t between T + 1 and T + . As discussed in Remark 8, we consider the aggregate outcome Y agg (·) = g{Y T +1 (·), · · · , Y T + (·)} = T +1≤t≤T + Y t (·). We further invoke Assumption 1 so that Y agg (·) depends on Z t 0 :T = z t 0 :T only via CD(z t 0 :T ; z * t 0 :T , W lag ). Our primary analysis tests the following causal null hypothesis for the 1, 211×2 = 2, 422 counties in our matched samples:
t 0 :T ) = 0, ∀z t 0 :T , ∀i = 1, · · · , I = 1211, j = 1, 2.
This null hypothesis states that the treatment dose trajectory from April 27th to June 28th had no effect whatsoever on the COVID-19 related death toll from June 29th to August 2nd.
The top left panel of Figure 6 plots F min (CDF of the better social distancing counties' observed outcomes) and F tr max (CDF of the worse social distancing counties' observed outcomes under H 0,primary ); we calculate the test statistic t KS = 0.735 and contrast it to a reference distribution generated using 1, 000, 000 samples from all possible 2 1211 randomizations; see the top right panel of Figure 6 . In this way, an exact p-value equal to 2.06 × 10 −4 is obtained and the causal null hypothesis H 0,primary is rejected at 0.05 level. We further conducted a sensitivity analysis to assess the no unmeasured confounding assumption we made in the primary analysis. In a sensitivity analysis, we allowed the dose trajectory assignment probability π i1 and π i2 as in (9) to be biased from the randomization probability and then generated the reference distribution with this biased randomization probability; specifically, we considered a biased treatment assignment model where log(Γ i ) = log{π i1 /π i2 } in each matched pair i was proportional to the absolute difference in the cumulative doses of two units in the pair (π i1 = π i2 = 1/2 and Γ i = 1 in a randomized experiment for all i). We found that our primary analysis conclusion would hold up to Γ i having a median as large as 3.82. See Supplementary Material G.3 for details.
Let Y ij,case,agg denote the cumulative COVID-19 cases per 100, 000 people from June 29th to July 12th, as specified in our pre-analysis plan. We test the following null hypothesis concerning the secondary outcome Y ij,case,agg : H 0,secondary : Y ij,case,agg (z t 0 :T ) − Y ij,case,agg (z * t 0 :T ) = 0, ∀z t 0 :T , ∀i = 1, · · · , I = 1211, j = 1, 2. (14) The exact p-value is less than 10 −5 ; see the bottom panels of Figure 6 . Our result suggests strong evidence that social distancing from April 27th to June 28th had an effect on cumulative COVID-19 cases per 100, 000 people from June 29th to July 12th in our matched samples. cumulative dose (defined w.r.t. z * t 0 :T and W lag ) is less than a threshold τ ; after the cumulative dose exceeds this threshold, the COVID-19 case number increases exponentially at a rate proportional to how much the cumulative dose exceeds the threshold. We tested (15) for different τ = τ 0 and β = β 0 combinations; the maximum p-value obtained at (τ 0 , β 0 ) = (4.8, 10.0) is equal to 0.417 and hence the kink model (15) cannot be rejected. The left panel of Figure 7 plots the level-0.1 and level-0.05 confidence sets of (τ, β). We see the confidence set of τ is tightly centered around 0.5, suggesting that as a county's average percentage change in total distance traveled from April 27th to June 28th increases from −50% to around −5% to 5%, the potential COVID-19 cases number from June 29th to July 12th would remain unchanged; once beyond this threshold, the case number rises exponentially. To further illustrate this doseresponse relationship, the right panel of Figure 7 plots dose-response curves with baseline (i.e., 50% reduction) case number equal to 1 per 100, 000 people for three selected (τ 0 , β 0 ) pairs in the level-0.1 confidence set.
We also conducted subgroup analyses by repeating the primary and secondary analyses described in Section 6.1 to 6.3 on 462 matched pairs of 2 non-rural counties and 749 matched pairs of 2 rural counties. P-values when testing the primary analysis hypothesis H 0,primary concerning the death toll and secondary analysis hypothesis H 0,secondary concerning the case number are 0.004 and less than 10 −5 , respectively, in the non-rural subgroup, and 0.008 and 0.009, respectively, in the rural subgroup. We also allowed a differential dose-response relationship between social distancing and case numbers in rural and non-rural counties and constructed confidence sets for (τ, β) separately; see Figure 8 . A comparison of confidence sets for the non-rural counties and rural counties seemed to suggest that the activation dose required to trigger an exponential increase in case numbers in rural counties might be larger than that in non-rural counties. 
We studied in detail the effect of social distancing during the early phased reopening in the United States on COVID-19 related death toll and case numbers using our compiled county-level data.
To address the statistical challenge brought by a time-dependent, continuous treatment dose trajectory, we developed a design-based framework based on nonbipartite matching to embed observational data with time-dependent, continuous treatment dose trajectory into a randomized controlled experiment. This embedding induces a randomization scheme that we then used to conduct randomization-based, model-free statistical inference for causal relationships, including testing a causal null hypothesis and a structured dose-response relationship.
Upon applying the proposed design and testing procedures to the social distancing and COVID-19 data, we found very strong evidence against the causal null hypothesis and supporting a causal effect of social distancing during the early phases of reopening on both the COVID-19 related death toll and case numbers. We then tested a dose-response kink model featuring both an "activation dose threshold" such that there was no treatment effect when the social distancing level was below this threshold and a "slope" that allowed an exponential growth when the social distancing level exceeded the threshold. This dose-response kink model could not be rejected; in fact, the confidence set around the activation dose threshold was tight and its magnitude suggested that once the total distance traveled returned to or even superseded the pre-coronavirus level, it would have a devastating effect on the COVID-19 case numbers by contributing to an exponential growth.
Finally, in a subgroup analysis where we allowed a differential dose-response relationship and tested the dose-response kink model separately in the non-rural and rural counties, we found that in order to avoid an exponential increase in COVID-19 case numbers, more stringent social distancing would be needed for non-rural counties compared to rural counties.
Supplementary Material A provides details on the pilot study described in Section 1. 