# Robust inference for matching under rolling enrollment

Amanda K. Glazer and Samuel D. Pimentel Amanda K. Glazer is a doctoral
candidate and Samuel D. Pimentel is an Assistant Professor in the Statistics
Department at University of California, Berkeley, 367 Evans Hall, Berkeley, CA
94720. Correspondence should be addressed to<EMAIL_ADDRESS>
###### Abstract

Matching in observational studies faces complications when units enroll in
treatment on a rolling basis. While each treated unit has a specific time of
entry into the study, control units each have many possible comparison, or
“pseudo-treatment,” times. Valid inference must account for correlations
between repeated measures for a single unit, and researchers must decide how
flexibly to match across time and units. We provide three important
innovations. First, we introduce a new matched design, GroupMatch with
instance replacement, allowing maximum flexibility in control selection. This
new design searches over all possible comparison times for each treated-
control pairing and is more amenable to analysis than past methods. Second, we
propose a block bootstrap approach for inference in matched designs with
rolling enrollment and demonstrate that it accounts properly for complex
correlations across matched sets in our new design and several other contexts.
Third, we develop a permutation-based falsification test to detect violations
of the timepoint agnosticism assumption, which is needed to permit flexible
matching across time. We demonstrate the practical value of these tools via
simulations and a case study of the impact of short-term injuries on batting
performance in major league baseball.

Keywords: matching, block bootstrap, repeated measures, falsification test

##  1 Introduction

Matching methods attempt to estimate average causal effects by grouping each
treated unit with one or more otherwise similar controls and using paired
individuals to approximate the missing potential outcomes. Assuming that
paired individuals are sufficiently similar on observed attributes and that no
important unobserved attributes confound the comparison, the difference in
outcomes approximates the impact of treatment for individuals in the pair
[Stuart, 2010]. Despite matching’s transparency and intuitive appeal, it faces
complications in datasets containing repeated measures for the same
individuals over time. When only a single time of treatment is present, the
primary challenge is deciding how to construct matching distances from pre-
treatment repeated measures and assess outcomes using post-treatment repeated
measures [Haviland et al., 2008]. The situation is more complex under rolling
enrollment, or staggered adoption, when individuals opt into treatment at
different times [Ben-Michael et al., 2021]. Several authors [Li et al., 2001;
Lu, 2005; Witman et al., 2019; Imai et al., 2020] proceed by matching each
treated unit to the version of the control unit present in the data at the
time of treatment. For example, in Imai et al. [2020]’s reanalysis of data
from Acemoglu et al. [2019] on the impact of democratization on economic
growth, countries undergoing democratizing political reforms are matched to
similar control countries not undergoing such reforms in the same year.

Although this method is logical whenever strong time trends are present, in
other cases it may overemphasize similarity on time at the expense of other
variables. Bohl et al. [2010] study the impact of serious falls on subsequent
healthcare expenditures for elderly adults using patient data from a large
healthcare system. While patients who fall could be matched to patients who
appear similar based on recent health history on the calendar date of the
fall, the degree of similarity in health histories is likely much more
important than the similarity of the exact date at which each patient is
measured. Following this idea, the GroupMatch algorithm [Pimentel et al.,
2020] constructs matches optimally across time, prioritizing matching on
important covariates over ensuring that units are compared at the same point
in time.

Several challenges remain outstanding for matching methods under rolling
enrollment. GroupMatch’s flexible approach relies heavily on a strong
assumption that time itself is not a confounder, and discussion of checking
this assumption has been minimal so far. Even when flexible matching is
warranted, the presence of multiple copies of the same control individual
necessitates a constraint to ensure that a treated unit is not simply paired
to multiple slightly different copies of the same control; several choices of
this constraint exist permitting varying degrees of flexibility, and users
must choose among them. Most importantly, for both GroupMatch and methods that
match exactly in time there is substantial ambiguity about how to conduct
valid inference. When multiple copies of a control individual are forbidden
from appearing in the matched design, randomization inference may be used [Lu,
2005; Pimentel et al., 2020] but no strong guarantees exist outside this
special case.

In what follows, we present several innovations that greatly enhance the
toolkit for matching and treatment effect evaluation under rolling enrollment.
First, we introduce a new matched design called GroupMatch with instance
replacement, which has computational, analytical, and statistical advantages
over existing designs in many common settings. Second, we give a comprehensive
characterization of a new block-bootstrap-based method for inference that
applies broadly across existing methods for matching under rolling enrollment,
including our new design. The block-bootstrap approach was originally
suggested by Imai et al. [2020] and is based on related work in the cross-
sectional case by Otsu and Rai [2017], but until now has not carried any
formal guarantee. Finally, we introduce a falsification test to partially
check the assumption of timepoint agnosticism underpinning GroupMatch’s
validity, empowering investigators to extract evidence from the data about
this key assumption prior to matching. We prove the validity of our bootstrap
method under the most relevant set of constraints on reuse of controls, and we
demonstrate the effectiveness of both the placebo test and the bootstrap
inference approach through simulations and an analysis of injury data in major
league baseball. In particular, the bootstrap method shows improved
performance over linear-regression-based approaches to inference often applied
in similar settings, while making much weaker assumptions.

The paper is organized as follows. Section 2 presents the basic statistical
framework and reviews the GroupMatch framework, inference approaches for
matching designs, and other related literature. In Section 3 we introduce a
new constraint for use of controls in GroupMatch designs, leading to a new
design called GroupMatch with instance replacement. Section 4 presents a block
bootstrap inference approach for matching under rolling enrollment, and
Section 5 evaluates it via simulation. In Section 6 we present a falsification
test for the assumption that time is not a confounder. In Section 7 we apply
our methods to evaluate whether minor injuries impact short-term MLB
performance. Section 8 concludes.

##  2 Statistical framework

###  2.1 Setting and notation

We observe $n$ subjects. For each subject $i$ in the study, we observe
repeated measures $(Y_{i,t},Z_{i,t},\mathbf{X}_{i,t})$ for timepoints
$t=1,\ldots,T$, where $Y_{i,t}$ is an outcome of interest, $Z_{i,t}$ is equal
to the number of timepoints since subject $i$ entered treatment (inclusive of
$t$) or equal to zero if $i$ has not yet been treated, and $\mathbf{X}_{i,t}$
is a vector of covariates. We denote the collection of repeated measures for
each subject $i$ as the trajectory $O_{i}$. We also define $T_{i}$ as the time
$t$ for which $Z_{i,t}=1$ (or $\infty$ otherwise) and $D_{i}$ as an indicator
for $T_{i}<\infty$. We specify a burn-in period of length $L-1$ during which
no individuals are treated (or allow treatment at $t=1$ by setting $L=1$). Let
$Y_{i,t}(z)$ (with $z\leq\max\\{t-L+1,0\\}$) be the potential outcome for unit
$i$ at time $t$ if it had been enrolled in treatment for $z$ timepoints. For
clarity, we focus on outcomes observed immediately following treatment at the
same timepoint. Thus we will restrict attention to the two potential outcomes
$Y_{i,t}(1)$ and $Y_{i,t}(0)$ associated with just having entered treatment at
$t$ and not yet having entered treatment at $t$. The finite sample average
effect of treatment on the treated (ATT) is denoted by $\Delta$:

$\displaystyle\Delta$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}\sum_{t=1}^{T}1\\{t=T_{i}\\}[Y_{i,t}(1)-Y_{i,t}(0)]$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}D_{i}\left[Y_{i,t=T_{i}}(1)-Y_{i,t=T_{i}}(0)\right]\quad\quad\text{for
}z\in\\{0,1\\}.$

We assume that trajectories $O_{i}$ are sampled independently from some
infinite population, although we do not assume independence of observations
within the same trajectory. Defining expectation $E(\cdot)$ with respect to
sampling from this population, we define the population ATT as
$\Delta_{pop}=E(\Delta)$. For future convenience, we also introduce a concise
notation for conditional expectation (again, over the sampling distribution)
of potential outcomes given no treatment through time $t$ and the covariates
observed in the previous $L$ timepoints:

$\mu_{z}^{t}(\mathbf{X})=E[Y_{i,t}(z)|\\{X_{i,t^{\prime}}\\}_{t^{\prime}=t-L+1}^{t^{\prime}=t}=\mathbf{X},T_{i}>t]$

Throughout, we abuse notation slightly by writing $\mu_{0}(\mathbf{X}_{i,t})$
to indicate conditional expectation given the $L$ lagged values of
$\mathbf{X}_{i}$ directly preceding time $t$, inclusive.

###  2.2 Identification assumptions

Pimentel et al. [2020] studied the following difference-in-means estimator in
designs where each treated unit is matched to $C$ control observations.
$M_{it,jt^{\prime}}$ is an indicator for whether subject $i$ at time $t$ has
been matched to subject $j$ at time $t^{\prime}$:

$\displaystyle\hat{\Delta}=\frac{1}{N_{1}}\sum_{i=1}^{n}D_{i}[Y_{i,t=T_{i}}-\frac{1}{C}\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{T}M_{iT_{i},jt^{\prime}}Y_{j,t^{\prime}}]$

Pimentel et al. [2020] show that this estimator is unbiased for the population
ATT under the following conditions:

  1. 1.

Exact matching: matched units share identical values for covariates in the $L$
timepoints preceding treatment.

  2. 2.

$L$-ignorability: conditional on the covariate history over the previous $L$
timepoints and any treatment enrollment in or prior to baseline, an
individual’s potential outcome at a given time is independent of the
individual’s treatment status. Formally,

$D_{i}\perp\\!\\!\\!\perp
Y_{i,t}(0)|Z_{i,t},\\{X_{i,s}\\}_{s=t-L+1}^{t}\text{, }\forall i.$

  3. 3.

Timepoint agnosticism: mean potential outcomes under control do not differ for
any instances with identical covariate histories at different timepoints.
Formally, for any set of $L$ covariate values $\mathbf{X}$,

$\mu_{0}^{t}(\mathbf{X})=\mu_{0}^{t^{\prime}}(\mathbf{X})=\mu_{0}(\mathbf{X})\text{
for any }1\leq t,t^{\prime}\leq T.$

For clarity we drop the $t$ superscript when discussing the conditional
expectation $\mu_{0}(\mathbf{X})$ in what follows, with the exception of
Section 6 where we temporarily consider failures of this assumption.

  4. 4.

Covariate $L$-exogeneity: future covariates do not encode information about
the potential outcome at time $t$ given covariates and treatment status over
the previous L timepoints. Formally,

$(X_{i,1},...,X_{i,T})\perp\\!\\!\\!\perp
Y_{i,t}(0)|(Z_{i,t},\\{X_{i,s}\\}_{s=t-L+1}^{t})\text{, }\forall i.$

  5. 5.

Overlap: given that a unit is not yet treated at time $t-1\geq L$, the
probability of entering treatment at the next time point is neither 0 nor 1
for any choice of covariates over the $L$ timepoints at and preceding $t$.

$0<P(T_{i}=t\mid T_{i}>t-1,X_{i}^{t},\ldots,X_{i}^{t-L+1})<1\quad\quad\forall
t>L$

While not stated explicitly in Pimentel et al. [2020], we note that the
authors rely on an overlap assumption of this type in the proof of their main
result.

The exact matching assumption is no longer needed for asymptotic
identification of the population ATT if we modify the estimator by adding in a
bias correction term. As in Otsu and Rai [2017] and Abadie and Imbens [2011],
we first estimate the conditional mean function $\mu_{0}(\mathbf{X})$ of the
potential outcomes and use this outcome regression to adjust each matched pair
for residual differences in covariates not addressed by matching. As outlined
in Abadie and Imbens [2011], bias correction leads to asymptotic consistency
under regularity conditions on the potential outcome mean estimator
$\widehat{\mu}(\cdot)$ (for further discussion of regularity assumptions on
$\widehat{\mu}_{0}(\cdot)$ see the proof of Theorem 1 in Section A of the
supplemental appendix). Many authors have also documented benefits from
adjusting matched designs using outcome models [Rubin, 1979; Antonelli et al.,
2018]. The specific form of our bias-corrected estimator is as follows:

$\displaystyle\hat{\Delta}_{adj}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{n}D_{i}[(Y_{i,t=T_{i}}-\hat{\mu}_{0}(\mathbf{X}_{i,T_{i}}))-\frac{1}{C}\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{T}M_{iT_{i},jt^{\prime}}(Y_{j,t^{\prime}}-\hat{\mu}_{0}(\mathbf{X}_{i,t^{\prime}}))]$

Large datasets with continuous variables ensure that exact matching is rarely
possible in practice, and in light of this we focus primarily on estimator
$\widehat{\Delta}_{adj}$ in what follows.

##  3 GroupMatch with instance replacement

Before discussing our method for inference in general matched designs under
rolling enrollment, we introduce a new type of GroupMatch design. Pimentel et
al. [2020] described two different designs produced by GroupMatch denoted
Problems A and B, designs we refer to as GroupMatch without replacement and
GroupMatch with trajectory replacement respectively.

  1. 1.

GroupMatch without replacement: each control unit can be matched to at most
one treated unit. If a treated unit is matched to an instance of a control
unit, no other treated unit can match to (any instance) of that control unit.

  2. 2.

GroupMatch with trajectory replacement: each control instance can be matched
to at most one treated unit. Each treated unit can match to no more than one
instance from the same control trajectory. However, different treated units
can match to different instances of the same control trajectory, so a single
control trajectory can contribute multiple distinct instances to the design.

As our chosen names for these designs suggest, their relative costs and
benefits reflect the choice between matching without and with replacement in
cross-sectional settings. As discussed by Hansen [2004], matching without
replacement (in which each control may appear in at most one matched set),
leads to less similar matches compared to matching with replacement (in which
controls can reappear in many matched sets) since in cases where two treated
units both share the same nearest control only one can use it. On the other
hand, matching without replacement frequently leads to estimators with lower
variance than those from matching with replacement, where an individual
control unit may appear in many matched sets, making the estimator more
sensitive to random fluctuations in its response. Thus, one aspect of choosing
between these designs is a choice about how to strike a bias-variance
tradeoff. The other important aspect distinguishing these designs is that
randomization inference, which is based on permuting treatment assignments in
each matched set independently of others, generally requires matching without
replacement.

These same dynamics play out in comparing GroupMatch without replacement and
GroupMatch with trajectory replacement. GroupMatch without replacement ensures
that responses in distinct matched sets are statistically independent (under a
model in which trajectories are sampled independently), allowing for
randomization inference, and ensures that the total weight on observations
from any one control trajectory can sum only to $1/C$, ensuring that the
estimator’s variance cannot be too highly inflated by a single trajectory with
large weight. On the other hand, GroupMatch with trajectory replacement leads
to higher-quality matches and reduced bias in matched pairs.

We suggest a third GroupMatch design which leans even further towards
expanding the potential control pool and reducing bias.

  1. 3.

GroupMatch with instance replacement: Each treated unit can match to no more
than one instance from the same control unit, but control instances can be
matched to more than one treated unit.

GroupMatch with instance replacement is identical to GroupMatch without
trajectory replacement except that it allows repetition of individual
instances within the matched design as well as non-identical instances from
the same trajectory. As such, it is guaranteed to produce higher-quality
matches than GroupMatch without trajectory replacement, but may lead to
higher-variance estimators since individual instances may receive weights
larger than $1/C$. Figure 1 illustrates the these three GroupMatch methods
with a toy example that matches injured baseball players to non-injured
players based on on-base percentage (OBP).

Figure 1: Toy example illustrating the three GroupMatch matching methods. Two
injured baseball players (T1 and T2) are matched 1-1 to non-injured baseball
players (C1a/b and C2a/b) based on player on-base percentage (OBP). Each non-
injured player has two pseudo-injury times or instances. Under GroupMatch
without replacement, T2 must match to an instance in Trajectory 2 because at
most one instance from Trajectory 1 can participate in the match. Under
GroupMatch with trajectory replacement, T2 can match to C1b but not to C1a,
since multiple control instances can be chosen from the same trajectory as
long as they are distinct. Under GroupMatch with instance replacement, both T1
and T2 are able to match to C1a. However, if each treated instance were
matched to two control instances instead of one, GroupMatch with instance
replacement would still forbid either T1 or T2 to match to a second instance
in Trajectory 1.

In practice we view GroupMatch with instance replacement as a more attractive
approach than GroupMatch with trajectory replacement almost without exception.
One reason is that while the true variance of estimators from GroupMatch with
instance replacement may often exceed that of estimators from GroupMatch with
trajectory replacement by a small amount, our recommended approach for
_estimating_ the variance and conducting inference are not able to capture
this difference. As we describe in Section 4, in the absence of a specific
parametric model for correlations within a trajectory, inference proceeds in a
conservative manner by assuming arbitrarily high correlations within a
trajectory (much like the clustered standard error adjustment in linear
regression). Since the variance advantage for GroupMatch with trajectory
replacement arises only when correlations between instances within a
trajectory are lower than one, the estimation strategy is not able to take
advantage of them. This disconnect means that GroupMatch with trajectory
replacement will not generally lead to narrower empirical confidence
intervals, much as variance gains associated with paired randomized trials
relative to less-finely-stratified randomized trials may not translate into
reduced variance estimates [Imbens, 2011].

A second important advantage of GroupMatch with instance replacement is its
computational and analytical tractability relative to the other GroupMatch
designs. One way to implement GroupMatch with instance replacement as a
network flow optimization problem is to remove a set of constraints in
Pimentel et al. [2020]’s Network B (specifically the upper capacity on the
directed edges connected to the sink node), and in Sections 5 and 7 we use
this implementation for its convenient leveraging of the existing groupmatch
package in R. However, much more computationally efficient algorithms are also
possible. Crucially, the removal of the constraint forbidding instance
replacement means that matches can be calculated for each treated instance
without reference to the choices made for other treated units; the $C$ best
matches for a given treated unit are simply the $C$ nearest neighbor instances
such that no two such control instances within the matched set come from the
same trajectory. In principle, this allows for complete parallelization of the
matching routine. On the analytical side, this aspect of the design makes it
possible to characterize the matching algorithm as a generalized form of
nearest neighbor matching, a strategy we adopt in the proof of Theorem 1 to
leverage proof techniques used by Abadie and Imbens [2006] for cross-sectional
nearest neighbor matching. In light of these considerations, we focus
primarily on GroupMatch with instance replacement in what follows, although
the methods derived appear to perform well empirically for other GroupMatch
designs too.

##  4 Block Bootstrap Inference

###  4.1 Inference methods for matched designs

Broadly speaking, there are two schools of thought in conducting inference for
matched designs. One approach, spearheaded by Abadie and Imbens [2006, 2008,
2011, 2012], views the raw data as samples from an infinite population and
demonstrates that estimators based on matched designs (which in this framework
are considered to be random variables, as functions of random data) are
asymptotically normal. Inferences are based on the asymptotic distributions of
matched estimators. A second approach, described in detail in Rosenbaum
[2002a, b] and Fogarty [2020], adopts the perspective of randomization
inference in controlled experiments. Conditional on the structure of the match
and the potential outcomes, the null distribution of a test statistic over all
possible values of the treatment vector is obtained by permuting values of
treatment within matched sets. When matches are exact and unobserved
confounding is absent, strong finite sample guarantees hold for testing sharp
null hypotheses without further assumptions on outcome variables. Asymptotic
guarantees for weak null hypotheses may be obtained too, assuming a sequence
of successively larger finite populations [Li and Ding, 2017]. Well-developed
methods of sensitivity analysis are also available.

As described in Pimentel et al. [2020], while standard methods of inference
may be applied to GroupMatch without replacement, in which control individuals
contribute at most one unit to any part of the match, none have been
adequately developed for GroupMatch with trajectory replacement, in which
distinct matched sets may contain different versions of the same control
individual. For randomization inference, the barrier appears to be quite
fundamental, because permutations of treatment within one matched set can no
longer be considered independently for different matched sets. In GroupMatch
with trajectory replacement, a treated unit receives treatment at one time and
appears in a match only once; if treatment is permuted among members of a
matched set so that a former control now attains treatment status, what is to
be done about other versions of this control unit that are present in distinct
matched sets? We note that similar issues arise when contemplating
randomization inference for general cross-sectional matching designs with
replacement, and we are aware of no solutions for randomization inference even
in this simpler case.

In contrast, the primary issue in applying sampling-based inference to
GroupMatch designs with trajectory replacement is the unknown correlation
structure for repeated measures from a single control individual. The
literature on matching with replacement provides estimators for pairs that are
fully independent [Abadie and Imbens, 2012] and for cases in which a single
observation appears identically in multiple pairs [Abadie and Imbens, 2006],
but not for the intermediate case of GroupMatch with trajectory replacement
where distinct but correlated observations appear in distinct matched sets.
These issues extend beyond the GroupMatch family to any matched design under
rolling enrollment in which control trajectories contribute to multiple
matched sets, including those of Witman et al. [2019] and Imai et al. [2020].

In what follows we give formal guarantees for a sampling-based inference
method appropriate for general matching designs under rolling enrollment
suggested by Imai et al. [2020], which generalizes a recent proposal of Otsu
and Rai [2017] for valid sampling-based inference of cross-sectional matched
studies using the bootstrap. Although the bootstrap often works well for
matched designs without replacement [Austin and Small, 2014], naïve
applications of the bootstrap in matched designs with replacement have been
shown to produce incorrect inferences as a consequence of the failure of
certain regularity conditions [Abadie and Imbens, 2008]. Intuitively, if
matching is performed after bootstrapping the original data, multiple copies
of a treated unit will necessarily all match to the same control unit,
creating a clumping effect not present in the original data. However, Otsu and
Rai [2017] arrived at an asymptotically valid bootstrap inference method for
matching by bootstrapping weighted and bias-corrected functions of the
original observations _after_ matching rather than repeatedly matching from
scratch in new bootstrap samples. We show that a similar bootstrap approach,
applied to entire trajectories of repeated measures in a form of the block
bootstrap, provides valid inference for matched designs under rolling
enrollment. Note that in our formal results we focus on GroupMatch with
instance replacement as the most difficult case, since the designs of Witman
et al. [2019] and Imai et al. [2020] may be understood as restricted special
cases in which matching on time is exact.

###  4.2 Block Bootstrap

In order to conduct inference under GroupMatch with trajectory or instance
replacement we propose a weighted block bootstrap approach. We rearrange the
GroupMatch ATT estimator from Section 2 as follows, letting $K_{M}(i,t)$ be
the number of times the instance at trajectory $i$ and time $t$ is used as a
match.

$\displaystyle\hat{\Delta}_{adj}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}D_{i}[(Y_{i,T_{i}}-\hat{\mu}_{0}(\mathbf{X}_{i,T_{i}}))-\frac{1}{C}\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{T}M_{iT_{i},jt^{\prime}}(Y_{j,t^{\prime}}-\hat{\mu}_{0}(\mathbf{X}_{i,t^{\prime}}))]$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}D_{i}[(Y_{i,T_{i}}-\hat{\mu}_{0}(\mathbf{X}_{i,T_{i}}))-(1-D_{i})\sum_{t=1}^{T}\frac{K_{M}(i,t)}{C}(Y_{i,t}-\hat{\mu}_{0}(\mathbf{X}_{i,t}))]=\frac{1}{N_{1}}\sum_{i=1}^{N}\hat{\Delta}_{i}.$

Because different instances of the same control unit are correlated, we
resample the _trajectory_ -level quantities $\widehat{\Delta}_{i}$ rather than
the instance-level quantities. Since the $\widehat{\Delta}_{i}$ are functions
of the $K_{M}(i,t)$ weights in the original match, we do not repeat the
matching process within bootstrap samples. In particular, we proceed as
follows:

  1. 1.

Fit an outcome regression $\widehat{\mu}_{0}(\cdot)$ for outcomes based on
covariates in the previous $L$ timepoints using only control trajectories.

  2. 2.

Match treated instances to control instances using GroupMatch with instance
replacement. Calculate matching weights $K_{M}(i,t)$ equal to the number of
times the instance at time $t$ in trajectory $i$ appears in the matched
design.

  3. 3.

Calculate the bias-corrected ATT estimator $\widehat{\Delta}_{adj}$.

  4. 4.

Repeat $B$ times:

    1. (a)

Randomly sample $N$ elements $\widehat{\Delta}^{*}_{i}$ with replacement from
$\\{\widehat{\Delta}_{1},\ldots,\widehat{\Delta}_{N}\\}$.

    2. (b)

Calculate the bootstrap bias-corrected ATT estimator
$\widehat{\Delta}_{adj}^{*}$ for this sample of trajectories as follows:

$\widehat{\Delta}^{*}_{adj}=\frac{1}{N_{1}}\sum_{i=1}^{N}\widehat{\Delta}^{*}_{i}$

  5. 5.

Construct a (1 - $\alpha$) confidence interval based on the $\alpha/2$ and
$1-\alpha/2$ percentile of the $\widehat{\Delta}_{adj}^{*}$\- values
calculated from the bootstrap samples.

This method is essentially a block bootstrap, very similar to the method
proposed in Imai et al. [2020]. Our main result below shows the asymptotic
validity of this approach. Several assumptions, in addition to Assumptions 2-5
in Section 2.2, are needed to prove this result. We summarize these
assumptions verbally here, deferring formal mathematical statements to Section
A.1 of the supplemental appendix. First, we require the covariates $X_{i}$ to
be continuous with compact and convex support and a density both bounded and
bounded away from zero. Second, we require that the conditional mean functions
are smooth in $\mathbf{X}$, with bounded fourth moments. In addition, we
require that conditional variances of the treated potential outcomes and
conditional variances of nontrivial linear combinations of control potential
outcomes from the same trajectory are smooth and bounded away from zero. We
also require that conditional fourth moments of potential outcomes under
treatment and linear combinations of potential outcomes under control are
uniformly bounded in the support of the covariates. Finally, we make
additional assumptions related to the conditional outcome mean estimator
$\widehat{\mu}_{0}(\cdot)$, specifically that the $kL$th derivative of the
true conditional mean functions $\mu_{1}^{t}(\cdot)$ and $\mu_{0}(\cdot)$
exist and have finite suprema, and that the $\widehat{\mu}_{(}\cdot)$
converges to $\mu_{0}(\cdot)$ at a sufficiently fast rate. To state the
theorem, we also define

$\sqrt{N_{1}}U^{*}=\frac{1}{N_{1}}\sum_{i=1}^{N}\left(\widehat{\Delta}_{i}^{*}-\widehat{\Delta}_{adj}\right).$

######  Theorem 1.

Under assumptions Mt, W, and R presented in Section A.1 of the supplemental
appendix,

$sup_{r}|Pr\\{\sqrt{N_{1}}U^{*}\leq
r|(\mathbf{Y,D,X})\\}-Pr\\{\sqrt{N_{1}}(\hat{\Delta}_{adj}-\Delta)\leq
r\\}|\xrightarrow{p}0$

as $N\rightarrow\infty$ with fixed control:treated ratio C.

######  Remark 1.

While we focus on the nonparametric bootstrap, this result holds for a wide
variety of other bootstrap approaches including the wild bootstrap and the
Bayesian bootstrap. For required conditions on the bootstrap algorithm see the
proof in Section A.1 of the supplemental appendix.

Assumptions M and R are modeled closely on those of Abadie and Imbens [2006]
and later Otsu and Rai [2017], and our proof technique is very similar to
arguments in Otsu and Rai [2017]. Briefly, $U^{*}$ is decomposed into three
terms which correspond to deviations of the potential outcome variables around
their conditional means, approximation errors for
$\widehat{\mu}_{0}(\mathbf{X})$ terms as estimates of $\mu_{0}(\mathbf{X})$
terms, and deviations of conditional average treatment effects
$\mu_{1}^{t}(\mathbf{X})-\mu_{0}(\mathbf{X})$ around the population ATT
$\Delta$. Regularity conditions from Assumption M ensure that the conditional
average treatment effects converge quickly to the population ATT. Assumption
R, combined with Assumption-M-reliant bounds on the largest nearest-neighbor
discrepancies in $\mathbf{X}$ vectors due originally to Abadie and Imbens
[2006] and adapted to the GroupMatch with instance replacement design, show
that the deviation between $\widehat{\mu}_{0}(\cdot)$ and $\mu_{0}(\cdot)$
disappears at a fast rate. Finally, a central limit theorem applies to the
deviations of the potential outcomes. For details, see Section A of the
supplemental appendix.

###  4.3 Difference-in-Differences Estimator

While we have focused so far on the difference-in-means estimator, Imai et al.
[2020] recommend a difference-in-differences estimator for matched designs
with rolling enrollment in the context of designs that match exactly on time.
We can easily adapt the results of the previous section to show that the block
bootstrap gives valid inference for this setting too.

$\displaystyle\hat{\Delta}_{DiD}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}D_{i}[((Y_{i,T_{i}}-\hat{\mu}_{0}(\mathbf{X}_{i,T_{i}}))-(Y_{i,T_{i}-1}-\hat{\mu}_{0}(\mathbf{X}_{i,T_{i}-1})))-$
$\displaystyle\frac{1}{C}\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{T}M_{iT_{i},jt^{\prime}}((Y_{j,t^{\prime}}-\hat{\mu}_{0}(\mathbf{X}_{i,t^{\prime}}))-(Y_{j,t^{\prime}-1}-\hat{\mu}_{0}(\mathbf{X}_{i,t^{\prime}-1})))]$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}D_{i}[((Y_{i,t=T_{i}}-\hat{\mu}_{0}(\mathbf{X}_{i,T_{i}}))-(Y_{i,t=T_{i}-1}-\hat{\mu}_{0}(\mathbf{X}_{i,T_{i}-1})))-$
$\displaystyle(1-D_{i})\sum_{t=1}^{T}\frac{K_{M}(i,t)}{C}((Y_{i,t}-\hat{\mu}_{0}(\mathbf{X}_{i,t}))-(Y_{i,t-1}-\hat{\mu}_{0}(\mathbf{X}_{i,t-1})))]=\frac{1}{N_{1}}\sum_{i=1}^{N}\hat{\Delta}^{DiD}_{i}$

This estimator requires $L$ lags to be measured at time $T_{i}-1$, so a burn-
in period of length $L$ rather than $L-1$ is needed. Imai et al. [2020] assume
exact matching, which eliminates the need for a bias correction term,
$\hat{\mu}_{0}(\mathbf{X}_{i,t})$, and simplifies our proof of Theorem 1
further.

As described in the previous section, valid inference is possible if we
resample the $\hat{\Delta}^{DiD}_{i}$. Our proof of Theorem 1 presented in
Section A of the supplemental appendix requires mild modification to work for
this difference-in-differences estimator. In particular, the variance
estimators include additional covariance terms. For more details, see Section
A.3 of the supplemental appendix.

##  5 Simulations

We now explore the performance of weighted block bootstrap inference via
simulation. In particular, we investigate coverage and length of confidence
intervals compared to those obtained by conducting parametric inference for
weighted least squares estimators with and without cluster-robust error
adjustment for controls from the same trajectory.

###  5.1 Data Generation

We generate eight covariates, four of them uniform across time for each
individual $i$, (i.e., they take on the same value at every timepoint):

$\displaystyle X_{1,i},X_{3,i},X_{4,i}$ $\displaystyle\sim N(0,1)$
$\displaystyle X_{2,i}\sim N(0,1)\text{ for control units and }$
$\displaystyle X_{2,i}\sim N(0.25,1)\text{ for treated units}$

Additionally, for treated units:

$\displaystyle X_{5},X_{7},X_{8}\sim N(0,1)\text{ and }X_{6}\sim N(0.5,1)$

Four of the covariates are time-varying for control units. For each control
unit, three instances are generated from a random walk process to correlate
their values across time. Formally, for instance $t$ in trajectory $i$,
covariate $j$ is generated as follows:

$\displaystyle X_{j,i,1}$ $\displaystyle\sim N(0,1)$ $\displaystyle X_{j,i,t}$
$\displaystyle=X_{j,i,(t-1)}+\epsilon_{j,i,(t-1)}\text{ for }t=2,3$
$\displaystyle\epsilon_{j,i,1},\epsilon_{j,i,2}$ $\displaystyle\sim
N(0,0.5^{2})$

Fixing $a_{L}=log(1.25)$, $a_{M}=log(2)$, $a_{H}=log(4)$ and $a_{VH}=log(10)$,
and drawing the $\epsilon_{i,t}$ terms independently from a standard normal
distribution, we define out outcome as:

$\displaystyle Y_{i,t}$
$\displaystyle=a_{L}\sum_{j=1}^{4}X_{j,i,t}+a_{VH}X_{5,i,t}+a_{M}(X_{6,i,t}+X_{8,i,t})+a_{H}(X_{7,i,t})+\Delta
D_{i}+\epsilon_{i,t}$ (1)

The outcome for a unit is correlated across time as it is generated from some
time-varying covariates. Each simulation consists of 400 treated and 600
control individuals. We consider 1:2 matching. The true treatment effect,
$\Delta$, is 0.25.

We consider two alternative ways of generating the continuous outcome variable
besides model (1). First, we add correlation to the error terms within
trajectories. Specifically, the $\epsilon_{i,t}$s for a given trajectory $i$
are generated from a normal distribution with mean 0 and covariance matrix
with off diagonal values of 0.8. Second, in addition to the correlated error
terms, we square the $X_{2,i,t}$ term in the model, so it is no longer linear.

We compare the bias-corrected block bootstrap approach outlined in Section 4
to the confidence intervals obtained from weighted least squares (WLS)
regression and WLS with clustered standard errors. We choose to compare to WLS
because this is commonly recommended in matching literature [Ho et al., 2007;
Stuart et al., 2011]. However, Abadie and Spiess [2021] pointed out that
standard errors from regression may be incorrect due to dependencies among
outcomes of matched units, and identified matching with replacement as a
setting in which these dependencies are particularly difficult to correct for.
Our simulation results suggest that these difficulties carry over into the
case of repeated measures. It is worth noting that the standard functions in R
used to compute WLS with matching weights such as lm and Zelig (which calls
lm), compute biased standard error estimates in most settings. See Section B
of the supplemental appendix for details.

###  5.2 Results

Tables 1 and 2 show the coverage and average 95% confidence interval (CI)
length, respectively, of WLS, WLS cluster, and bootstrap bias-corrected
methods for each of our three simulation settings under 10,000 simulations. As
misspecification increases the bootstrap method is substantially more robust
(although under substantial misspecification the bias-corrected method also
fails to achieve nominal coverage). While the bootstrap confidence intervals
are generally slightly wider than the WLS and WLS cluster confidence
intervals, this is to be expected as the wider confidence intervals lead to
improved coverage. In settings where strong scientific knowledge about the
exact form of the outcome model is absent, the bootstrap approach appears more
reliable than its chief competitors.

Coverage | WLS | WLS Cluster | Bootstrap Bias Corrected
---|---|---|---
Linear DGP | 93.2% | 94.8% | 94.8%
Linear DGP, Correlated Errors | 89.4% | 91.5% | 94.5%
Nonlinear DGP, Correlated Errors | 83.4% | 86.0% | 89.8%
Table 1: Coverage of the WLS, WLS cluster and bootstrap bias corrected methods of inference for our three simulation set-ups. Average CI Length | WLS | WLS Cluster | Bootstrap Bias Corrected
---|---|---|---
Linear DGP | 0.25 | 0.27 | 0.27
Linear DGP, Correlated Errors | 0.25 | 0.27 | 0.30
Nonlinear DGP, Correlated Errors | 0.26 | 0.28 | 0.31
Table 2: Average 95% confidence interval length for the WLS, WLS cluster and
bootstrap bias corrected methods of inference for our three simulation set-
ups.

Results in tables 1 and 2 are for GroupMatch with instance replacement,
however matching with trajectory replacement performed very similarly in our
simulations. Computation time was similar for GroupMatch with instance
replacement and with trajectory replacement. In principle, GroupMatch with
instance replacement should be substantially faster, however in its current
form GroupMatch does not implement the most computationally efficient
algorithm for with instance replacement. Over 100 iterations, the average
matching computation time was 4.63 seconds for matching with instance
replacement and 4.71 seconds for matching with trajectory replacement. The
average block bootstrap computation time was 2.51 seconds. Computation time
was calculated on an Apple M1 Max 10-core CPU with 3.22 GHz processor and 64
GB RAM running on macOS Monterey.

##  6 Testing for Timepoint Agnosticism

The key advantage of GroupMatch relative to other matching techniques designed
for rolling enrollment settings is its ability to consider and optimize over
matches between units at different timepoints, which leads to higher quality
matches on lagged covariates. This advantage comes with a price in additional
assumptions, notably the assumption of timepoint agnosticism. Timepoint
agnosticism means that mean potential outcomes under control for any two
individual timepoints in the data should be identical; in particular, this
rules out time trends of any kind in the outcome model that cannot be
explained by covariates in the prior $L$ timepoints.

While in many applications scientific intuition about the data generating
process suggests this assumption may be reasonable, it is essential that we
consider any information contained in the observed data about whether it holds
in a particular case. Accordingly, we present a falsification test for
timepoint agnosticism. Falsification tests are tests “for treatment effects in
places where the analyst knows they should not exist,” [Keele, 2015] and are
useful in a variety of settings in observational studies [Rosenbaum, 1999]. In
particular, our test is designed to detect violations of timepoint
agnosticism, or “treatment effects of time” when they should be absent;
rejections indicate settings in which GroupMatch is not advisable and other
rolling enrollment matching techniques that do not rely on timepoint
agnosticism are likely more suitable. While failure to reject may not
constitute proof positive of timepoint agnosticism’s validity, it rules out
gross violations, thereby limiting the potential for bias.

To test the timepoint agnosticism assumption we use _control-control time
matching_ : matching control units at different timepoints and testing if the
average difference in outcomes between the two timepoint groups, conditional
on relevant covariates, is significantly different from zero using a
permutation test. Specifically, restricting attention to trajectories $i$ from
the control group, we select two timepoints $t_{0}$ and $t_{1}$ and match each
instance at one timepoint to one at the other timepoint using the GroupMatch
optimization routine, based on similarity of covariate histories over the
previous $L$ timepoints. Since this match compares instances at two fixed time
points, any optimal method of matching without replacement may be used. One
practical issue arises: GroupMatch and related matching routines expect one
group to be designated “treated,” all members of which are generally retained
in the match, and the other “control,” some members of which will be included,
but both matching groups are controls in this case. We label whichever of the
two groups has fewer instances as treated; without loss of generality, we will
assume there are fewer instances at time $t_{1}$ and use these instances as
the reference group to be retained.

The test statistic for the falsification test is motivated by the ATT
estimator in section 2.2. Let $N_{c}$ be the total number of control units and
let $N_{t_{1}}$ be the number of control instances at time $t_{1}$. Let
$\hat{\mu}^{t_{0}}_{0}$ be a bias correction model fit on our new control
group (i.e., control instances at time $t_{0}$). In addition, let
$D^{\prime}_{i}=1$ if unit $i$ is present at time $t_{1}$. We define the test
statistic as follows:

$\hat{\Delta}_{cc}=\frac{1}{N_{t_{1}}}\sum_{i=1}^{N_{c}}D^{\prime}_{i}((Y_{i,t=t_{1}}-\hat{\mu}_{0}^{t_{0}}(\mathbf{X}_{i,t=t_{1}}))-\sum_{j=1}^{N_{c}}M_{it_{1},jt_{0}}(Y_{j,t=t_{0}}-\hat{\mu}_{0}^{t_{0}}(\mathbf{X}_{j,t=t_{0}})))$

We use a permutation test to test the following null hypothesis, where
$E_{0}^{t_{1}}\left\\{\cdot\right\\}$ indicates expectation over the
distribution of the covariates in control instances at time $t_{1}$.

$E_{0}^{t_{1}}\left\\{\mu_{0}^{t_{0}}(\mathbf{X})\right\\}=E_{0}^{t_{1}}\left\\{\mu_{0}^{t_{1}}(\mathbf{X})\right\\}$

In words, this null hypothesis says that, accounting for differences in the
covariate distribution at times 0 and 1, the difference in the average
outcomes of control instances at the two timepoints is zero. The test
considers the tail probability of the distribution of the following test
statistic $\widehat{\Delta}_{perm}$ under many draws of the random vector $R$,
where $R=(R_{1},...,R_{N_{c}})$ and $R_{i}$ are independent Rademacher random
variables:

$\widehat{\Delta}_{perm}=\frac{1}{N_{t_{1}}}\sum_{i=1}^{N_{c}}R_{i}D^{\prime}_{i}((Y_{i,t=t_{1}}-\hat{\mu}_{0}^{t_{0}}(\mathbf{X}_{i,t=t_{1}}))-\sum_{j=1}^{N_{c}}M_{it_{1},jt_{0}}(Y_{j,t=t_{0}}-\hat{\mu}_{0}^{t_{0}}(\mathbf{X}_{j,t=t_{0}})))$

This permutation test is identical in implementation to the standard test of a
sharp null hypothesis that outcomes are unchanged by group assignment for
matched designs without replacement, discussed in detail in [Rosenbaum, 2002b,
§2-§3]. In steps:

  1. 1.

Randomly partition the set of control trajectories into two groups. Label
control instances from the first group of trajectories at timepoint $t_{1}$
the new “treated” units, and control instances from the second group of
trajectories at timepoint $t_{0}$ the new “control” units.

  2. 2.

Fit a bias correction model on the new control units.

  3. 3.

Match the new treated units to the new control units and calculate the test
statistic.

  4. 4.

Repeat $B$ times:

    1. (a)

Independently across pairs, switch the treatment and control labels with
probability 0.5.

    2. (b)

Calculate $\hat{\Delta}_{cc}$ on the resulting data.

  5. 5.

Calculate the proportion of permutations with an absolute value of the test
statistic greater than or equal to the absolute observed value calculated in
Step 1. This is the $P$-value.

  6. 6.

If the $P$-value is smaller than $\alpha$, reject the null hypothesis of no
difference between groups.

We do not allow the same unit to appear in both the new control and new
treated group, because this would lead to dependence across matches. We employ
1-1 matching because we are testing a weak null hypothesis using a permutation
test. As Wu and Ding [2020] demonstrate, issues can arise when using a
randomization test to test a weak null hypothesis when treatment and control
sample sizes are not equal. However, 1-1 matching allows us to avoid this
issue by balancing the treatment and control sample sizes. We recommend the
use of caliper matching to ensure high quality matches, especially in the case
where all control units are present at both timepoints. Note that it is
important to permute treatment after matching (indeed, conditional on the
matched pairs chosen). If covariates are strongly related to time, then
covariate distributions at the two timepoints may differ substantially;
permuting labels for the units prior to matching as though they were
exchangeable will not preserve these differences and may lead the test to
reject falsely.

We choose to use a permutation test here rather than the bootstrap because the
data split and 1:1 matching ratio ensure matches are independent under the
original sampling model, making for a tractable permutation distribution. If
desired, the bootstrap approach of Section 4 could be applied instead, and we
expect that results would be similar given the fundamental similarity between
bootstrap and randomization inference where both are viable [Romano, 1989].

See Section C of the supplemental appendix for simulations illustrating this
method.

##  7 Application: Baseball Injuries

We study the impact of short-term injury on hitting performance in
observational data from major league baseball (MLB) during 2013-2017.
Quantitative studies of major league hitting performance [Baumer, 2008] and of
injury trends and impact in athletics [Conte et al., 2016] have been performed
repeatedly, but only a few studies so far have evaluated the impact of injury
on position players’ hitting performance. These have focused on specific
injury types, and have not found strong evidence that injury is associated
with a decline in performance [Begly et al., 2018; Frangiamore et al., 2018;
Wasserman et al., 2015].

We use GroupMatch to match baseball players injured at certain times to
similar players at other points in the season that were not injured. We
evaluate whether players see a decline in offensive performance immediately
after their return from injury. In contrast to other studies, we pool across
injury types to see if there is a more general effect of short term injury on
hitter performance.

###  7.1 Data and Methodology

We use publicly-available MLB player data from Retrosheet.org and injury data
scraped from ProSportsTransactions.com for the years 2013-2017. Our dataset is
composed of player height, weight and age, quantities that remain constant
over a single season of play, as well as on-base percentage (OBP), plate
appearances (PAs) at different points in the season, and dates of short-term
injuries, in which the player’s team designated him for a 7-10 day stay on the
team’s official injured list, for each year. OBP is a common measure of hitter
performance and is approximately equal to the number of times a player reaches
base divided by their number of plate appearances.111OBP = (Hits + Walks + Hit
By Pitch) / (At Bats + Walks + Hit by Pitch + Sacrifice Flies)

For each non-injured player, we generate three pseudo-injury dates evenly
spaced over their PAs. In each season, we match injured players to four non-
injured players. Matches were formed using GroupMatch with instance
replacement, matching on age, weight, height, number of times previously
injured, recent performance measured by OBP over the previous 100 PAs, and
performance over the entire previous year as measured by end-of-year OBP after
James-Stein shrinkage222See https://chris-said.io/2017/05/03/empirical-bayes-
for-multiple-sample-sizes/ for discussion of James-Stein shrinkage to
estimators with variable sample sizes. We choose to shrink the OBP using
James-Stein to limit the impact of sampling variability for players with a
relatively small number of PAs the previous season [Efron and Morris, 1975].

Table 3 shows the balance for each of the covariates before and after
matching. For each covariate, matching shrinks the standardized difference
between the treated and control means.

| Treated | Control Mean | Standardized Difference
---|---|---|---
Variable | Mean | Before | After | Before | After
Height | 73.7 | 73.1 | 73.4 | 0.26 | 0.14
Weight | 213 | 209 | 212 | 0.24 | 0.07
2016 OBP (JS Shrunk) | .324 | .328 | .323 | -0.09 | 0.02
Lag OBP | .336 | .341 | .338 | -0.07 | -0.02
Birth Year | 1988 | 1988 | 1988 | -0.08 | -0.06
Number Previous Injuries | 2.73 | 1.91 | 2.16 | 0.30 | 0.21
Table 3: Balance table for MLB injury analysis before and after matching each
injured player to four non-injured players.

###  7.2 Results

We compare the results for bias-corrected block bootstrap inference, WLS, and
WLS with clustered standard errors. The ATT estimates are positive (0.010),
but the 95% confidence intervals cover zero for all methods, indicating that
there is not strong evidence that short term injury impacts batter
performance. We present the results for 2017 in Figure 2. Results from each of
2013 - 2016 were substantively the same, as were results obtained by pooling
the matched data across years. The data pass the timepoint agnosticism test,
comparing the first and last pseudo-injury dates.

Figure 2: Estimates and 95% confidence intervals for block bootstrap, WLS and
cluster WLS inference methods for the ATT in our 2017 baseball injury
analysis.

##  8 Discussion

The introduction of GroupMatch with instance replacement, a method for block
bootstrap inference, and a test for timepoint agnosticism provide substantial
new capabilities for matching in settings with rolling enrollment. We now
discuss a number of limitations and opportunities for improvement.

Our proof of the block bootstrap approach assumes the use of GroupMatch with
instance replacement. The large-sample properties of matched-pair
discrepancies are substantially easier to analyze mathematically in this
setting than GroupMatch with trajectory replacement or GroupMatch without
replacement, designs in which different treated units may compete for the same
control units, and the technical argument must be altered to account for this
complexity. However, Abadie and Imbens [2012] successfully characterized
similar large-sample properties in cross-sectional settings for matching
without replacement. While beyond the scope of our work here, we believe it is
likely that this approach could provide an avenue for extending Theorem 1 to
cover the other two GroupMatch designs. Empirically, we have found that the
block bootstrap performs well when matches are calculated using any of the
three GroupMatch designs.

Setting aside the technical barriers associated with extending the theory to
GroupMatch without replacement, our new approach provides a competitor method
to the existing randomization inference framework described by Pimentel et al.
[2020] available for GroupMatch without replacement. The randomization
inference framework offers the advantage of closely-related methods of
sensitivity analysis and freedom from making assumptions about the sampling
distribution of the response variables; on the other hand, the block bootstrap
method avoids the need to assume a sharp null hypothesis. In general these
same considerations arise in choosing between sampling-based inference and
randomization-based inference for a cross-sectional matched study, although
such choices have received surprisingly little direct and practical attention
in the literature thus far.

A key consideration for the falsification test proposed in Section 6 is which
timepoints to choose as $t_{0}$ and $t_{1}$. The choice of timepoint
comparison depends largely on what a plausible time trend would be for the
problem at hand. For example, if you suspect a linear time trend, it makes
sense to look at the first and last timepoints. If the trend is linear, this
test should have high power to detect a problem in moderate to large samples.
If one is uncertain about the specific shape of the time trend that is most
likely to occur and wants to test for all possible trends, we recommend
testing each sequential pair of timepoints (i.e., timepoints 1 and 2, 2 and 3,
3 and 4, and so on) and combining the tests via a nonparametric combination of
tests [Pesarin and Salmaso, 2010].

The falsification test is subject to several common criticisms levied at
falsification tests, particularly their ineffectiveness in settings with low
power. One possible approach is to reconfigure the test to assume violation of
timepoint agnosticism as a null hypothesis and seek evidence in the data to
reject it; Hartman and Hidalgo [2018] recommend a similar change for
falsification tests used to assess covariate balance. However, even in the
absence of such a change the test may prove useful in concert with a
sensitivity analysis. Sensitivity analysis, already widely studied in causal
inference as a way to assess the role of ignorability assumptions, places a
nonzero bound on the degree of violation of an assumption and reinterprets the
study’s results under this bound, often repeating the process for larger and
larger values of the bound to gain insight. Such a procedure, which focuses
primarily on assessing the impact of small or bounded violations of an
assumption, naturally complements our falsification test, which can
successfully rule out large violations but is more equivocal about minor
violations.

Unfortunately, no sensitivity analysis appropriate for block bootstrap
inference has been developed yet, either for timepoint agnosticism or other
strong assumptions such as ignorability. The many existing methods for
sensitivity analysis (developed primarily with ignorability assumptions in
mind) are unsatisfying in our framework for a variety of reasons: some rely on
randomization inference [Rosenbaum, 2002b], others focus on weighting methods
rather than matching [Zhao et al., 2019; Soriano et al., 2021], and others are
limited to specific outcome measures [Ding and VanderWeele, 2016] or specific
test statistics [Cinelli and Hazlett, 2020]. We view the development of
compelling sensitivity analysis approaches to be an especially important
methodological objective for matching under rolling enrollment.

## Acknowledgments

The major league baseball performance data was obtained from and is
copyrighted by Retrosheet (www.retrosheet.org). We thank the author and
maintainer of the GitHub repository https://github.com/robotallie/baseball-
injuries for making the injury data easily available. We also thank Eli Ben-
Michael, Peng Ding, Avi Feller, Lauren Forrow, Shirshendu Ganguly, and Jiaqi
Li for helpful conversations and feedback.

## Authors’ Statements

Amanda Glazer acknowledges support from the National Science Foundation (DMS
RTG #1745640).

## References

  * Abadie and Imbens [2006] Abadie, A. and Imbens, G. W. (2006), “Large sample properties of matching estimators for average treatment effects,” Econometrica, 74, 235–267.
  * Abadie and Imbens [2008] — (2008), “On the failure of the bootstrap for matching estimators,” Econometrica, 76, 1537–1557.
  * Abadie and Imbens [2011] — (2011), “Bias-corrected matching estimators for average treatment effects,” Journal of Business & Economic Statistics, 29, 1 – 11.
  * Abadie and Imbens [2012] — (2012), “A Martingale representation for matching estimators,” Journal of the American Statistical Association, 107, 833 – 843.
  * Abadie and Spiess [2021] Abadie, A. and Spiess, J. (2021), “Robust post-matching inference,” Journal of the American Statistical Association, 1–13.
  * Acemoglu et al. [2019] Acemoglu, D., Naidu, S., Restrepo, P., and Robinson, J. A. (2019), “Democracy does cause growth,” Journal of Political Economy, 127, 47–100.
  * Antonelli et al. [2018] Antonelli, J., Cefalu, M., Palmer, N., and Agniel, D. (2018), “Doubly robust matching estimators for high dimensional confounding adjustment,” Biometrics, 74, 1171–1179.
  * Austin and Small [2014] Austin, P. C. and Small, D. S. (2014), “The use of bootstrapping when using propensity-score matching without replacement: a simulation study,” Statistics in medicine, 33, 4306–4319.
  * Baumer [2008] Baumer, B. S. (2008), “Why on-base percentage is a better indicator of future performance than batting average: An algebraic proof,” Journal of Quantitative Analysis in Sports, 4.
  * Begly et al. [2018] Begly, J. P., Guss, M. S., Wolfson, T. S., Mahure, S. A., Rokito, A. S., and Jazrawi, L. M. (2018), “Performance outcomes after medial ulnar collateral ligament reconstruction in Major League Baseball positional players,” Journal of Shoulder and Elbow Surgery, 27, 282–290.
  * Ben-Michael et al. [2021] Ben-Michael, E., Feller, A., and Rothstein, J. (2021), “Synthetic controls and weighted event studies with staggered adoption,” arXiv preprint arXiv:1912.03290.
  * Bohl et al. [2010] Bohl, A. A., Fishman, P. A., Ciol, M. A., Williams, B., LoGerfo, J., and Phelan, E. A. (2010), “A longitudinal analysis of total 3-year healthcare costs for older adults who experience a fall requiring medical care,” Journal of the American Geriatrics Society, 58, 853–860.
  * Cinelli and Hazlett [2020] Cinelli, C. and Hazlett, C. (2020), “Making sense of sensitivity: Extending omitted variable bias,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82, 39–67.
  * Conte et al. [2016] Conte, S., Camp, C. L., and Dines, J. S. (2016), “Injury trends in Major League Baseball over 18 seasons: 1998-2015,” Am J Orthop, 45, 116–123.
  * Ding and VanderWeele [2016] Ding, P. and VanderWeele, T. J. (2016), “Sensitivity analysis without assumptions,” Epidemiology (Cambridge, Mass.), 27, 368.
  * Efron and Morris [1975] Efron, B. and Morris, C. (1975), “Data analysis using Stein’s estimator and its generalizations,” Journal of the American Statistical Association, 70, 311–319.
  * Fogarty [2020] Fogarty, C. B. (2020), “Studentized sensitivity analysis for the sample average treatment effect in paired observational studies,” Journal of the American Statistical Association, 115, 1518–1530.
  * Frangiamore et al. [2018] Frangiamore, S. J., Mannava, S., Briggs, K. K., McNamara, S., and Philippon, M. J. (2018), “Career Length and Performance Among Professional Baseball Players Returning to Play After Hip Arthroscopy,” The American Journal of Sports Medicine, 46, 2588–2593.
  * Hansen [2004] Hansen, B. B. (2004), “Full matching in an observational study of coaching for the SAT,” Journal of the American Statistical Association, 99, 609–618.
  * Hartman and Hidalgo [2018] Hartman, E. and Hidalgo, F. D. (2018), “An Equivalence Approach to Balance and Placebo Tests,” American Journal of Political Science, 62, 1000–1013.
  * Haviland et al. [2008] Haviland, A., Nagin, D. S., Rosenbaum, P. R., and Tremblay, R. E. (2008), “Combining group-based trajectory modeling and propensity score matching for causal inferences in nonexperimental longitudinal data,” Developmental psychology, 44, 422–436.
  * Ho et al. [2007] Ho, D. E., Imai, K., King, G., and Stuart, E. A. (2007), “Matching as nonparametric preprocessing for reducing model dependence in parametric causal inference,” Political analysis, 15, 199–236.
  * Imai et al. [2020] Imai, K., Kim, I. S., and Wang, E. (2020), “Matching methods for causal inference with time-series cross-section data,” Tech. rep., Harvard University.
  * Imbens [2011] Imbens, G. W. (2011), “Experimental design for unit and cluster randomid trials,” in Conference International Initiative for Impact Evaluation, Cuernavaca.
  * Keele [2015] Keele, L. (2015), “The statistics of causal inference: A view from political methodology,” Political Analysis, 313–335.
  * Li and Ding [2017] Li, X. and Ding, P. (2017), “General forms of finite population central limit theorems with applications to causal inference,” Journal of the American Statistical Association, 112, 1759–1769.
  * Li et al. [2001] Li, Y. P., Propert, K. J., and Rosenbaum, P. R. (2001), “Balanced risk set matching,” Journal of the American Statistical Association, 96, 870–882.
  * Lu [2005] Lu, B. (2005), “Propensity score matching with time-dependent covariates,” Biometrics, 61, 721–728.
  * Otsu and Rai [2017] Otsu, T. and Rai, Y. (2017), “Bootstrap inference of matching estimators for average treatment effects,” Journal of the American Statistical Association, 112, 1720–1732.
  * Pesarin and Salmaso [2010] Pesarin, F. and Salmaso, L. (2010), Permutation Tests for Complex Data: Theory, Applications and Software, Wiley.
  * Pimentel et al. [2020] Pimentel, S. D., Forrow, L. V., Gellar, J., and Li, J. (2020), “Optimal matching approaches in health policy evaluations under rolling enrollment,” Journal of the Royal Statistical Society: Series A (Statistics in Society), 183, 1411–1435.
  * Romano [1989] Romano, J. P. (1989), “Bootstrap and randomization tests of some nonparametric hypotheses,” The Annals of Statistics, 141–159.
  * Rosenbaum [1999] Rosenbaum, P. R. (1999), “Choice as an alternative to control in observational studies,” Statistical science, 259–278.
  * Rosenbaum [2002a] — (2002a), “Covariance adjustment in randomized experiments and observational studies,” Statistical Science, 17, 286–327.
  * Rosenbaum [2002b] — (2002b), Observational Studies, New York, NY: Springer.
  * Rubin [1979] Rubin, D. B. (1979), “Using Multivariate Matched Sampling and Regression Adjustment to Control Bias in Observational Studies,” Journal of the American Statistical Association, 318–328.
  * Soriano et al. [2021] Soriano, D., Ben-Michael, E., Bickel, P. J., Feller, A., and Pimentel, S. D. (2021), “Interpretable sensitivity analysis for balancing weights,” arXiv preprint arXiv:2102.13218.
  * Stuart [2010] Stuart, E. (2010), “Matching methods for causal inference: A review and a look forward,” Statistical science: a review journal of the Institute of Mathematical Statistics, 25.
  * Stuart et al. [2011] Stuart, E. A., King, G., Imai, K., and Ho, D. (2011), “MatchIt: nonparametric preprocessing for parametric causal inference,” Journal of statistical software.
  * Wasserman et al. [2015] Wasserman, E. B., Abar, B., Shah, M. N., Wasserman, D., and Bazarian, J. J. (2015), “Concussions are associated with decreased batting performance among Major League Baseball players,” The American Journal of Sports Medicine, 43, 1127–1133.
  * Witman et al. [2019] Witman, A., Beadles, C., Liu, Y., Larsen, A., Kafali, N., Gandhi, S., Amico, P., and Hoerger, T. (2019), “Comparison group selection in the presence of rolling entry for health services research: Rolling entry matching,” Health services research, 54, 492–501.
  * Wu and Ding [2020] Wu, J. and Ding, P. (2020), “Randomization tests for weak null hypotheses in randomized experiments,” Journal of the American Statistical Association, 1–16.
  * Zhao et al. [2019] Zhao, Q., Small, D. S., and Bhattacharya, B. B. (2019), “Sensitivity analysis for inverse probability weighting estimators via the percentile bootstrap,” Journal of the Royal Statistical Society: Series B (Statistical Methodology).

## Supplement to “Robust inference for matching under rolling enrollment”

This appendix provides supplemental material to the manuscript “Robust
inference for matching under rolling enrollment.” We organize this appendix as
follows. In Section A we present the proof of the principle result, Theorem 1,
of the main manuscript. In order, we list the additional assumptions necessary
for the proof, lemmas and their proofs, and finally the proof of the theorem.
Section B discusses the implementation of weighted least squares in common
software packages and its impact on our simulation results. Finally, in
Section C we illustrate the use of the falsification test, described in
Section 6 of the main manuscript, via simulation.

##  Appendix A Proof of Theorem

###  A.1 Assumptions

We begin by rewriting $\sqrt{N_{1}}U^{*}$ as follows:

$\displaystyle\sqrt{N_{1}}U^{*}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}\left(\widehat{\Delta}_{i}^{*}-\widehat{\Delta}_{adj}\right)=\sum_{i=1}^{N}W^{*}_{i}(\hat{\Delta}_{i}-D_{i}\hat{\Delta}_{adj})$

Here the $W_{i}^{*}$, quantities we denote as the bootstrap weights, are
random variables of the form $Q_{i}/N_{1}$ where $Q_{i}$ is a count of the
number of times observation $i$ is selected to appear in the bootstrap sample.
This new form enables us to work separately with stochasticity arising from
the bootstrap and stochasticity arising from the original data-generating
process, and it also makes it easy to generalize our results to other
bootstrap approaches as discussed below.

Our proof relies on the assumptions on sampling described in Section 2.1 and
the GroupMatch identification assumptions described in Section 2.2 of the main
manuscript, with the exception of the exact matching assumption. In addition
we invoke two additional sets of assumptions, which we denote M and R
following similar labeling in Otsu and Rai [2017], from whom we adapt our
proof strategy.

Assumption M. (Conditions for $\hat{\Delta}_{adj}$)

  1. 1.

Let the population distribution of $\mathbf{X}_{i}=(X_{i,1},\ldots,X_{i,T})$
be continuous on $\mathbb{R}^{kT}$ with compact and convex support
$\mathbb{X}^{T}$. In addition, let the densiity of $\mathbf{X}_{i}$ be
bounded, and bounded away from zero on its support.

  2. 2.

For some $r\geq 1$, $\frac{N_{1}^{r}}{N_{0}}\rightarrow\theta$ for
$\theta\in(0,\infty)$.

  3. 3.

For $z=0,1$, let $\mu_{z}^{t}(\mathbf{X})$ be Lipschitz in $\mathbb{X}^{L}$
for all $t=L+1,\ldots,T$.

  4. 4.

For all $t\in\\{L+1,\ldots,T\\}$, and $z=0,1$,
$E[Y_{i,t}^{4}(z)|Z_{i,t-1}=z,\mathbf{X}_{i}]$ and
$Cov(Y_{it},Y_{it^{\prime}}|D_{i}=d,\mathbf{X}_{i,t}=x,\mathbf{X}_{i,t^{\prime}}=x^{\prime})$
are bounded uniformly on $\mathbb{X}^{T}$, and
$\text{Var}(Y_{i,t}(z)|Z_{i,t-1}=z,\mathbf{X}_{i})$ is Lipschitz in
$\mathbb{X}^{T}$ and bounded away from zero.

Assumption R. (Conditions for $\mu_{d}(x)$)

For $d=0,1$ and $\lambda$ satisfying $\sum_{l=1}^{kL}\lambda_{l}={kL}$, the
derivative $\partial^{kL}\mu^{t}_{d}(x)$ exists and satisfies
$\text{sup}_{x\in\mathbb{X}}|\partial^{kL}\mu^{t}_{d}(x)|\leq R$ for some
$R>0$ and for all $t\in\\{L+1,\ldots T\\}$. Furthermore, $\hat{\mu}_{0}(x)$
satisfies
$|\hat{\mu}_{0}(\cdot)-\mu_{0}(\cdot)|_{kL-1}=o_{p}(N^{-1/2+1/(kL)})$.

While not necessary to prove Theorem 1, we mention one additional set of
conditions on the bootstrap weights $W_{i}^{*}$. These are satisfied trivially
by the nonparametric bootstrap we adopt, but in fact the proof goes through
for any bootstrap algorithm that can be represented by a set of bootstrap
weights satisfying these conditions (for instance, the wild bootstrap). Again
following Otsu and Rai [2017], we denote this set of assumptions as Assumption
W and refer to it in our proof to make the path for generalization clear.

Assumption W. (Conditions for $W_{i}^{*}$)

  1. 1.

$(W_{1}^{*},...,W_{N}^{*})$ is exchangeable and independent of
$\mathbf{O}=(\mathbf{Y,D,X})$.

  2. 2.

$\sum_{i=1}^{N}(W_{i}^{*}-\bar{W}^{*})^{2}\xrightarrow[]{p}1$ where
$\bar{W}^{*}=\frac{1}{N}\sum_{i=1}^{N}W_{i}^{*}$

  3. 3.

$\text{max}_{i=1,...,N}|W_{i}^{*}-\bar{W}^{*}|\xrightarrow[]{p}0$

  4. 4.

$E[W_{i}^{*2}]=O(N^{-1})$ for all $i=1,...,N$

###  A.2 Lemmas

To bound the size of matching discrepancies, we compare those obtained by
GroupMatch with instance replacement to nearest-neighbor matching at a fixed
timepoint, in which only one instance from each control unit can potentially
be used in a match. Nearest-neighbor matching matches each treated unit to the
control instance that is most similar. While GroupMatch also uses nearest-
neighbor matching, nearest-neighbor matching at a fixed timepoint considers a
smaller pool of control instances since there is only one instance for each
control unit that can potentially be used in a match. Comparing matching
discrepancies between GroupMatch with instance replacement and nearest-
neighbor matching is important to apply Lemma A.2, used to prove Theorem 2, in
Abadie and Imbens [2011]. For each treated unit $i$, let $j_{m}(i)$ and
$j_{m}^{gm}(i)$ represent the $m$th instance used as a match for nearest-
neighbors at a fixed timepoint and GroupMatch with instance replacement
respectively. Let $U_{m,i}=\mathbf{X}_{j_{m}(i)}-\mathbf{X}_{i,T_{i}}$ and
$U^{gm}_{m,i}=\mathbf{X}_{j^{gm}_{m}(i)}-\mathbf{X}_{i,T_{i}}$ be the matching
discrepancies under nearest neighbors at a fixed timepoint and GroupMatch with
instance replacement matching respectively.

Lemma 1. $||U^{gm}_{m,i}||\leq||U_{m,i}||$ for all $m,i$.

Proof. Because nearest neighbors (NN) matching at a fixed timepoint only
considers one instance from each control trajectory whereas GroupMatch with
instance replacement (GM) considers multiple, the set of control instances
that can be used in a match in NN matching, $\mathcal{C}$, is a subset of the
set of control instances that can be used in a match in GM,
$\mathcal{C}^{gm}$: $\mathcal{C}\subseteq\mathcal{C}^{gm}$. Both NN and GM
match treated units to the control instance that minimizes $U_{m,i}$ and
$U^{gm}_{m,i}$, so we have that

$||U^{gm}_{m,i}||=min_{m,j\in\mathcal{C}^{gm}}|\mathbf{X}_{j}-\mathbf{X}_{i,T_{i}}|\leq
min_{m,j\in\mathcal{C}}|\mathbf{X}_{j}-\mathbf{X}_{i,T_{i}}|=||U_{m,i}||$

where $min_{m}$ denotes the $m$th minimum.

Lemma 2. $E[K_{M}(i,t)^{q}]$ is bounded uniformly in $N$.

Proof. This proof follows closely with the proof of Lemma 3 in Abadie and
Imbens [2006] with modifications described below. Let $f^{t}$ be the density
of $\mathbf{X}_{i,t}$ and define $\underline{f}=inf_{x,w,t}f^{t}_{w}(x)$ and
$\overline{f}=sup_{x,w,t}f^{t}_{w}(x)$. We define the catchment area,
$\mathbb{A}_{M}(i,t)$ as the subset of $\mathbb{X}$ such that control unit $i$
at time $t$ is matched to each observation $j$ at time $t^{\prime}$ with
$W_{j}=1-W_{i}$ and $\mathbf{X}_{j,t^{\prime}}\in\mathbb{A}_{M}(i,t)$:

$\mathbb{A}_{M}(i,t)=\\{x|\sum_{l|W_{l}=W_{i},l\neq
i}1\\{(min_{t^{\prime}}||\mathbf{X}_{l,t^{\prime}}-x||)\leq||\mathbf{X}_{i,t}-x||\\}1\\{min_{t^{\prime}}||\mathbf{X}_{i,t^{\prime}}-x||\geq||\mathbf{X}_{i,t}-x||\\}\leq
M\\}$

Ultimately, we want to bound the volume of the catchment area. In order to do
this, we need to bound the probability that the distance to a match exceeds
some value. To derive this bound, we must account for our trajectory structure
by showing the following inequality:

$\displaystyle
Pr(||\mathbf{X}_{j,t}-\mathbf{X}_{i,t^{\prime}}||>uN_{1-W_{i}}^{-1/k}|W_{1},...,W_{N},\mathbf{X}_{i,t^{\prime}}=x,j\in\mathcal{J}_{M}(i))$
$\displaystyle\leq
Pr(||\mathbf{X}_{j,t}-\mathbf{X}_{i,t^{\prime}}||>uN_{1-W_{i}}^{-1/k}|W_{1},...,W_{N},\mathbf{X}_{i,t^{\prime}}=x,(j,t)=j_{M}(i,t^{\prime}))$
$\displaystyle=\sum_{m=0}^{M-1}\binom{N_{1-W_{i}}}{m}Pr(min_{t}||\mathbf{X}_{j,t}-\mathbf{X}_{i,t^{\prime}}||>uN_{1-W_{i}}^{-1/k}|W_{1},...,W_{N},W_{j}=1-W_{i},\mathbf{X}_{i,t^{\prime}}=x)^{N_{1-W_{i}}-m}\times$
$\displaystyle Pr(min_{t}||\mathbf{X}_{j,t}-\mathbf{X}_{i,t^{\prime}}||\leq
uN_{1-W_{i}}^{-1/k}|W_{1},...,W_{N},W_{j}=1-W_{i},\mathbf{X}_{i,t^{\prime}}=x)^{m}$
$\displaystyle\leq\sum_{m=0}^{M-1}\binom{N_{1-W_{i}}}{m}Pr(||\mathbf{X}_{j,t}-\mathbf{X}_{i,t^{\prime}}||>uN_{1-W_{i}}^{-1/k}|W_{1},...,W_{N},W_{j}=1-W_{i},\mathbf{X}_{i,t^{\prime}}=x)^{N_{1-W_{i}}-m}\times$
$\displaystyle Pr(||\mathbf{X}_{j,t}-\mathbf{X}_{i,t^{\prime}}||\leq
uN_{1-W_{i}}^{-1/k}|W_{1},...,W_{N},W_{j}=1-W_{i},\mathbf{X}_{i,t^{\prime}}=x)^{m}$

The second inequality follows from the fact that the probability that the
minimal distance over a trajectory is less than or equal to any particular
instance. The rest of the proof follows directly from the proof of Abadie and
Imbens [2006]’s Lemma 3 after substituting in our catchment area, this
inequality, and indexing over time.

Lemma 3. Given that $E[K_{M}(i,t)^{q}]$ is uniformly bounded,
$\sum_{i}\sum_{t}K_{M}(i,t)=N_{1}$, and $K_{M}(i,t)\geq 0$, we have that
$\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime})]\leq
cN_{1}^{1+o(1)}$ for some constant $c$.

Proof. For all $\epsilon>0$, we have that

$\displaystyle\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime})]$
$\displaystyle=\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t),K_{M}(i,t^{\prime}),K_{M}(i,t^{\prime\prime}),K_{M}(i,t^{\prime\prime\prime})\leq
N_{1}^{\epsilon}]+$
$\displaystyle\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});max_{t_{1}=t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}}K_{M}(i,t_{1})>N_{1}^{\epsilon}]$
$\displaystyle\leq\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t),K_{M}(i,t^{\prime}),K_{M}(i,t^{\prime\prime}),K_{M}(i,t^{\prime\prime\prime})\leq
N_{1}^{\epsilon}]+$ $\displaystyle
4\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t_{1})>N_{1}^{\epsilon}]$

This upper bound is derived from the fact that

$\displaystyle
E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});max_{t_{1}=t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}}K_{M}(i,t_{1})>N_{1}^{\epsilon}]$
$\displaystyle\leq\sum_{t_{1}\in\\{t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}\\}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t_{1})>N_{1}^{\epsilon}]$

and since we are summing over all
$t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}$ we are able to replace
the summation above with multiplying by four. We bound each of the terms
separately. First,

$\displaystyle\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t),K_{M}(i,t^{\prime}),K_{M}(i,t^{\prime\prime}),K_{M}(i,t^{\prime\prime\prime})\leq
N_{1}^{\epsilon}]$ $\displaystyle\leq
N_{1}^{3\epsilon}\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)]$
$\displaystyle=T^{3}N_{1}^{3\epsilon+1}$

For the next term we use proof by contradiction to show

$\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]\leq
N_{1}.$

Suppose
$\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]>N_{1}$,
then:

$\displaystyle\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)^{\frac{r}{\epsilon}+1}K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]$
$\displaystyle\geq(N_{1}^{\epsilon})^{\frac{r}{\epsilon}}\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]$
$\displaystyle=N_{1}^{r}\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]$
$\displaystyle>N_{1}^{r}N_{1}=N_{1}^{r+1}$ $\displaystyle\implies
NT^{4}sup_{i,t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}}E[K_{M}(i,t)^{\frac{r}{\epsilon}+1}K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]>N_{1}^{r+1}$
$\displaystyle\implies
sup_{i,t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}}E[K_{M}(i,t)^{\frac{r}{\epsilon}+1}K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]>\frac{N_{1}^{r+1}}{NT^{4}}=cN_{1}$

for some constant $c$. This follows from Assumption M.

So, then we have that:

$\displaystyle sup_{i,t}E[K_{M}(i,t)^{\frac{r}{\epsilon}+4}]$
$\displaystyle\geq
sup_{i,t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}}E[K_{M}(i,t)^{\frac{r}{\epsilon}+1}K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime})]$
$\displaystyle>sup_{i,t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}}E[K_{M}(i,t)^{\frac{r}{\epsilon}+1}K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]$
$\displaystyle>cN_{1}$

The first inequality holds by applying Cauchy-Schwarz twice. But then
$E[K_{M}(i,t)^{q}]$ is not uniformly bounded for all $q$, so by contradiction

$\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]\leq
N_{1}.$

So we have that:

$\displaystyle\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime})]$
$\displaystyle=\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});$
$\displaystyle
K_{M}(i,t),K_{M}(i,t^{\prime}),K_{M}(i,t^{\prime\prime}),K_{M}(i,t^{\prime\prime\prime})\leq
N_{1}^{\epsilon}]+$
$\displaystyle\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});max_{t_{1}=t,t^{\prime},t^{\prime\prime},t^{\prime\prime\prime}}K_{M}(i,t_{1})>N_{1}^{\epsilon}]$
$\displaystyle\leq\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});$
$\displaystyle
K_{M}(i,t),K_{M}(i,t^{\prime}),K_{M}(i,t^{\prime\prime}),K_{M}(i,t^{\prime\prime\prime})\leq
N_{1}^{\epsilon}]+$ $\displaystyle
4\sum_{i}\sum_{t}\sum_{t^{\prime}}\sum_{t^{\prime\prime}}\sum_{t^{\prime\prime\prime}}E[K_{M}(i,t)K_{M}(i,t^{\prime})K_{M}(i,t^{\prime\prime})K_{M}(i,t^{\prime\prime\prime});K_{M}(i,t)>N_{1}^{\epsilon}]$
$\displaystyle\leq T^{3}N_{1}^{3\epsilon+1}+4c^{\prime}N_{1}\leq
cN_{1}^{1+o(1)}$

for some constant $c$.

###  A.3 Proof

We follow the proof of Theorem 1 in Otsu and Rai [2017] closely, adapting
where necessary to address the possible presence of multiple correlated
potential outcomes from the same trajectory multiple control instances in the
matched design. In the case where every control unit has only one instance,
our argument reduces to exactly that presented in Otsu and Rai [2017].

We decompose $\sqrt{N_{1}}U^{*}$ as follows:

$\displaystyle\sqrt{N_{1}}U^{*}$
$\displaystyle=\sum_{i=1}^{N}W^{*}_{i}(\hat{\Delta}_{i}-D_{i}\hat{\Delta}_{adj})$
$\displaystyle=\sum_{i=1}^{N}(W^{*}_{i}-\bar{W}^{*})(\hat{\Delta}_{i}-D_{i}\hat{\Delta}_{adj})$
$\displaystyle=\sum_{i=1}^{N}(W^{*}_{i}-\bar{W}^{*})(D_{i}(\hat{\Delta}_{i}-\hat{\Delta}_{adj})+(1-D_{i})\hat{\Delta}_{i})$
$\displaystyle=\sum_{i=1}^{N}(W^{*}_{i}-\bar{W}^{*})[D_{i}(Y_{i,T_{i}}-\hat{\mu}_{0}(\mathbf{X}_{i,T_{i}})-\hat{\Delta}_{adj})+(1-D_{i})\sum_{t=1}^{T}\frac{K_{M}(i,t)}{C}(Y_{i,t}-\hat{\mu}_{0}(\mathbf{X}_{i,t}))]$
$\displaystyle=\sqrt{N_{1}}(T^{*}+R^{*}_{1N_{1}}+R^{*}_{2N_{1}})$

We define the following:

$\displaystyle e_{i,t}$ $\displaystyle=Y_{i,t}-\mu_{D_{i}}(\mathbf{X}_{i,t})$
$\displaystyle\xi_{i,t}$
$\displaystyle=(2D_{i}-1)(\mu_{D_{i}}(\mathbf{X}_{i,t})-\mu_{1-D_{i}}(\mathbf{X}_{i,t}))-\Delta$

We can now rewrite the three components as follows.

$\displaystyle\sqrt{N_{1}}T^{*}$
$\displaystyle=\sum_{i=1}^{N}(W^{*}_{i}-\bar{W^{*}})(D_{i}(e_{i,t=T_{i}}+\xi_{i,t=T_{i}})-(1-D_{i})\sum_{t=1}^{T}\frac{K_{M}(i,t)}{C}e_{i,t})$
$\displaystyle=\sum_{i=1}^{N}(W^{*}_{i}-\bar{W^{*}})[D_{i}((Y_{i,t=T_{i}}-\mu_{1}^{T_{i}}(\mathbf{X}_{i,t}))+(\mu_{1}^{T_{i}}(\mathbf{X}_{i,t})-\mu_{0}(\mathbf{X}_{i,t})-\Delta))$
$\displaystyle-(1-D_{i})\sum_{t=1}^{T}\frac{K_{M}(i,t)}{C}(Y_{i,t}-\mu_{0}(\mathbf{X}_{i,t}))]$
$\displaystyle=\sum_{i=1}^{N}(W^{*}_{i}-\bar{W}^{*})[D_{i}(Y_{i,t=T_{i}}-\mu_{0}(\mathbf{X}_{i,t})-\Delta)-(1-D_{i})\sum_{t=1}^{T}\frac{K_{M}(i,t)}{C}(Y_{i,t}-\mu_{0}(\mathbf{X}_{i,t}))]$

$\displaystyle\sqrt{N_{1}}R^{*}_{1N_{1}}$
$\displaystyle=\sum_{i=1}^{N}(W^{*}_{i}-\bar{W^{*}})(D_{i}(\mu_{0}(\mathbf{X}_{i,t})-\hat{\mu}_{0}(\mathbf{X}_{i,t}))-(1-D_{i})\frac{K_{M}(i,t)}{C}\sum_{t=1}^{T}(\mu_{0}(\mathbf{X}_{i,t})-\hat{\mu}_{0}(\mathbf{X}_{i,t})))$

$\displaystyle\sqrt{N_{1}}R^{*}_{2N_{1}}$
$\displaystyle=\sum_{i=1}^{N}(W^{*}_{i}-\bar{W^{*}})D_{i}(\Delta-\hat{\Delta}_{adj})$

We have that
$Pr\\{\sqrt{N_{1}}R^{*}_{1N_{1}}>\epsilon|\mathbf{O}\\}\xrightarrow[]{p}0$ and
$Pr\\{\sqrt{N_{1}}R^{*}_{2N_{1}}>\epsilon|\mathbf{O}\\}\xrightarrow[]{p}0$ for
any $\epsilon>0$, by the same argument as Otsu and Rai [2017] which utilizes
our Assumptions W, R, Lemma 2 and the Markov Inequality. This part of the
argument also relies on Lemma 1, although the reliance is not explicit in Otsu
and Rai [2017]; see Abadie and Imbens [2011] for a similar derivation where
the role of Lemma 1 is more clear.

Next, to show that that $\text{sup}_{r}|Pr\\{\sqrt{N_{1}}T^{*}\leq
r|\mathbf{O}\\}-Pr\\{\sqrt{N_{1}}(\hat{\Delta}_{adj}-\Delta)\leq
r\\}|\xrightarrow[]{p}0$, we define:

$\displaystyle\eta_{i}$
$\displaystyle=[D_{i}(Y_{i,t=T_{i}}-\mu_{0}(\mathbf{X}_{i,t})-\Delta)-(1-D_{i})\sum_{t=1}^{T}\frac{K_{M}(i,t)}{C}(Y_{i,t}-\mu_{0}(\mathbf{X}_{i,t}))]/\sqrt{N_{1}}$
$\displaystyle=D_{i}(e_{i,t=T_{i}}+\xi_{i,t=T_{i}})-(1-D_{i})\sum_{t=1}^{T}\frac{K_{M}(i,t)}{C}e_{i,t}$

We have the following:

$\displaystyle\sigma_{N}^{2}$ $\displaystyle=\sigma_{1N}^{2}+\sigma_{2}^{2}$
$\displaystyle\sigma_{1N}^{2}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}Var(\sum_{t=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))Y_{it}|\mathbf{D},\mathbf{X})$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}Cov\Big{\\{}\sum_{t=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))Y_{it},$
$\displaystyle\sum_{t=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))Y_{it}|\mathbf{D},\mathbf{X}\Big{\\}}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(D_{i}1\\{t^{\prime}=T_{i}\\}+\frac{K_{M}(i,t^{\prime})}{C}(1-D_{i}))Cov(Y_{it},Y_{it^{\prime}})$
$\displaystyle\sigma_{2}^{2}$
$\displaystyle=E[((\mu_{1}^{T_{i}}(\mathbf{X}_{i,t})-\mu_{0}(\mathbf{X}_{i,t}))-\Delta)^{2}|D_{i}=1]$

Within $\sigma_{1N}^{2}$, note that
$Var(\sum_{t=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))Y_{it}|\mathbf{D},\mathbf{X})$
reduces to $Var(Y_{it})$ for treated units. However, for control units, we
have this variance term plus extra covariance terms between that control unit
instance and the other instances in its trajectory. From here, we are able to
follow the same proof strategy as in Otsu and Rai [2017], with our modified
assumptions, and some modifications (detailed below) to Lemmas (i)-(iii) in
Otsu and Rai [2017] to account for the extra covariance terms resulting from
the trajectory structure of the control units. Lemmas (i)-(iii) in Otsu and
Rai [2017] show that the sampling variance of the $\eta_{i}$’s converge to the
population variance, $\sigma_{N}^{2}$, and that other random variables
converge in probability to 0. They show this by leveraging the boundedness of
higher order moments of $\eta_{i}$ and use of the Markov inequality. We are
able to utilize the same arguments, with changes laid out more explicitly for
lemma (i), below, and generalized for lemmas (ii) and (iii).

To apply Lemmas (i)-(iii) in Otsu and Rai [2017] we define:

$\displaystyle\hat{\sigma}_{1N}^{2}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(D_{i}1\\{t^{\prime}=T_{i}\\}+\frac{K_{M}(i,t^{\prime})}{C}(1-D_{i}))e_{i,t}e_{i,t^{\prime}}$

Then note that:

$\displaystyle E[(\hat{\sigma}_{1N}^{2}-\sigma_{1N}^{2})^{2}]$
$\displaystyle=E[\\{(\frac{1}{N_{1}}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(D_{i}1\\{t^{\prime}=T_{i}\\}+\frac{K_{M}(i,t^{\prime})}{C}(1-D_{i}))e_{i,t}e_{i,t^{\prime}})-$
$\displaystyle(\frac{1}{N_{1}}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(D_{i}1\\{t^{\prime}=T_{i}\\}+\frac{K_{M}(i,t^{\prime})}{C}(1-D_{i}))Cov(Y_{it},Y_{it^{\prime}}))\\}^{2}]$
$\displaystyle=\frac{1}{N_{1}^{2}}E[\\{\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(D_{i}1\\{t^{\prime}=T_{i}\\}+\frac{K_{M}(i,t^{\prime})}{C}(1-D_{i}))$
$\displaystyle\times(e_{i,t}e_{i,t^{\prime}}-Cov(Y_{it},Y_{it^{\prime}}))\\}^{2}]$
$\displaystyle=\frac{1}{N_{1}^{2}}E[\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}\sum_{j=1}^{N}\sum_{t_{1}=1}^{T}\sum_{t_{1}^{\prime}=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(D_{i}1\\{t^{\prime}=T_{i}\\}+\frac{K_{M}(i,t^{\prime})}{C}(1-D_{i}))$
$\displaystyle\times(D_{j}1\\{t_{1}=T_{j}\\}+\frac{K_{M}(j,t_{1})}{C}(1-D_{j}))(D_{j}1\\{t_{1}^{\prime}=T_{j}\\}+\frac{K_{M}(j,t_{1}^{\prime})}{C}(1-D_{j}))$
$\displaystyle\times(e_{i,t}e_{i,t^{\prime}}-Cov(Y_{it},Y_{it^{\prime}}))(e_{j,t_{1}}e_{j,t_{1}^{\prime}}-Cov(Y_{jt_{1}},Y_{jt_{1}^{\prime}}))]$
$\displaystyle=\frac{1}{N_{1}^{2}}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}\sum_{t_{1}=1}^{T}\sum_{t_{1}^{\prime}=1}^{T}E[(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(D_{i}1\\{t^{\prime}=T_{i}\\}+\frac{K_{M}(i,t^{\prime})}{C}(1-D_{i}))$
$\displaystyle\times(D_{i}1\\{t_{1}=T_{i}\\}+\frac{K_{M}(i,t_{1})}{C}(1-D_{i}))(D_{i}1\\{t_{1}^{\prime}=T_{i}\\}+\frac{K_{M}(i,t_{1}^{\prime})}{C}(1-D_{i}))$
$\displaystyle\times
E[(e_{i,t}e_{i,t^{\prime}}-Cov(Y_{it},Y_{it^{\prime}}))(e_{i,t_{1}}e_{i,t_{1}^{\prime}}-Cov(Y_{it_{1}},Y_{it_{1}^{\prime}}))|D_{i},\mathbf{X}_{i,t},\mathbf{X}_{i,t^{\prime}},\mathbf{X}_{i,t_{1}},\mathbf{X}_{i,t^{\prime}_{1}}]]$
$\displaystyle\leq\frac{1}{N_{1}^{2}}(N_{1}+cN_{1}^{1+o(1)})\times$
$\displaystyle 3\times max(sup_{d,x}E[e_{it}^{4}|D_{i}=d,\mathbf{X}_{i,t}=x],$
$\displaystyle-2sup_{d,x,x^{\prime}}E[e^{2}_{it}|D_{i}=d,\mathbf{X}_{i,t}=x]sup_{d,x,x^{\prime}}Cov(Y_{it},Y_{it^{\prime}}|D_{i}=d,\mathbf{X}_{i,t}=x,\mathbf{X}_{i,t^{\prime}}=x^{\prime}),$
$\displaystyle
sup_{d,x,x^{\prime}}Cov(Y_{it},Y_{it^{\prime}}|D_{i}=d,\mathbf{X}_{i,t}=x,\mathbf{X}_{i,t^{\prime}}=x^{\prime})^{2})\rightarrow
0$

The inequality follows from Lemma 3. The convergence follows from Assumption
M(4) and Lemma 2.

####  A.3.1 Difference-in-differences Estimator

This proof extends easily to the difference-in-differences ATT estimator
described in Section 4.3 of the main manuscript with the modifications
described below.

We replace $\hat{\Delta}_{adj}$ with $\hat{\Delta}_{DiD}$ and the outcome
$Y_{i,t}$ with the difference-in-differences outcome $Y_{i,t}-Y_{i,t-1}$

We modify

$e_{i,t}=Y_{i,t}-Y_{i,t-1}-(\mu_{D_{i}}(\mathbf{X}_{i,t})-\mu_{D_{i}}(\mathbf{X}_{i,t-1}))$

and

$\xi_{i,t}=(2D_{i}-1)((\mu_{D_{i}}(\mathbf{X}_{i,t})-\mu_{1-D_{i}}(\mathbf{X}_{i,t}))-(\mu_{D_{i}}(\mathbf{X}_{i,t-1})-\mu_{1-D_{i}}(\mathbf{X}_{i,t-1})))$

.

We also update our variance formulas to reflect the additional covariance
terms introduced by the difference-in-differences outcome. In particular, our
formula for $\sigma_{1N}^{2}$ includes additional terms:

$\displaystyle\sigma_{1N}^{2}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}Var(\sum_{t=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(Y_{it}-Y_{i,t-1})|\mathbf{D},\mathbf{X})$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}Cov\Big{\\{}\sum_{t=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(Y_{it}-Y_{i,t-1}),$
$\displaystyle\sum_{t=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(Y_{it}-Y_{i,t-1})|\mathbf{D},\mathbf{X}\Big{\\}}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}(D_{i}1\\{t=T_{i}\\}+\frac{K_{M}(i,t)}{C}(1-D_{i}))(D_{i}1\\{t^{\prime}=T_{i}\\}+\frac{K_{M}(i,t^{\prime})}{C}(1-D_{i}))\times$
$\displaystyle(Cov(Y_{it},Y_{it^{\prime}})-Cov(Y_{it},Y_{i,t^{\prime}-1})-Cov(Y_{i,t-1},Y_{it^{\prime}})+Cov(Y_{i,t-1},Y_{i,t^{\prime}-1}))$

To show that $E[(\hat{\sigma}_{1N}^{2}-\sigma_{1N}^{2})^{2}]\rightarrow 0$ we
use a similar argument to before, where we are able to bound the difference
between the $e_{i,t}e_{i,t^{\prime}}$ and $Cov(Y_{it},Y_{it^{\prime}})$ by a
fourth order polynomial in $e_{it}$ terms and use the assumption that the
expectation of these fourth moments is bounded.

##  Appendix B Weighted Least Squares

Weighted least squares (WLS) is commonly used with matching weights to
calculate the ATT and its corresponding confidence interval [Ho et al., 2007;
Stuart et al., 2011]. For WLS, we have the following estimator for $\beta$:

$\displaystyle\hat{\beta}=(X^{T}WX)^{-1}X^{T}WY$

Where $W$ is a diagonal matrix with entries corresponding to the matching
weights. Note that $W=W^{T}$. We can compute the variance as follows:

$\displaystyle
Var(\hat{\beta})=(X^{T}WX)^{-1}X^{T}W[Var(Y)]((X^{T}WX)^{-1}X^{T}W)^{T}$

Software to compute WLS estimates, such as lm and Zelig (which calls on lm) in
R, often assumes that $Var(Y)=\sigma^{2}W^{-1}$. If this is true we get a
cancellation that yields a nicer formula for variance:

$\displaystyle Var(\hat{\beta})$
$\displaystyle=(X^{T}WX)^{-1}X^{T}W\sigma^{2}W^{-1}((X^{T}WX)^{-1}X^{T}W)^{T}$
$\displaystyle=(X^{T}WX)^{-1}X^{T}WW^{-1}W^{T}X\sigma^{2}(W^{T}WX)^{-1}$
$\displaystyle=(X^{T}WX)^{-1}(X^{T}WX)\sigma^{2}(W^{T}WX)^{-1}$
$\displaystyle=\sigma^{2}(W^{T}WX)^{-1}$

In R, lm (and Zelig) use this formula333Last verified 04-11-2022. in order to
compute standard errors and confidence intervals for WLS regression. However,
when this assumption is not true, as is the case in our simulations (and often
in practice), this formula is incorrect.

Assume $Var(Y)=\sigma^{2}I$, as in some of our simulations, then:

$\displaystyle Var(\hat{\beta})$
$\displaystyle=(X^{T}WX)^{-1}X^{T}W\sigma^{2}I((X^{T}WX)^{-1}X^{T}W)^{T}$
$\displaystyle=\sigma^{2}(X^{T}WX)^{-1}X^{T}W^{2}X(X^{T}WX)^{-1}$

To get an unbiased estimator of $\sigma^{2}$ we note that:

$\displaystyle E[e^{T}Ie]=tr(I\Sigma_{ee})$

Where $e=Y-\hat{Y}=Y-X\hat{\beta}=Y-X(X^{T}WX)^{-1}X^{T}WY$. Now,

$\displaystyle\Sigma_{ee}$
$\displaystyle=Cov(Y-X(X^{T}WX)^{-1}X^{T}WY,Y-X(X^{T}WX)^{-1}X^{T}WY)$
$\displaystyle=Var(Y)-2Cov(Y,X(X^{T}WX)^{-1}X^{T}WY)+Var(X(X^{T}WX)^{-1}X^{T}WY)$
$\displaystyle=\sigma^{2}I-2X(X^{T}WX)^{-1}X^{T}WVar(Y)+(X(X^{T}WX)^{-1}X^{T}W)^{T}(X(X^{T}WX)^{-1}X^{T}W)Var(Y)$
$\displaystyle=\sigma^{2}(I-2X(X^{T}WX)^{-1}X^{T}W+WX(X^{T}WX)^{-1}X^{T}X(X^{T}WX)^{-1}X^{T}W)$

Thus, to get an unbiased estimator of $\sigma^{2}$ we must divide $e^{T}e$ by
the trace of
$I-2X(X^{T}WX)^{-1}X^{T}W+WX(X^{T}WX)^{-1}X^{T}X(X^{T}WX)^{-1}X^{T}W$. When we
use this formula to create confidence intervals for WLS, instead of using lm,
our simulations cover in the linear DGP, uncorrelated errors case. However, we
do not expect that this correction addresses all issues with model-based
parametric standard errors after matching under rolling enrollment — in
particular, many of the problems identified by Abadie and Spiess [2021] likely
persist in some form — and we recommend use of the block bootstrap procedure
discussed in the main manuscript instead.

##  Appendix C Falsification Test Simulations

In this section, we illustrate the timepoint agnosticism falsification test
presented in Section 6 of the main minuscript via simulation. We generate a
dataset of 1000 control units each with 4 covariates and 2 instances occuring
at times $t_{0}$ and $t_{1}$. Two of the covariates vary with time, and two
are uniform across time:

$\displaystyle X_{1,i,t},X_{2,i,t},X_{3,i,t_{0}},X_{4,i,t_{0}}\sim
N(0,1)\text{ for }t=t_{0},t_{1}$ $\displaystyle
X_{j,i,t_{1}}=X_{j,i,t_{0}}+\epsilon_{j,i}\text{ for }j=3,4$
$\displaystyle\epsilon_{j,i}\sim N(0,0.5^{2})\text{ for }j=3,4$

The outcome variable is a linear combination of the four covariates, a time
trend controlled by parameter $\gamma$, and an error term ($\epsilon_{i,t}\sim
N(0,1)$):

$Y_{i,t}=log(4)(X_{1,i,t}+X_{4,i,t})+log(10)(X_{3,i,t}+X_{4,i,t})+1\\{t=t_{1}\\}\gamma+\epsilon_{i,t}$

We generate data for the setting $\gamma=0$, which does not include a time
varying component, $1000$ times. On each dataset, we perform the test for
timepoint agnosticism outlined in this section, with 1-1 matching and bias
correction. Our simulated data only contains two timepoints and the sample
size is balanced so we choose the first timepoint as $t_{0}$, our new control
group, and the second timepoint as $t_{1}$, our new treated group. In 0.049 of
the simulations the $P$-value is less than 0.05, which shows that type I error
is controlled. Now, we add in a time trend where there is an additional term,
$\gamma$, added to the second timepoint. For $\gamma=0.1$ for the second
timepoint (resulting in a time trend of 0.1), our simulations result in a
$P$-value less than 0.05 in 0.327 of the 1000 simulations. For $\gamma=0.25$
for the second timepoint (resulting in a time trend of 0.25), our simulations
result in a $P$-value less than 0.05 in 0.981 of the 1000 simulations. Table 4
summarizes these results.

Figure 3 shows two simulated datasets with $\gamma=0.1$. The test for
timepoint agnosticism detects the trend in one of the two datasets. Overall,
the simulations show that the test is not a panacea for issues with timepoint
agnosticism, sometimes failing to detect small violations. However, it still
adds substantial value to the analysis pipeline, detecting moderate violations
of timepoint agnosticism not especially obvious to the eye in visualization
plots with a very high rate of success.

Time Trend, $\gamma$ | 0 | 0.1 | 0.25
---|---|---|---
Proportion $P$-values $<0.05$ | 0.049 | 0.327 | 0.981
Table 4: Summary of timepoint agnosticism simulation results. Proportion of
1000 simulations where the $P$-value from the test is less than 0.05, for time
trends, $\gamma=0,0.1,0.25$.

Figure 3: Simulated datasets with time trend, $\gamma=0.1$. The figures show
outcome data for a dataset where the timepoint agnosticism test detects the
time trend ($P=0.04$) and does not detect the time trend ($P=0.37$)
respectively.