# Generalized difference-in-differences

Denis Agniel, Max Rubinstein, Jessie Coe, Maria DeYoreo

###### Abstract

We propose a new method for estimating causal effects in longitudinal/panel
data settings that we call generalized difference-in-differences. Our approach
unifies two alternative approaches in these settings: ignorability estimators
(e.g., synthetic controls) and difference-in-differences (DiD) estimators. We
propose a new identifying assumption – a stable bias assumption – which
generalizes the conditional parallel trends assumption in DiD, leading to the
proposed generalized DiD framework. This change gives generalized DiD
estimators the flexibility of ignorability estimators while maintaining the
robustness to unobserved confounding of DiD. We also show how ignorability and
DiD estimators are special cases of generalized DiD. We then propose
influence-function based estimators of the observed data functional, allowing
the use of double/debiased machine learning for estimation. We also show how
generalized DiD easily extends to include clustered treatment assignment and
staggered adoption settings, and we discuss how the framework can facilitate
estimation of other treatment effects beyond the average treatment effect on
the treated. Finally, we provide simulations which show that generalized DiD
outperforms ignorability and DiD estimators when their identifying assumptions
are not met, while being competitive with these special cases when their
identifying assumptions are met.

##  1 Introduction

Researchers often use longitudinal or panel data to assess causal questions in
public policy, economics, health services research, medicine, and other
fields. These studies typically include units (e.g., employees,
municipalities, patients) observed over time (e.g., months, quarters, years),
where some of the units are exposed to an intervention or a treatment during
the course of the study. For example, employees may participate in a training
program, municipalities may enact a policy change, or patients may initiate a
new treatment.

To identify and estimate an average treatment effect with observational, non-
experimental data, two main approaches exist. The first invokes an
ignorability assumption: that the treatment was effectively randomized within
units sharing pre-treatment outcome and covariate values. The ignorability
assumption allows for non-random treatment in that treatment can be related to
all observed pre-treatment variables. The ignorability assumption often
underpins various standard estimation approaches that effectively control for
observed variables in estimation of a treatment effect, such as regression,
propensity score weighting/matching, or synthetic controls. While this
assumption allows for treated and untreated units to systematically differ,
these differences are assumed only to occur among observed variables. In other
words, we assume the researcher observes all confounding variables.

The other main approach is a difference-in-differences (DiD) or event study
design. Instead of relying on an ignorability assumption, DiD designs allow
for some unobserved confounding between treated and untreated units.
Specifically, DiD designs invoke a parallel trends assumption that allows for
average outcomes between treated and untreated units to have different levels,
as long as they would have the same changes over time absent treatment. When
observed variables are thought to influence the evolution of outcomes, DiD
designs hinge on a conditional parallel trends assumption. The conditional
parallel trends assumption states that, within subgroups defined by levels of
the observed variables, the average outcome of the treated units would have
evolved in parallel to the average outcome of the untreated units had the
treated units not received treatment. While the strain of literature related
to estimation of treatment effects under an ignorability assumption has
increasingly developed and adopted sophisticated machine learning techniques
to account for high-dimensional control variables and non-linear functional
forms, the DiD literature has largely relied on time-invariant controls and
using generalized linear models for estimation.

We propose a novel approach to causal identification and estimation of average
treatment effect parameters from observational, longitudinal data that unifies
these two strategies. We call our design “generalized difference-in-
differences” (gDiD). The validity of gDiD relies on a stable bias assumption
which generalizes conditional parallel trends. While the conditional parallel
trends assumption is generally stated as an equivalence between two groups,
the stable bias assumption reformulates the assumption as an equivalence
between two time periods. Conditional parallel trends is generally stated as:
the conditional average difference over time in the outcomes of the treated
units had they not been treated equals the conditional average difference over
time in the outcomes of the untreated units, where the conditioning set is the
same for both groups. The stable bias assumption states that the difference in
conditional average outcomes of the treated units and untreated units after
treatment had no one been treated equals the difference in conditional average
outcomes of the treated units and untreated units before treatment.
Reformulating the assumption as a comparison of time periods instead of groups
clarifies how all pre-period information (for example, pre-period outcomes and
other time-varying variables) can be accounted for in the conditioning. In
particular, controlling for prior outcomes is straightforward. Like
ignorability, the stable bias assumption allows for a rich dependence between
the post-treatment outcomes and pre-period outcomes. Like parallel trends, the
stable bias assumption allows for the presence of unmeasured confounding with
stable effects over time. We further show that conditional and unconditional
DiD estimators are special cases of gDiD, as are ignorability estimators under
additional assumptions.

The stable bias assumption states that the bias of an estimator in the post-
treatment period can be estimated by examining the same estimator in the pre-
treatment period. The identification result thus suggests debiasing an
ignorability estimand in the post-treatment period by subtracting the
corresponding estimand in the pre-treatment period. In this paper, we focus on
estimation of the average treatment effect on the treated (ATT). We propose an
influence-function-based estimation strategy and provide the asymptotic
distribution of the estimator. To implement our estimator, we make use of
modern double/debiased machine learning techniques. In settings where the
usual conditional parallel trends assumption holds, making use of these
nonparametric estimation techniques makes standard DiD estimators robust to
model mis-specification, and our simulation results suggest that this non-
parametric approach remains competitive with DiD when parallel trends holds.
Additionally, we provide simulation evidence that when conditional parallel
trends holds but not all conditioning variables are observed to the
researcher, a setting plausible in practice, our gDiD approach which includes
prior outcomes in the conditioning set has lower bias than usual DiD
estimators.

Our contributions are thus threefold: first, we propose the gDiD design and
show how the stable bias assumption provides a unified identification
framework to think about both DiD estimators and ignorability estimators;
second, we provide simulation evidence showing that the gDiD estimator
performs well under various settings, including where ignorability and
parallel trends hold; and third, our proposed estimation strategy, which can
be implemented for standard DiD estimators as well as our gDiD estimator,
employs the latest best practices in non-parametric machine learning
techniques. Our accompanying R package implements these tools using
SuperLearner to allow for flexible machine learning methods for estimation.

###  1.1 Related work

This work builds on recent advances in the difference-in-differences
literature [23]. Many recent estimators have generalized the classical
parallel trends assumption [8, 25, 6]. These conditional parallel trends
assumptions allow the relationship between trends in the counterfactual
outcome under no intervention to have a complex dependence on covariates. We
further generalize these assumptions to depend on pre-period covariates – in
such a way that they’re no longer “parallel trends” assumptions, but related
“stable bias” assumptions. This allows for more complex relationships between
post-period (counterfactual) outcomes and pre-period outcomes. Furthermore,
this additional conditioning on pre-period outcomes can minimize the effect of
missing a covariate that would otherwise be required to ensure conditional
parallel trends (see Section 4 for simulation evidence in this setting).

We also build on the application of synthetic control techniques [1, 5] to the
setting where many units may be treated. While the assumptions of synthetic
control estimators are often written in terms of an assumed model (e.g., a
linear factor model), they can alternatively be justified by an ignorability
assumption conditional on the pre-treatment outcomes. Augmented synthetic
controls [5] can furthermore be seen as a special case of an augmented inverse
probability weighted estimator [3, 10] where the pre-treatment outcomes are
treated as confounders. Generalized difference-in-differences relaxes the
ignorability assumption these estimators make, while still allowing
practitioners to take advantage of the full suite of ignorability tools –
i.e., synthetic control estimation or double machine learning.

Finally, this work addresses the concerns raised by Daw and Hatfield [13] that
controlling for pre-period outcomes can lead to bias in DiD settings because
of regression to the mean. They specifically consider a caliper matching
strategy, where units in the treated sample that do not have a suitable match
in the control sample are discarded. Related work by Ham and Miratrix [16]
suggests that matching on pre-period outcomes can improve on traditional DiD
and provide conditions under a simple model where this improvement can be
ensured. The proposed gDiD framework explicitly conditions on pre-period
outcomes to achieve causal identification and establishes conditions under
which this is justified. We argue that the required stable bias assumption is
much more flexible than parallel trends (and thus likely in general more
plausible), but we also show how this additional flexibility requires a
stricter positivity (or overlap) assumption. The setting of Daw and Hatfield
is an example of a failure of our positivity assumption, which may be
investigated using standard causal inference tools.

The rest of the paper proceeds as follows. In Section 2, we outline the basic
problem setup and the assumptions we require to identify the target causal
effect, including the stable bias assumption. We compare our identifying
assumptions to the corresponding assumptions (ignorability, parallel trends,
and conditional parallel trends) for alternative estimators. We also give the
influence function for the identifying functional of the ATT under our
assumptions, use it to propose an estimator, and derive the limiting normal
distribution of the estimator. In Section 3, we outline extensions of our
approach to staggered adoption and clustered treatment assignment settings,
and we discuss extensions to other types of treatment effects beyond the ATT.
In Section 4, we give extensive simulation results for our approach in
comparison to competing DiD and ignorability estimators. We apply gDiD to
estimate the causal effect of health system affiliation on health care quality
among Medicare Fee-for-service beneficiaries in Section 5. Finally, we give
concluding remarks in Section 6.

##  2 Generalized difference-in-differences

###  2.1 Preliminaries and basic assumptions

We observe $n$ units and assume ${\bf O}=({\bf
O}_{i})_{i=1,...,n}=(Y_{i1},A_{i},{\bf X}_{i},\bar{{\bf Y}}_{i0})_{i=1,...,n}$
are a set of iid random variables. For simplicity, we consider a single post-
treatment time period $t=1$ and a binary treatment $A_{i}=\\{0,1\\}$ such that
no units are treated before period 1 and some units are treated in period 1.
$Y_{it}$ denotes the outcomes at time $t$, which we group into pre-treatment
outcomes up to time 0, $\bar{{\bf Y}}_{i0}=(Y_{it})_{t\leq 0}$, and post-
treatment outcomes, $Y_{i1}$. The ${\bf X}_{i}$ are a set of time-invariant
covariates.111While our notation implies that the covariates are time-
invariant, ${\bf X}_{i}$ may also include vectors of time-varying covariates
measured prior to treatment, similar to $\bar{{\bf Y}}_{i0}$. Let ${\bf
W}_{it}=({\bf X}_{i},\bar{{\bf Y}}_{it})$ for $t\leq 0$ be the set of
covariates and pre-treatment outcomes at time $t$.

We make use of potential outcomes notation [24] and define the potential
outcome $Y_{i1}(a)$ as the outcome at $t=1$ that one would observe for unit
$i$ if, possibly contrary to fact, the treatment took value $A_{i}=a$ for
$a=\\{0,1\\}$.

Our primary causal estimand of interest is the average treatment effect on the
treated (ATT) at time-period $t=1$:

$\displaystyle\tau=E\\{Y_{i1}(1)-Y_{i1}(0)|A_{i}=1\\}$

This is a conventional estimand in policy evaluation settings where we observe
longitudinal or panel data.222We later extend this to staggered adoption
settings in the vein of [8] in Section 3.1 as well as to clustered treatment
assignment in Section 3.2. We will also later consider other estimands like
the average treatment effect (not limited to the treated group) in Section
3.3.

To identify this estimand from the observed data, we require the following
standard assumptions:

######  Assumption 1 (Consistency).

$Y_{i1}=A_{i}Y_{i1}(1)+(1-A_{i})Y_{i1}(0)$.

######  Assumption 2 (No anticipation).

$Y_{it}(a)=Y_{it},a=0,1;t<1.$

######  Assumption 3 (Positivity).

$P\left\\{\pi({\bf W}_{i0})<1-\epsilon\right\\}=1$ for the propensity score
$\pi({\bf w})=P(A_{i}=1|{\bf W}_{i0}={\bf w})$ and some $\epsilon>0$.

Assumption 1 allows us to connect the observed data with counterfactual
outcomes and states that we observe the counterfactual outcome under treatment
$A_{i}=a$ when $A_{i}$ takes value $a$. Assumption 2 ensures no anticipatory
effects of treatment prior to the intervention and states that counterfactuals
in the pre-treatment period are observed for treatment and control. Finally,
Assumption 3 ensures that there is sufficient overlap in the treated and
control groups and states that there is at least some probability of finding
units like the treated units among the controls.

###  2.2 Stable bias assumption and identification

The crucial identifying assumption for our generalized difference-in-
differences approach is the stable bias assumption:

######  Assumption 4 (Stable bias).

$\displaystyle E\left[E\left\\{Y_{i1}(0)|{\bf
W}_{i0},A_{i}=1\right\\}-E\left\\{Y_{i1}|{\bf
W}_{i0},A_{i}=0\right\\}|A_{i}=1\right]=E\left[E\left\\{Y_{i0}|{\bf
W}_{i-1},A_{i}=1\right\\}-E\left\\{Y_{i0}|{\bf
W}_{i-1},A_{i}=0\right\\}|A_{i}=1\right].$

This assumption states that for counterfactual outcomes under no treatment,
the difference in conditional average outcomes of the treated units and
untreated units at $t=1$ (after treatment) equals the corresponding difference
at $t=0$ (before treatment).

Under Assumptions 1, 2, 3, and 4, $\tau$ may be identified in the observed
data via the following expression:

$\displaystyle\begin{aligned}
\tau^{\text{gDiD}}=&E(Y_{i1}|A_{i}=1)-E\left\\{E(Y_{i1}|{\bf
W}_{i0},A_{i}=0)|A_{i}=1\right\\}-\\\
&\qquad\left\\{E(Y_{i0}|A_{i}=1)-E\left\\{E(Y_{i0}|{\bf
W}_{i-1},A_{i}=0)|A_{i}=1\right\\}\right\\}\end{aligned}$ (1)

This observed data functional has a difference-in-differences structure. The
first differences is the difference in the mean post-treatment observed
outcomes among the treated units and the mean post-treatment outcomes among
the control units, if the control units had a distribution of ${\bf W}_{i0}$
similar to the treated units. The second difference is the difference in the
mean observed outcomes at $t=0$ among the treated units and the control units,
if the control units had a distribution of ${\bf W}_{i-1}$ similar to the
treated units. The first difference is the identifying functional of the ATT
under the assumption that ${\bf W}_{i0}$ is sufficient to control for
confounding (Assumption 5 below) – see $\tau^{\text{ign.}}$ given in (2). The
second component is equivalent to a similar ignorability estimand using only
pre-treatment data. The form of the expression - a difference in these two
estimands - takes the form of a DiD estimator but one which uses pre-period
outcomes in the conditioning set.

###  2.3 Assumptions and identification of previous approaches

Versions of Assumptions 1-3 have been required by essentially all estimators
in this setting. Differences arise between estimators in the additional
assumption required for identifying $\tau$ from the observed data. We briefly
review the two standard assumptions (justifying ignorability and DiD
estimators, respectively) and discuss how our stable bias assumption relates
to previous approaches.

Ignorability estimators of the ATT require the following assumption

######  Assumption 5 (Conditional ignorability).

$Y_{i1}(0)\perp\\!\\!\\!\perp A_{i}|{\bf W}_{i0}$.

The assumption states that conditional on everything known prior to treatment,
the treatment is as if randomized with respect to the counterfactual outcome
under no treatment at $t=1$. This assumption is alternatively known as
selection on observables or no unmeasured confounding. Under Assumptions 1, 2,
3, and 5, $\tau$ is identified from the observed data by

$\displaystyle\begin{aligned} \tau^{\text{ign.}}&=E\left\\{E(Y_{i1}|{\bf
W}_{i0},A_{i}=1)-E(Y_{i1}|{\bf W}_{i0},A_{i}=0)|A_{i}=1\right\\}\\\
&=E(Y_{i1}|A_{i}=1)-E\left\\{E(Y_{i1}|{\bf
W}_{i0},A_{i}=0)|A_{i}=1\right\\}.\end{aligned}$ (2)

Compare $\tau^{\text{ign.}}$ to $\tau^{\text{gDiD}}$. Note that
$\tau^{\text{gDiD}}$ can be considered as a difference in ignorability
estimators - one applied to the post-treatment period and one applied to the
pre-treatment period. While the conditional ignorability assumption
(Assumption 5) assumes that the bias of $\tau^{\text{ign.}}$ is zero, the
proposed stable bias assumption (Assumption 4) assumes that the bias of
$\tau^{\text{ign.}}$ may be non-zero and may be estimated in the pre-period
data. Note that ignorability implies stable bias as long as there is no pre-
period treatment effect: $E(Y_{i0}|A_{i}=1)-E\left\\{E(Y_{i0}|{\bf
W}_{i-1},A_{i}=0)|A_{i}=1\right\\}=0$, which in turn is equivalent to there
being no unmeasured pre-treatment confounding.

All DiD estimators require some form of the parallel trends assumption. The
typlical conditional parallel trends assumption states that the mean
counterfactual trend in the outcome absent treatment in the treated group is
equal to the observed trend in the control group within strata defined by
time-invariant covariates:

######  Assumption 6 (Conditional parallel trends).

$E\left\\{Y_{i1}(0)-Y_{i0}(0)|{\bf
X}_{i},A_{i}=1\right\\}=E\left\\{Y_{i1}(0)-Y_{i0}(0)|{\bf
X}_{i},A_{i}=0\right\\}.$

Consider a toy example of health status and gender. The traditional
conditional parallel trends assumption (Assumption 6) assumes that the health
status of treated women if they hadn’t been treated would evolve similarly to
the health status of untreated women. Our stable bias assumption (Assumption
4) allows the evolution of health status to depend on initial health status
and accounts for effects of past outcomes on present outcomes. In parallel
trends language, our stable bias assumption assumes that in a counterfactual
setting without treatment, the health status of treated women with low pre-
treatment ($t<0$) health would evolve in parallel to the health status of
untreated women with low pre-treatment health, while accounting for potential
associations between $Y_{i0}$ and $Y_{i1}$.

Further, because the usual identification requires only conditioning on time-
invariant ${\bf X}_{i}$, the version of the positivity assumption required by
classic DiD estimators is slightly weaker:

######  Assumption 7 (Time-invariant covariate positivity).

$P\left\\{\pi_{X}({\bf X}_{i})<1-\epsilon\right\\}=1$ for the propensity score
$\pi_{X}({\bf x})=P(A_{i}=1|{\bf X}_{i}={\bf x})$ and some $\epsilon>0$.

Under Assumptions 1, 2, 6, and 7, $\tau$ is identified from the observed data
by

$\displaystyle\tau^{\text{cDiD}}$ $\displaystyle=E\left\\{E(Y_{i1}|{\bf
X}_{i},A_{i}=1)-E(Y_{i0}\mid{\bf X}_{i},A_{i}=1)-(E(Y_{i1}|{\bf
X}_{i},A_{i}=0)-E(Y_{i0}\mid{\bf X}_{i},A_{i}=0))|A_{i}=1\right\\}$
$\displaystyle=E(Y_{i1}|A_{i}=1)-E\left\\{E(Y_{i1}|{\bf
X}_{i},A_{i}=0)|A_{i}=1\right\\}-E(Y_{i0}|A_{i}=1)+E\left\\{E(Y_{i0}\mid{\bf
X}_{i},A_{i}=0)|A_{i}=1\right\\}.$

Consider a further comparison of these three approaches to identifying causal
effects in this setting. All of these estimators of the ATT are imputation
estimators of the form:

$\displaystyle\tau=E\\{Y_{i1}(1)-Y_{i1}(0)|A_{i}=1\\}=E(Y_{i1}|A_{i}=1)-E\\{\eta({\bf
W}_{i0})|A_{i}=1\\}$

where by assumption the counterfactual $E\\{Y_{i1}(0)|A_{i}=1\\}$ equals
$E\\{\eta({\bf W}_{i0})|A_{i}=1\\}$, which is identified from the observed
data.

Assumption 5 (ignorability) leads to $\eta^{\text{ign.}}({\bf
W}_{i0})=E(Y_{i1}|{\bf W}_{i0},A_{i}=0)$, while under parallel trends
(Assumption 6), the imputation function looks like $\eta^{\text{cDiD}}({\bf
W}_{i0})=Y_{i0}+\Delta({\bf X}_{i})$ for $\Delta({\bf
X}_{i})=E(Y_{i1}-Y_{i0}|{\bf X}_{i},A_{i}=0)$, the conditional trend in the
untreated group. Under stable bias (Assumption 4), the imputation function is
$\eta^{\text{gDiD}}({\bf W}_{i0})=Y_{i0}+\Delta^{*}({\bf W}_{i0})$ for
$\Delta^{*}({\bf W}_{i0})=E(Y_{i1}|{\bf W}_{i0},A_{i}=0)-E(Y_{i0}|{\bf
W}_{i-1},A_{i}=0)=\eta^{\text{ign.}}({\bf W}_{i0})-E(Y_{i0}|{\bf
W}_{i-1},A_{i}=0)$.

Comparing these imputation functions reveals that the complexity of the
ignorability imputation function is essentially unbounded – it can depend
nonlinearly on $Y_{i0}$ and may encode interactions between $Y_{i0},\bar{{\bf
Y}}_{i-1},$ and ${\bf X}_{i}$ – while the DiD imputation function is in
comparison quite simple. It depends strictly linearly on $Y_{i0}$ and can
encode no interactions between $Y_{i0}$ and ${\bf X}_{i}$ and ignores
information in previous pre-period outcomes completely. While the generalized
DiD imputation function includes a term that is linear in $Y_{i0}$, it also
includes $\Delta^{*}({\bf W}_{i0})$, a term that can capture nonlinearity and
interactions in the relationship between $\bar{{\bf Y}}_{i0}$ and $Y_{i1}$ and
is of similar complexity to $\eta^{\text{ign.}}$.

######  Remark 1.

We argue that the gDiD positivity and stable bias assumptions (Assumptions 3
and 4) in comparison to the positivity and parallel trends assumptions of DiD
(Assumptions 7 and 6) represent a kind of bias-variance trade-off. The
differences between these sets of assumptions lie in the conditioning set. If
pre-period outcomes are ignored and ${\bf W}_{it}$ is taken to be ${\bf
X}_{i}$, then Assumption 3 is equivalent to Assumption 7, and Assumption 4 is
equivalent to Assumption 6.

We have further just argued that by using Assumption 4, gDiD is able to use a
much more complex imputation function for the counterfactual
$E\\{Y_{i1}(0)|A_{i}=1\\}$. The price of this additional complexity in the
imputation function – which should in general lead to lower bias – is to make
the positivity assumption stronger. If there is no or only weak overlap in the
pre-period outcomes, Assumption 3 may be less likely to hold than Assumption
7. When Assumption 3 is nearly violated, this can lead to inflated variance in
gDiD (as well as in related ignorability estimators). We demonstrate this
numerically in Section 4. Note also that there is a similar trade-off required
to move from the unconditional parallel trends assumption to Assumption 6: the
bias reduction due to conditioning on ${\bf X}_{i}$ comes at the cost of the
assumption of overlap on ${\bf X}_{i}$. There is, as they say, no free lunch.

###  2.4 Estimation and inference

Having identified $\tau$ in terms of the observed data distribution,
estimation of $\tau$ amounts to efficient estimation of the expectations
$E\\{\mu_{1}({\bf W}_{i0})|A_{i}=1\\}=E\left\\{E(Y_{i1}|{\bf
W}_{i0},A_{i}=0)|A_{i}=1\right\\}$ and $E\\{\mu_{0}({\bf
W}_{i-1})|A_{i}=1\\}=E\left\\{E(Y_{i0}|{\bf
W}_{i-1},A_{i}=0)|A_{i}=1\right\\}$. We propose using an influence-function
based approach. The influence functions for these functionals [9] can be
expressed as

$\displaystyle E(A_{i})^{-1}\left[\frac{(1-A_{i})\pi({\bf
W}_{i0})Y_{i1}+\\{A_{i}-\pi({\bf W}_{i0})\\}\mu_{1}({\bf W}_{i0})}{1-\pi({\bf
W}_{i0})}-E\\{\mu_{1}({\bf W}_{i0})|A_{i}=1\\}\right]\text{ and }$
$\displaystyle E(A_{i})^{-1}\left[\frac{(1-A_{i})\pi_{0}({\bf
W}_{i-1})Y_{i0}+\\{A_{i}-\pi_{0}({\bf W}_{i-1})\\}\mu_{0}({\bf
W}_{i-1})}{1-\pi_{0}({\bf W}_{i-1})}-E\\{\mu_{0}({\bf
W}_{i-1})|A_{i}=1\\}\right]$

where $\pi_{0}({\bf W}_{i-1})=P(A_{i}=1|{\bf X}_{i}={\bf x},\bar{{\bf
Y}}_{i-1}=\bar{{\bf y}}_{-1}).$ Because the efficient influence function for a
sum is the sum of the influence functions, we have therefore that the
influence function for (1) can be written as

$\displaystyle\begin{aligned} \phi({\bf
O}_{i})&=E(A_{i})^{-1}\left[Y_{i1}A_{i}-\frac{(1-A_{i})\pi({\bf
W}_{i0})Y_{i1}+\\{A_{i}-\pi({\bf W}_{i0})\\}\mu_{1}({\bf W}_{i0})}{1-\pi({\bf
W}_{i0})}-\right.\\\
&\qquad\qquad\left.Y_{i0}A_{i}+\frac{(1-A_{i})\pi_{0}({\bf
W}_{i-1})Y_{i0}+\\{A_{i}-\pi_{0}({\bf W}_{i-1})\\}\mu_{0}({\bf
W}_{i-1})}{1-\pi_{0}({\bf W}_{i-1})}\right]-\tau^{\text{gDiD}}.\end{aligned}$
(3)

Thus, a doubly-robust estimator for $\tau$ may be constructed as

$\displaystyle\begin{aligned}
{\widehat{\tau}}&=n_{1}^{-1}\sum_{i=1}^{n}\left[Y_{i1}A_{i}-\frac{(1-A_{i}){\widehat{\pi}}({\bf
W}_{i0})Y_{i1}+\\{A_{i}-{\widehat{\pi}}({\bf
W}_{i0})\\}{\widehat{\mu}}_{1}({\bf W}_{i0})}{1-{\widehat{\pi}}({\bf
W}_{i0})}-\right.\\\
&\qquad\qquad\left.Y_{i0}A_{i}+\frac{(1-A_{i}){\widehat{\pi}}_{0}({\bf
W}_{i-1})Y_{i0}+\\{A_{i}-{\widehat{\pi}}_{0}({\bf
W}_{i-1})\\}{\widehat{\mu}}_{0}({\bf W}_{i-1})}{1-{\widehat{\pi}}_{0}({\bf
W}_{i-1})}\right]\end{aligned}$ (4)

where
${\widehat{\pi}},{\widehat{\pi}}_{0},{\widehat{\mu}}_{1},{\widehat{\mu}}_{0}$
are estimates of the nuisance functions $\pi,\pi_{0},\mu_{1},\mu_{0}$ and
$n_{1}=\sum_{i=1}^{n}A_{i}$. There are two predominant strategies for
estimating the nuisance functions. When the dataset is small, parametric
models may be specified for the nuisance functions, as in, for example, [8]
and [25]. However, parametric models are difficult to specify and often incur
bias due to model misspecification, so when the dataset is not small we
propose to use non-parametric machine learning methods to estimate the
nuisance functions ($\pi,\pi_{0},\mu_{1},\mu_{0}$) in combination with cross-
fitting [10]. Cross-fitting separates the estimation of the nuisance functions
from their evaluation in (4), allowing only weak assumptions on the complexity
of the underlying nuisance functions to achieve root-$n$ consistent and
asymptotically normal estimators of the causal effect [18]. Intuitively,
cross-fitting prevents a kind of overfitting that might otherwise occur when
using non-parametric estimation methods.

While cross-fitting has been described extensively elsewhere (see, e.g.,
[18]), we briefly review its implementation. Specifically, we first randomly
splitting the data into $K$ disjoint subsets (indexed from $k=1,...,K$) with
$N_{k}$ observations, respectively. Next, for each subset, one estimates the
nuisance functions on all data except the $k$-th subset and constructs the
estimator $\hat{\tau}_{k}$ on the the $k$th subset. Finally, the cross-fit
estimator is constructed as the sample-size weighted average of each of the
$K$ estimators, i.e. $\sum_{k=1}^{K}\frac{N_{k}}{N}\hat{\tau}_{k}$.

Under minimal assumptions on the nuisance function estimates,
${\widehat{\tau}}$ will have an asymptotic normal distribution. This means
that we may do standard inference and construct confidence intervals for
$\tau$, even when using machine learning methods for estimation. Specifically,
we employ the following assumption on the rates of convergence of the nuisance
function estimates:

######  Assumption 8 (Nuisance function convergence).

$\displaystyle E\left[\left\\{{\widehat{\mu}}_{1}({\bf W}_{i0})-\mu_{1}({\bf
W}_{i0}\right\\}^{2}\right]^{1/2}\times E\left[\left\\{{\widehat{\pi}}({\bf
W}_{i0})-\pi({\bf W}_{i0}\right\\}^{2}\right]^{1/2}=o_{p}(n^{-\frac{1}{2}})$
(5) $\displaystyle E\left[\left\\{{\widehat{\mu}}_{0}({\bf
W}_{i-1})-\mu_{0}({\bf W}_{i-1}\right\\}^{2}\right]^{1/2}\times
E\left[\left\\{{\widehat{\pi}}_{0}({\bf W}_{i-1})-\pi_{0}({\bf
W}_{i-1}\right\\}^{2}\right]^{1/2}=o_{p}(n^{-\frac{1}{2}})$ (6)

The cross-fitting approach we have proposed in Section 2.4 along with most
off-the-shelf machine learning estimators is generally sufficient to ensure
Assumption 5, as has been noted in other contexts previously [10], such as for
regression trees/random forests[30] or deep learning[14]. This allows us to
state the limiting distribution of the proposed estimator

######  Proposition 1.

Under Assumptions 1, 2, 3, 4, 5, and additional regularity conditions in
Appendix A,

$\displaystyle n^{-1/2}\sigma^{-1}({\widehat{\tau}}-\tau)\longrightarrow
N(0,1),\qquad\sigma^{2}=E\\{\phi({\bf O}_{i})^{2}\\}.$ (7)

The proof of this proposition is ensured by previous results for the ATT [10,
18], as our proposed estimator is the difference between two ATTs. We give
details in Appendix A. And, as with previous results, Proposition 1 continues
to hold if $\sigma$ is replaced by ${\widehat{\sigma}}$ defined by

$\displaystyle{\widehat{\sigma}}^{2}$
$\displaystyle=n_{1}^{-1}\sum_{i=1}^{n}\left(\left[Y_{i1}A_{i}-\frac{(1-A_{i}){\widehat{\pi}}({\bf
W}_{i0})Y_{i}+\\{A_{i}-{\widehat{\pi}}({\bf
W}_{i0})\\}{\widehat{\mu}}_{1}({\bf W}_{i0})}{1-{\widehat{\pi}}({\bf
W}_{i0})}-\right.\right.$ (8)
$\displaystyle\qquad\qquad\left.\left.Y_{i0}A_{i}+\frac{(1-A_{i}){\widehat{\pi}}_{0}({\bf
W}_{i-1})Y_{i}+\\{A_{i}-{\widehat{\pi}}_{0}({\bf
W}_{i-1})\\}{\widehat{\mu}}_{0}({\bf W}_{i-1})}{1-{\widehat{\pi}}_{0}({\bf
W}_{i-1})}\right]-{\widehat{\tau}}\right)^{2}$ (9)

The immediate consequence of Proposition 1 is that a (1-$\alpha$)% confidence
interval can therefore be constructed as ${\widehat{\tau}}\pm
z_{1-\alpha/2}{\widehat{\sigma}}$ where $z_{\gamma}$ is the $\gamma$th
quantile of the standard normal distribution. We review an alternative
multiplier bootstrap variance estimator in Appendix B.

######  Remark 2.

Assumption 5 will not be ensured when parametric models are used to estimate
nuisance functions without sample splitting unless both models are correctly
specified. In Appendix C we give a general form for the asymptotic variance
under no assumptions about correct model specification. Furthermore, we also
note that it would be straightforward to compute design-based variances in the
style of [21], which may be of interest in settings where units cannot be
thought of as arising from a super-population of units.

######  Remark 3.

The use of non-parametric estimation methods in combination with influence-
function based estimators has frequently been discouraged in the DiD
literature. For example, [25] propose using generalized linear models in
combination with their influence-function based DiD estimator, arguing that
the “curse of dimensionality” makes non-parametric methods less useful in
practice. While we agree with this caution in a small-sample setting, non-
parametric estimation methods can yield asymptotically normal and root-n
consistent estimates under relatively weak conditions, even in a high-
dimensional setting. For example, as noted by [18], if we assume that all true
model propensity-score and regression models lie within some $s$-sparse model
class, lasso-based estimates of the nuisance parameters can achieve root-n
consistent and asymptotically normal estimates provided $s=o(\sqrt{n}/\log d)$
and that we’re able to achieve $O_{p}(\sqrt{\frac{s\log d}{n}})$ nuisance
estimation rates in the $L_{2}$ norm. By contrast, when using parametric
models in combination with influence-function based parametric estimators, all
models are likely misspecified in practice leading to inconsistent estimators.
Influence-based estimators using parametric models for the nuisance parameters
may be relatively less useful than non-parametric models. We therefore
recommend using modern non-parametric and machine methods in combination with
these estimators when the sample size is not small. Our accompanying R package
implements our proposed method using SuperLearner, an ensemble of machine
learning methods; however, our package can also accommodate generalized linear
models for estimation if desired.

##  3 Extensions

###  3.1 Staggered adoption and group-time treatment effects

In settings where units may adopt or be exposed to the treatment at different
times (the so-called staggered adoption setting), the basic framework of
generalized DiD may be extended to estimate ATTs that capture potentially
heterogeneous treatment effects over time. We first establish some additional
notation.

Let ${\bf A}_{i}=(A_{it})_{t=1,...,T}$ be a vector of treatment indicators
over time, where now time is indexed from $t=1$ to $t=T$. Let ${\cal A}_{i}$
be the time of adoption ${\cal A}_{i}=\min_{A_{it}=1}t$ if any $A_{it}=1$ and
is $\infty$ otherwise, and let ${\cal P}_{i}={\cal A}_{i}-1$ be the last pre-
period time point. We also redefine the counterfactuals so that they
potentially depend on the entire treatment history prior to the current time,
$Y_{it}(\bar{{\mbox{\boldmath$a$}}}_{t})$ for
$\bar{{\mbox{\boldmath$a$}}}_{t}=(a_{s})_{s=1,..,t}$, and let ${\cal
A}_{\bar{{\mbox{\boldmath$a$}}}_{t}}$ be the implied time of adoption for a
given treatment history $\bar{{\mbox{\boldmath$a$}}}_{t}$: ${\cal
A}_{\bar{{\mbox{\boldmath$a$}}}_{t}}=\min_{a_{s}=1}s$. Define $\bar{{\bf
A}}_{it}=(A_{is})_{s=1,...,t}$, and let ${\cal
P}_{\bar{{\mbox{\boldmath$a$}}}_{t}}={\cal
A}_{\bar{{\mbox{\boldmath$a$}}}_{t}}-1$.

Then we may define the group-time treatment effect

$\displaystyle\tau_{\bar{{\mbox{\boldmath$a$}}}_{t}}=E\\{Y_{it}(\bar{{\mbox{\boldmath$a$}}}_{t})-Y_{it}({\bf
0}_{t})|\bar{{\bf A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}\\}.$

Notice that $\tau_{\bar{{\mbox{\boldmath$a$}}}_{t}}$ actually generalizes the
estimands considered in [8] as it places no restrictions on the treatment
history. That is, we do not require that once a unit is treated it remains
treated.

The assumptions required to identify $\tau_{\bar{{\mbox{\boldmath$a$}}}_{t}}$
are as follows. First, we update the consistency assumption for the new
counterfactuals,

######  Assumption 9 (Consistency, staggered adoption setting).

$Y_{it}(\bar{{\mbox{\boldmath$a$}}}_{t})=Y_{it}\text{ when }\bar{{\bf
A}}_{t}=\bar{{\mbox{\boldmath$a$}}}_{t}.$

######  Assumption 10 (No anticipation, staggered adoption setting).

$Y_{is}(\bar{{\mbox{\boldmath$a$}}}_{t})=Y_{is},\text{ when }s<{\cal
A}(\bar{{\mbox{\boldmath$a$}}}_{t})$.

######  Assumption 11 (Positivity, staggered adoption setting).

$P\left\\{\pi_{\bar{{\mbox{\boldmath$a$}}}_{t}}({\bf W}_{i{\cal
P}_{\bar{{\mbox{\boldmath$a$}}}_{t}}})<1-\epsilon\right\\}=1$, for
$\pi_{\bar{{\mbox{\boldmath$a$}}}_{t}}({\bf w})=P(\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}|{\bf W}_{i{\cal
P}_{\bar{{\mbox{\boldmath$a$}}}_{t}}}={\bf w})$.

######  Assumption 12 (Stable bias, staggered adoption setting).

$\displaystyle\begin{aligned} &E\left[E\left\\{Y_{it}({\bf 0}_{t})|{\bf
W}_{i{\cal P}_{\bar{{\mbox{\boldmath$a$}}}_{t}}},\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}\right\\}-E\left\\{Y_{it}|{\bf
W}_{i{\cal P}_{\bar{{\mbox{\boldmath$a$}}}_{t}}},\bar{{\bf A}}_{it}={\bf
0}_{t}\right\\}|\bar{{\bf A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}\right]=\\\
&\qquad E\left[E\left\\{Y_{i0}|{\bf W}_{i{\cal
P}_{\bar{{\mbox{\boldmath$a$}}}_{t}}-1},\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}\right\\}-E\left\\{Y_{i0}|{\bf
W}_{i{\cal P}_{\bar{{\mbox{\boldmath$a$}}}_{t}}-1},\bar{{\bf A}}_{it}={\bf
0}_{t}\right\\}|\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}\right].\end{aligned}$

These assumptions directly generalize Assumptions 1, 2, 3, and 4 and serve the
same purposes. Under Assumptions 9, 10, 11, and 12, the group-time ATT
$\tau_{\bar{{\mbox{\boldmath$a$}}}_{t}}$ may be estimated from the observed
data via the following functional

$\displaystyle\begin{aligned}
E\\{Y_{it}(\bar{{\mbox{\boldmath$a$}}}_{t})-Y_{it}({\bf 0}_{t})|\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}\\}&=E(Y_{it}|\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t})-E\left\\{E(Y_{it}|{\bf
W}_{i0},A_{i}=0)|\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}\right\\}-\\\ &\qquad
E(Y_{i0}|\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t})+E\left\\{E(Y_{i0}|{\bf
W}_{i-1},A_{i}=0)|\bar{{\bf
A}}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t}\right\\}\end{aligned}$ (10)

The influence function and estimation for (10) follow directly from Section
2.4, and the asymptotics follow precisely the same arguments as in that
section.

We can also aggregate these functionals over time, following [8]. Let
$\Xi_{t}\subset\\{0,1\\}^{t}$ be a set of treatment histories such that

$\displaystyle\Xi_{t}=\left\\{\bar{a}_{t}:P(\bar{A}_{it}=\bar{{\mbox{\boldmath$a$}}}_{t})>0\text{
and Assumptions \ref{asmp:consistency-sa}, \ref{asmp:anticipation-sa},
\ref{asmp:positivity-sa}, and \ref{stable-bias-sa} hold.}\right\\}$

Further let
$\boldsymbol{\omega}_{t}=\\{\omega_{\bar{{\mbox{\boldmath$a$}}}_{t}}\\}_{\bar{{\mbox{\boldmath$a$}}}_{t}\in\Xi_{t}}$
be a set of weights defined over the elements $\Xi_{t}$. Then we can define
aggregated treatment effects as

$\displaystyle\tau_{\boldsymbol{\omega}_{t}}=\sum_{\bar{{\mbox{\boldmath$a$}}}_{t}\in\Xi_{t}}\omega_{\bar{{\mbox{\boldmath$a$}}}_{t}}\tau_{\bar{{\mbox{\boldmath$a$}}}_{t}}$
(11)

Examples of weights leading to interpretable treatment effects include, the
average effect of the intervention in year $t$ among units that are treated at
year $t$, given by
$\omega_{\bar{{\mbox{\boldmath$a$}}}_{t}}=I\\{a_{t}=1\\}/\sum_{\bar{{\mbox{\boldmath$a$}}}_{t}\in\Xi_{t}}I\\{a_{t}=1\\}$.
Another example is the effect of the intervention in time $t$ among units that
adopted treatment at time $s$, given by
$\omega_{\bar{{\mbox{\boldmath$a$}}}_{t}}=I\\{{\cal
A}(\bar{{\mbox{\boldmath$a$}}}_{t})=s\\}/\sum_{\bar{{\mbox{\boldmath$a$}}}_{t}\in\Xi_{t}}I\\{{\cal
A}(\bar{{\mbox{\boldmath$a$}}}_{t})=s\\}$. A final example is the effect of
the intervention at time $t$ among units that adopted treatment at time $s$
and were also treated at time $t$, given by
$\omega_{\bar{{\mbox{\boldmath$a$}}}_{t}}=I\\{{\cal
A}(\bar{{\mbox{\boldmath$a$}}}_{t})=s,a_{t}=1\\}/\sum_{\bar{{\mbox{\boldmath$a$}}}_{t}\in\Xi_{t}}I\\{{\cal
A}(\bar{{\mbox{\boldmath$a$}}}_{t})=s,a_{t}=1\\}$.

######  Remark 4.

Just as in Section 2.3, if the pre-period outcomes are omitted from the
covariate vector, then Assumption 12 reduces to an assumption similar to the
parallel trends assumption of [8]. However, we are using a generalized version
of their estimand allowing identification of treatment effects when units may
switch between treatment and control group unrestrictedly.

###  3.2 Clustered treatment assignment

In many settings we observe clusters $j=1,...,C$; and units $i=1,...,n_{j}$
within each group, where the $j$th of $C$ clusters has $n_{j}$ observations
indexed by $i$. Treatment $A_{j}$ is then assigned to all units within each
cluster. We assume that the data in the $j$th cluster ${\bf
O}_{j}=\left\\{Y_{ji1}(0),Y_{ji1}(1),{\bf X}_{ji},A_{j},\bar{{\bf
Y}}_{ji0}\right\\}$ is independent of the data in the $k$th cluster ($k\neq
j$), ${\bf O}_{k}=\left\\{Y_{ki1}(0),Y_{ki1}(1),{\bf X}_{ki},A_{k},\bar{{\bf
Y}}_{ki0}\right\\}$, and, for simplicity, we assume that unit-cluster
assignment is fixed (i.e. units cannot change clusters). We define the ATT at
time-period $t=1$ to be:

$\displaystyle\tau=E\\{Y_{ji1}(1)-Y_{ji1}(0)|A_{j}=1\\}.$ (12)

The outcome in this setting may depend on characteristics of the individual
observation captured in $\bar{{\bf W}}_{ji0}=({\bf X}_{ji},\bar{{\bf
Y}}_{ji0}),$ and it may also depend on characteristics of the other members of
the cluster $\\{\bar{{\bf W}}_{jk0}:k\neq i\\}$. Define a set of cluster-level
summary variables ${\widetilde{\bf W}}_{j0},{\widetilde{\bf W}}_{j-1}$ – e.g.,
${\widetilde{\bf W}}_{jt}=n_{j}^{-1}\sum_{i=1}^{n_{j}}{\bf W}_{jit}$ – and
collect both the individual- and cluster-level information into the vector
${\bf V}_{jit}=({\bf W}_{jit},{\widetilde{\bf W}}_{jt})$. Let
$N_{1}=\sum_{j=1}^{C}n_{j}A_{j}$ be the total number of treated observations.
Then we can write an influence-function based estimator for (12) as

$\displaystyle\begin{aligned}
{\widehat{\tau}}&=N_{1}^{-1}\sum_{j=1}^{C}\sum_{i=1}^{n_{j}}\left[Y_{ji1}A_{j}-\frac{(1-A_{j}){\widehat{p}}({\widetilde{\bf
W}}_{j0})Y_{ji1}+\\{A_{j}-{\widehat{p}}({\widetilde{\bf
W}}_{j0})\\}{\widehat{m}}_{1}({\bf V}_{ji0})}{1-{\widehat{p}}({\widetilde{\bf
W}}_{j0})}-\right.\\\
&\qquad\qquad\left.Y_{ji0}A_{j}+\frac{(1-A_{j}){\widehat{p}}_{0}({\widetilde{\bf
W}}_{j-1})Y_{ji0}+\\{A_{j}-{\widehat{p}}_{0}({\widetilde{\bf
W}}_{j-1})\\}{\widehat{m}}_{0}({\bf
V}_{ji-1})}{1-{\widehat{p}}_{0}({\widetilde{\bf
W}}_{j-1})}\right]&=C^{-1}\sum_{j=1}^{C}{\widetilde{\phi}}({\bf
O}_{j})\end{aligned}$ (13)

where ${\widehat{p}}({\bf v})$ is an estimate of $P(A_{j}=1|{\bf V}_{ji0}={\bf
v})$, ${\widehat{p}}_{0}({\bf v})$ is an estimate of $P(A_{j}=1|{\bf
V}_{ji-1}={\bf v})$, and ${\widetilde{\phi}}({\bf O}_{j})$ is an uncentered
version of the influence function

$\displaystyle\phi({\bf O}_{j})$
$\displaystyle=\frac{C}{N_{1}}\sum_{i=1}^{n_{j}}\left[Y_{ji1}A_{j}-\frac{(1-A_{j})p({\bf
V}_{ji0})Y_{ji1}+\\{A_{j}-p({\bf V}_{ji0})\\}m_{1}({\bf V}_{ji0})}{1-p({\bf
V}_{ji0})}-\right.$
$\displaystyle\qquad\qquad\left.Y_{ji0}A_{j}+\frac{(1-A_{j})p_{0}({\bf
V}_{ji-1})Y_{ji0}+\\{A_{j}-p_{0}({\bf V}_{ji-1})\\}m_{0}({\bf
V}_{ji-1})}{1-p_{0}({\bf V}_{ji-1})}-\tau\right],$

######  Proposition 2.

Under Assumptions 18, 19, 20, 21, 22, and 24 in Appendix D,

$\displaystyle C^{-1/2}({\widehat{\tau}}-\tau)\longrightarrow N(0,\zeta^{2})$

for $\zeta^{2}=E\\{\phi({\bf O}_{j})^{2}\\}$.

Proof given in Appendix D. This proposition states that the estimator under
clustered treatment assignment has a similar asymptotic distribution to (7),
the distribution for the estimator in the iid case. A plug-in version of
$\zeta^{2}$ may be computed just as in (8), but we also show in Appendix B
that a multiplier bootstrap estimator may equivalently be used.

###  3.3 Other treatment effects

Traditional DiD and the parallel trends assumption have been engineered
specifically to estimate the ATT, but our alternative formulation of the
stable bias assumption opens the door to estimating many other, equally or
more relevant effects. For example, the Centers of Medicare & Medicaid
Services may implement a voluntary intervention that is adopted by only a
subset of eligible insurance plans participating in the Medicare Advantage
program. The ATT estimates the expected effect of the program if other similar
plans also adopted the intervention. However, policymakers may also be
interested in the effect if all plans were to implement the intervention,
which would be captured by the so-called average treatment effect (ATE):

$\displaystyle\tau_{\text{ATE}}=E\\{Y_{i1}(1)-Y_{i1}(0)\\}.$

Similarly, in the example in Section 5, we examine the effect of health system
affilliation on quality of care delivered by physician organizations (POs) to
Medicare beneficiaries. POs that affiliate with health systems tend to be
larger and tend to have fewer racial and ethnic minority beneficiaries.
Policymakers may not be interested particularly in this subset of POs but may
prefer to estimate effects of health system affiliation in all POs.

The assumptions needed to identify the ATE in the context of unmeasured
confounding and the tools needed to estimate it are nearly identical as to the
ATT, which we explored in detail above. The positivity assumption (Assumption
3 must be updated to be a little stricter, requiring positive probability of
being either treated or untreated for all covariate strata, and the stable
bias assumption must also be updated:

######  Assumption 13 (Positivity, ATE).

$P\\{\epsilon<\pi({\bf W}_{i0})<1-\epsilon\\}=1$

######  Assumption 14 (Stable bias, ATE).

$\displaystyle E\left\\{E(Y_{i1}|A=1,{\bf W}_{i0})-E(Y_{i1}|A=0,{\bf
W}_{i0})\right\\}-\tau_{\text{ATE}}$ $\displaystyle=E\left\\{E(Y_{i0}|A=1,{\bf
W}_{i-1})-E(Y_{i0}|A=0,{\bf W}_{i-1})\right\\}.$

Under these two assumptions, as well as Assumptions 1 and 2, we may estimate
the ATE in a manner analogous to (4).

######  Remark 5.

While we argue that alternative estimands – like the ATE – may have more
policy relevance in some cases, stable bias assumptions for these alternatives
may be harder to justify. Intuitively, pre-treatment data arguably provides
information about $Y_{i1}(0)$; however, it is harder to argue that these same
data provide information about $Y_{i1}(1)$. For causal estimands such as the
ATE that effectively require imputing $Y_{i1}(1)$ for untreated units, it
requires further study to establish that pre-treatment biases with respect to
$Y_{it}(0)$ (for $t<1$) can inform post-treatment biases with respect to
averages of these unobserved counterfactuals. Nevertheless, these sorts of
stable bias assumptions may be reasonable in some settings, and in these cases
we can follow the logic above to identify and estimate these causal
quantities.

##  4 Simulations

We conduct two complementary sets of simulation studies. The first closely
follows the setup from [25] and compares the performance of the proposed gDiD
estimator to state-of-the-art DiD estimation methods. However, the space of
possible data-generating mechanisms in the panel data setting is much broader
than considered in [25]. We therefore conduct a second set of simulation
studies that allows a more thorough exploration of the performance of
different estimators under a range of assumptions about the data-generating
processes. For all simulations, we fit all models on five-hundred random
datasets.

We consider the following estimators:

  * •

Three implementations of gDiD: one where we control for only one pre-period
outcome (gDiD-1) – ${\bf W}_{it}=({\bf X}_{i},Y_{it})$ – one where we control
for two pre-period outcomes (gDiD-2) – ${\bf W}_{it}=({\bf
X}_{i},Y_{it},Y_{it-1})$ – and one where we only use time-invariant covariates
(gDiD-0) – ${\bf W}_{it}={\bf X}_{i}$. The gDiD-0 estimator is similar to the
estimator in [8] and [25] but using machine learning models instead of
parametric models.

  * •

Two ignorability estimators, controlling for one (Ign-1) and two (Ign-2) pre-
period outcomes, using augmented inverse probability weighting.

  * •

Augmented synthetic controls [5, 4] controlling for three pre-period outcomes.
This is essentially an additional (regression-based) ignorability estimator
using regularized linear models [7].

  * •

A standard DiD estimator.

  * •

The doubly-robust conditional DiD estimator given by [25].

For the Ign and gDiD estimators, we use cross-fitting and a weighted
combination of GLMs, Random Forests, and GAMs, where the weights on each model
are learned using SuperLearner [27, 19]. For all estimators except for
augmented synthetic controls, we conduct inference using the empirical
variance of the estimated influence functions and standard normal quantiles.
We do not conduct inference for augmented synthetic controls because there are
not enough pre-treatment time-periods to use the default conformal inference
based approach [12] implemented in the augsynth package. Below, we report
bias, root mean squared error (RMSE), average 95% confidence interval length
(CIL), and 95% confidence interval coverage (Coverage).

Before getting into details about the simulations, we briefly give motivation
and high-level results. In the first simulation, we modify the simulations in
[25] such that a variable required for conditional parallel trends is not
available. Here we show that, while DiD estimators incur bias due to the
omitted variable, gDiD and ignorability estimators are able to perform well
because pre-period outcomes are able to proxy for the omitted variable. We
further show how positivity violations can degrade the performance of
ignorability estimators and how, in this setting, gDiD suffers less from
positivity issues than ignorability estimators.

In the second simulation setting, we start from a setting where parallel
trends holds but ignorability does not and show that gDiD is competitive with
DiD (ignorability estimators perform poorly). We then go on to vary the
simulation parameters away from parallel trends in various ways and show that
(generally speaking) gDiD outperforms all other estimators (often by a factor
of two or more in terms of RMSE) no matter how we vary the simulation.

###  4.1 First data-generating process

####  4.1.1 Setup

We first conduct a simulation study that closely follows [25]. However, we
make two key modifications to their specification: first, we change the data-
generating process to change the levels of overlap in pre-treatment outcomes;
second, we consider a setting where conditional parallel trends does not hold
given the observed covariates. To be precise, we consider a generic four-
dimensional covariate ${\bf W}=(W_{1},W_{2},W_{3},W_{4})$ and define
regression and propensity score functions: $f_{reg}({\bf
W})=205+27.4W_{1}+13.7(W_{2}+W_{3}+W_{4})$, $f_{ps}({\bf
W})=0.75(-W_{1}+0.5W_{2}-0.5W_{3}-0.25W_{4})$. We take $n=1,000$ independent
draws of a four-dimensional covariate vector ${\bf X}$ from $MVN(0,I_{4})$,
where $I_{4}$ is a 4x4 identity matrix. We then take nonlinear transformations
of ${\bf X}$, which we call $\tilde{{\bf Z}}$, and then standardize
$\tilde{{\bf Z}}$ to obtain ${\bf Z}$. Specifically, define
$\tilde{Z}_{1}=\exp(0.5X_{1})$, $\tilde{Z}_{2}=10+X_{2}/(1+\exp(X_{1}))$,
$\tilde{Z}_{3}=(0.6+X_{1}X_{3}/25)^{3}$, and
$\tilde{Z}_{4}=(20+X_{2}+X_{4})^{2}$. For each covariate $j=1,...,4$, we
define
$Z_{j}=(\tilde{Z}_{j}-\mathbb{E}[\tilde{Z}_{j}])/\sqrt{Var(\tilde{Z}_{j})}$.
For $t=1,...,4$, we build on [17] and [25] and consider the following data-
generating process: $Y_{t}(1)=Y_{t}(0)=0.1tf_{reg}({\bf X})+v(A,{\bf
X};\zeta)+\epsilon_{t}$, $A=1(\text{expit}(f_{ps}({\bf X}))\geq U)$, where
$\epsilon_{t}$ is an independent standard normal random variable, $U$ is an
independent random uniform variable, and $v(A,{\bf X};\zeta)$ is an
independent normal random variable with mean $f_{reg}({\bf X})(1+\zeta A)$ and
variance one. We change the parameter $\zeta$ to explore different levels of
overlap in the pre-treatment outcomes. When there is less overlap in the pre-
treatment outcomes, Assumption 3 becomes (nearly) violated.

We consider two sets of observed data ${\bf O}=({\bf O}_{i})_{i=1,...n}$,
which we call the linear and non-linear cases. In the linear case, we consider
data where the true covariates $X_{2},X_{3},X_{4}$ are available for analysis
${\bf O}_{i}=(Y_{i-2},...,Y_{i1i},A_{i},X_{i2},...,X_{i4})$. In the non-linear
case, only the transformed covariates $Z_{2},Z_{3},Z_{4}$ are available but
not the $X$s: ${\bf O}_{i}=(Y_{i-2},...,Y_{i1},A_{i},X_{i1},...,X_{i3})$. In
both cases the observed data do not include the first covariate $X_{i1}$ or
$Z_{i1}$. In contrast to the simulation in [25], this simulation includes
time-varying unobserved heterogeneity due to omission of $X_{1}$ and $Z_{1}$.
As a result, conditional parallel trends does not hold in our simulations,
though it would if we observed all four covariates.333Notice that $X_{1}$ has
the strongest confounding effect on the outcome as it has the strongest effect
on both the treatment assignment and outcome models. We allow
$\zeta\in\\{0,0.1,0.3,0.5\\}$, where the higher we set $\zeta$ the less
overlap in the pre-treatment outcome levels between the treatment and control
groups. When $\zeta=0$, both ignorability and conditional parallel trends
would hold were we to observe all covariates. We choose this specification
because we believe that in practice parallel trends is unlikely to hold
exactly, as we are not likely to measure all relevant confounders.

####  4.1.2 Results

Table 1 presents the results in both the linear and nonlinear settings when
$\zeta=0$. The results illustrate that the estimators that control for pre-
period outcomes (gDiD-1, gDiD-2, Ign-1, Ign-2, Augsynth) have bias close to
zero in both settings. This suggests that controlling for pre-period outcomes
is sufficient to control for confounding, while the estimators that rely
strictly on parallel trends suffer from bias. While the gDiD-1 and gDiD-2
estimators have low bias, their RMSEs are somewhat higher than the
ignorability-based estimation approaches (Ign-1, Ign-2, and AugSynth). The
additional bias correction from gDiD appears unnecessary in this simulation,
and the estimation of the bias therefore induces additional variability and no
bias. The worst performing estimators in terms of bias are the DiD, cDiD, and
gDiD-0 estimators. Interestingly, all three of these estimators perform
comparably, possibly a function of the relatively small confounding effects of
the included covariates. Overall, these results illustrate how gDiD and
ignorability estimators that take advantage of pre-period outcomes can
overcome omitted variables in some settings and how gDiD is competitive with
ignorability estimators when ignorability holds.

Table 1: Simulation study one results: $\zeta=0$ | Linear | Nonlinear
---|---|---
Estimator | Bias | RMSE | CIL | Coverage | Bias | RMSE | CIL | Coverage
DID | -0.21 | 0.312 | 0.966 | 88.2 | -0.224 | 0.344 | 0.965 | 84
cDID | -0.194 | 0.284 | 0.788 | 83.7 | -0.237 | 0.349 | 0.899 | 80.3
gDID-0 | -0.199 | 0.289 | 0.794 | 82.8 | -0.235 | 0.333 | 0.87 | 80.2
Ign-1 | -0.003 | 0.101 | 0.38 | 93.8 | -0.006 | 0.099 | 0.388 | 96
Ign-2 | 0.002 | 0.088 | 0.334 | 93.2 | -0.006 | 0.085 | 0.335 | 94.2
AugSynth | 0.006 | 0.185 | - | - | -0.002 | 0.149 | - | -
gDID-1 | -0.004 | 0.172 | 0.657 | 94 | 0.005 | 0.175 | 0.666 | 95.2
gDID-2 | 0.003 | 0.134 | 0.508 | 94.5 | -0.002 | 0.133 | 0.51 | 94.8

Table 2 presents results when setting $\zeta=0.1$. Note that in this setting
none of the key identifying assumptions (ignorability, conditional parallel
trends, or stable bias) hold. Changing $\zeta$ also causes a decreased level
of overlap between the pre-treatment outcomes in the treated and control
groups and causes a near violation of Assumption 3. These changes do not
effect the DiD estimators but have two implications for the performance of the
gDiD and ignorability estimators. First, because stable bias and ignorability
do not hold, these estimators incur much more bias than when $\zeta=0$.
Second, because Assumption 3 is closer to being violated, the variances and
associated confidence intervals get larger. These changes lead to unacceptably
large biases in the ignorability estimators, but the gDiD estimators remain
competitive with the DiD and cDiD estimators, with gDiD-1 performing the best
in both the linear and non-linear settings. In Appendix E.1 we display results
setting $\zeta=0.3,0.5$. The results are consistent with this trend, with
worsening performance among the gDiD and Ign estimators, and the performance
of the cDiD and DiD estimators remaining consistent.

Table 2: Simulation study one results: $\zeta=0.1$ | Linear | Nonlinear
---|---|---
Estimator | Bias | RMSE | CIL | Coverage | Bias | RMSE | CIL | Coverage
DID | -0.219 | 0.333 | 0.966 | 85.8 | -0.213 | 0.322 | 0.964 | 85.8
cDID | -0.194 | 0.278 | 0.788 | 86 | -0.229 | 0.322 | 0.887 | 83.5
gDID-0 | -0.197 | 0.285 | 0.794 | 85 | -0.225 | 0.317 | 0.87 | 82.8
Ign-1 | -1.549 | 1.553 | 0.489 | 0 | -1.553 | 1.556 | 0.487 | 0
Ign-2 | -2.381 | 2.386 | 0.656 | 0 | -2.392 | 2.396 | 0.652 | 0
AugSynth | -3.218 | 3.252 | - | - | -3.252 | 3.28 | - | -
gDID-1 | 0.128 | 0.227 | 0.707 | 89.8 | 0.13 | 0.226 | 0.698 | 87.8
gDID-2 | 0.221 | 0.317 | 0.788 | 78.7 | 0.215 | 0.324 | 0.771 | 78.8

###  4.2 Second data-generating process

####  4.2.1 Setup

The previous simulations make somewhat restrictive assumptions on the data-
generating processes. For example, it is unlikely that coefficients on the
unobserved covariates change linearly with time. Additionally, prior period
outcomes may directly enter the outcome model. However, the above
specification becomes complex when incorporating these nuances due to the
presence of covariates. We instead consider a setting where the data-
generating process in each time period is only a function of one unobserved
covariate and one prior-period outcome, with no other observed covariates
entering the outcome model. Specifically, we generate data according to the
following general model:
$Y_{t}(1)=Y_{t}(0)=\beta_{t}f(\theta)+\gamma_{t}g(Y_{t-1}(0))+(t-1)+\epsilon_{t},t=-2,-1,0,1$,
$A=1(\text{expit}(f(\theta))>U)$ where $\theta$, $Y_{-2}$, and $\epsilon_{t}$
are all independent standard normal random variables.

We consider several parameterizations of $(\beta_{t},\gamma_{t},f,g,\alpha)$.
First, consider the possible specification of $\gamma_{t}$, which govern the
association between the current outcome and previous outcomes:

  1. 1.

No prior-period effects: $\gamma_{t}=\gamma=0$

  2. 2.

Constant prior-period effects: $\gamma_{t}=\gamma\neq 0$

  3. 3.

Linear effect growth: $\gamma_{t}-\gamma_{t-1}=c>0$

  4. 4.

Random effect changes: $\gamma_{t}-\gamma_{t-1}=c_{t}\neq 0$

We have the same set of possible specifications for $\beta_{t}$, except that
we do not consider $\beta_{t}=0$. Additionally, we consider linear and
additive specifications (a third nonlinear specification is given in Appendix
E.2). Linear corresponds to $f,g$ being the identity function, and additive
corresponds to $f(x)=\text{sign}(x)\exp(\lvert\sin(x)\rvert)(\lvert
x\rvert>1)+(0.9x+x^{3})(\lvert x\rvert<=1)$ and $g(x)=(x-x^{2})(\lvert
x\rvert<=1.5)+2\text{sign}(x)\log(\lvert x\rvert)(\lvert x\rvert>1.5)$ and
$g(x)=(x-x^{2})(\lvert x\rvert<=1.5)+2\text{sign}(x)\log(\lvert
x\rvert)(\lvert x\rvert>1.5)$.

######  Remark 6.

These simulation settings were chosen to explore a wide range of potential
data-generating processes and not to simulate under the required assumptions
of any of the methods. Thus it is not expected that bias will be 0 for any of
the estimators necessarily in most settings.

####  4.2.2 Results

First consider the results for the linear specification, given in Table 3.
when previous outcomes do not directly effect the current outcome
($\gamma_{t}=\gamma=0$) and the effect of $\theta$ is constant over time
$\beta_{t}=\beta$ (the conventional DiD setting), traditional DiD performs the
best. Even in this setting the gDiD methods are competitive with DiD.
Augmented synthetic controls and the ignorability estimators have RMSEs that
are 2-4 times higher then gDiD. Outside of this setting, the gDiD-1 and gDiD-2
estimators outperform all other estimators by at least a factor of about 2 in
terms of RMSE, with two exceptions.

Results for the additive setting are similar and are given in Table 4. Results
for the nonlinear setting are also generally similar, though there are more
instances of violations of Assumption 3, which degrades the performance of
gDiD and ignorability estimators dramatically. See Appendix E.2.

Table 3: Simulation results showing the good performance of the proposed gDiD estimator in comparison to ignorability and DiD estimators. Performance is measured in terms of bias, root mean squared error (RMSE), confidence interval length (CIL), and 95% confidence interval coverage (“coverage”). Columns indicate whether the effect of the unobserved parameter $\theta$ is constant over time ($\beta_{t}=\beta$) or whether it changes over time ($\beta_{t}-\beta_{t-1}=d$). Results are also separated by the effect of the previous timepoint’s outcome: no effect of prior-period outcome ($\gamma_{t}=\gamma=0$), constant prior-period-outcome effect ($\gamma_{t}=\gamma=c$), and changing prior-period-outcome effect ($\gamma_{t}-\gamma_{t-1}=c)$. Results are aggregated over 500 simulations of the linear specification. | $\beta_{t}=\beta$ | $\beta_{t}-\beta_{t-1}=d$
---|---|---
Estimator | Bias | RMSE | CIL | Coverage | Bias | RMSE | CIL | Coverage
$\gamma_{t}=\gamma=0$
DID | -0.001 | 0.087 | 0.351 | 94.8 | 0.422 | 0.432 | 0.368 | 0.8
gDID-0 | -0.001 | 0.105 | 0.401 | 94 | 0.244 | 0.266 | 0.418 | 38
Ign-1 | 0.456 | 0.464 | 0.349 | 0 | 0.593 | 0.603 | 0.445 | 0
Ign-2 | 0.315 | 0.327 | 0.344 | 5.8 | 0.46 | 0.473 | 0.425 | 1.4
AugSynth | 0.24 | 0.255 | 0.963 | - | 0.425 | 0.439 | 1.701 | -
gDID-1 | 0.011 | 0.13 | 0.496 | 92.8 | -0.066 | 0.168 | 0.616 | 93.4
gDID-2 | 0.005 | 0.126 | 0.486 | 94.6 | -0.132 | 0.197 | 0.577 | 83.6
$\gamma_{t}=\gamma=c$
DID | 0.833 | 0.837 | 0.336 | 0 | 1.625 | 1.63 | 0.512 | 0
gDID-0 | 0.39 | 0.4 | 0.347 | 1.2 | 1.01 | 1.018 | 0.489 | 0
Ign-1 | 0.287 | 0.3 | 0.345 | 10.8 | 0.561 | 0.572 | 0.446 | 0.6
Ign-2 | 0.281 | 0.296 | 0.344 | 12.2 | 0.455 | 0.468 | 0.427 | 2.2
AugSynth | 0.273 | 0.285 | 1.091 | - | 0.434 | 0.448 | 1.736 | -
gDID-1 | -0.096 | 0.158 | 0.467 | 84.4 | -0.202 | 0.242 | 0.507 | 64.4
gDID-2 | -0.081 | 0.157 | 0.483 | 89 | -0.19 | 0.24 | 0.58 | 76.4
$\gamma_{t}-\gamma_{t-1}=c$
DID | 1.38 | 1.385 | 0.468 | 0 | 2.276 | 2.282 | 0.669 | 0
gDID-0 | 0.414 | 0.425 | 0.378 | 1.2 | 1.029 | 1.039 | 0.544 | 0
Ign-1 | 0.253 | 0.268 | 0.342 | 18.8 | 0.47 | 0.483 | 0.434 | 2.4
Ign-2 | 0.255 | 0.27 | 0.343 | 17 | 0.437 | 0.452 | 0.43 | 3.4
AugSynth | 0.252 | 0.265 | 1.008 | - | 0.44 | 0.452 | 1.759 | -
gDID-1 | -0.072 | 0.141 | 0.474 | 89.6 | -0.181 | 0.23 | 0.535 | 72.4
gDID-2 | -0.069 | 0.14 | 0.476 | 89.8 | -0.192 | 0.246 | 0.575 | 74.2
Table 4: Simulation results showing the good performance of the proposed gDiD estimator in comparison to ignorability and DiD estimators. Performance is measured in terms of bias, root mean squared error (RMSE), confidence interval length (CIL), and 95% confidence interval coverage (“coverage”). Columns indicate whether the effect of the unobserved parameter $\theta$ is constant over time ($\beta_{t}=\beta$) or whether it changes over time ($\beta_{t}-\beta_{t-1}=d$). Results are also separated by the effect of the previous timepoint’s outcome: no effect of prior-period outcome ($\gamma_{t}=\gamma=0$), constant prior-period-outcome effect ($\gamma_{t}=\gamma=c$), and changing prior-period-outcome effect ($\gamma_{t}-\gamma_{t-1}=c)$. Results are aggregated over 500 simulations of the additive (non-linear) specification. | $\beta_{t}=\beta$ | $\beta_{t}-\beta_{t-1}=d$
---|---|---
Estimator | Bias | RMSE | CIL | Coverage | Bias | RMSE | CIL | Coverage
$\gamma_{t}=\gamma=0$
DID | 0.019 | 0.093 | 0.357 | 94.6 | 0.387 | 0.398 | 0.375 | 1
gDID-0 | 0.017 | 0.101 | 0.337 | 90.8 | 0.188 | 0.212 | 0.344 | 45
Ign-1 | 0.401 | 0.409 | 0.328 | 0 | 0.493 | 0.503 | 0.406 | 0
Ign-2 | 0.285 | 0.297 | 0.321 | 6.4 | 0.385 | 0.398 | 0.388 | 3.4
AugSynth | 0.228 | 0.241 | 0.911 | - | 0.413 | 0.424 | 1.651 | -
gDID-1 | 0.023 | 0.126 | 0.489 | 94.8 | -0.083 | 0.169 | 0.593 | 93.8
gDID-2 | 0.023 | 0.121 | 0.467 | 94.8 | -0.133 | 0.193 | 0.545 | 81.6
$\gamma_{t}=\gamma=c$
DID | 0.103 | 0.148 | 0.427 | 84 | 0.623 | 0.634 | 0.474 | 0
gDID-0 | 0.218 | 0.242 | 0.376 | 38.4 | 0.573 | 0.584 | 0.414 | 0.4
Ign-1 | 0.368 | 0.378 | 0.372 | 1.8 | 0.631 | 0.641 | 0.474 | 0.2
Ign-2 | 0.337 | 0.35 | 0.366 | 4.8 | 0.572 | 0.585 | 0.464 | 0.4
AugSynth | 0.122 | 0.169 | 0.558 | - | 0.437 | 0.456 | 1.75 | -
gDID-1 | -0.034 | 0.134 | 0.528 | 96 | -0.098 | 0.178 | 0.589 | 89
gDID-2 | -0.055 | 0.138 | 0.497 | 93.2 | -0.077 | 0.168 | 0.577 | 91
$\gamma_{t}-\gamma_{t-1}=c$
DID | 0.283 | 0.307 | 0.451 | 32.6 | 0.771 | 0.781 | 0.51 | 0
gDID-0 | 0.192 | 0.227 | 0.405 | 52.6 | 0.5 | 0.517 | 0.441 | 2
Ign-1 | 0.333 | 0.348 | 0.391 | 7.6 | 0.527 | 0.539 | 0.488 | 0.8
Ign-2 | 0.303 | 0.319 | 0.385 | 12 | 0.498 | 0.511 | 0.477 | 0.6
AugSynth | 0.05 | 0.145 | 0.464 | - | 0.381 | 0.41 | 1.525 | -
gDID-1 | -0.02 | 0.145 | 0.539 | 93.8 | -0.089 | 0.177 | 0.591 | 91.4
gDID-2 | -0.033 | 0.139 | 0.505 | 94.2 | -0.072 | 0.174 | 0.569 | 90.4

##  5 Estimating the effect of health system affiliation on quality of care

Health systems, commonly known to be a group of affiliated health care
providers including hospitals and physician organizations, have comprised an
increasing share of the health care market in recent years[15]. Because health
systems are still relatively new, the effect of health system affiliation on
the quality of care delivered by physician organizations is still uncertain.

We use Medicare Fee-for-service claims data from 2013-2018 to study the effect
of health system affiliation on four measures of the quality of care received
by Medicare beneficiaries. These four quality measures come from the Health
Effectiveness Data and Information Set and are approved by the National
Quality Forum. Health system affiliation was determined by examining the
Medicare Provider Enrollment, Chain, and Ownership System and Internal Revenue
Service Form 990 database. Beneficiaries were attributed to physician
organizations (POs, identified by Tax Identification Numbers) at which they
received the plurality of their primary care visits during the calendar year.

To estimate the effect of affiliating in year $t,t=2015,2016,2017,2018$ on
outcomes in that same year, we let $A_{j}=1$ for POs who were first affiliated
in year $t$ and $A_{j}=0$ for POs who never affiliated during our study
period. This is a clustered data setting (Section 3.2) because beneficiaries
are clustered withiin POs and POs were the units affiliating with health
systems. Beneficiaries were included in the analysis of year $t$ if they were
eligible for the quality measure and were attributed to the same TIN in each
of the following years: $t-2,t-1,t$. For each year $t$, we estimated (12),
taking the outcome in year $t$ to be $Y_{ji1}$, and the outcome in year $t-1$
to be $Y_{ji0}$ for beneficiary $i$ in PO $j$.

We adjusted for the following beneficiary-level characteristics: HCC risk
score (a measure of health status), dual eligibility for Medicaid and Medicare
(a measure of low income), social deprivation index (a measure of
socioeconomic status in the beneficiary’s community), age, race, place of
residence (rural, urban, or small town), and an indicator of Medicare
eligibility due to disability prior to aging into eligibility. In addition to
these beneficiary-level covariates, we adjusted for PO-level aggregates of the
same characteristics, as well as the following PO-level characteristics:
number of attributed beneficiaries, number of physicians, and an indicator of
including both primary care and specialists in the PO. To manage computation
time, we subsampled within each affiliating and non-affiliating TIN, up to
1,000 beneficiaries in each affiliating TIN and up to 50 beneficiaries in each
non-affiliating TIN. We used the multiplier bootstrap outlined in Appendix B
for all standard errors. We also estimated aggregated effects, similar to
(11), combining the estimated effects across years of affiliation to arrive at
a single affiliation effect.

Aggregating across years, we estimate that affiliation has a sizable negative
effect on breast cancer screening but smaller effects on the other outcomes.
We estimate that the breast cancer screening rate is about 1.43 perecentage
points lower in the year of affiliation than it would have been in the absence
of affiliation (95% CI [-2.43, -0.42]). We find that the naive ATT for
affiliation (based on ignorability) in the year of affiliation is 3.71
percentage points (95% CI [2.82, 4.61]), but in the year prior to affiliation,
the ATT was even higher: 5.14 percentage points (95% CI [3.88, 6.40]). For the
other outcomes, we found smaller effects of -0.73 percentage points on
adherence to diabetes medication, -0.06 percentage points on adherence to
hypertension medication, and 0.36 percentage points on receipt of an eye exam
among diabetes patients. More detailed results broken out by year of first
affiliating with a health system are given in Figure 1. Effects on breast
cancer screening appear to have been stronger for POs who affiliated in 2015
and 2016.

Figure 1: Estimates of the effect of health system affiliation in the year of
affiliation (2015-2018) on four health care quality measures. Points
correspond to point estimates and line ranges correspond to 95% multiplier
bootstrap confidence intervals.

##  6 Discussion

We have proposed a very general approach for estimating treatment effects in
longitudinal/panel data settings. Generalized DiD uses a stable bias
assumption that is significantly weaker than either the no unmeasured
confounding assumption underpinning ignorability estimators and significantly
more flexible than the parallel trends assumption underpinning DiD estimators.
Our approach is fully non-parametric and is implemented in an R package
available at https://github.com/mrubinst757/gdid.

Generalized DiD also builds upon the “in-time” placebo tests frequently used
in the synthetic controls literature. Specifically, [2] propose using pre-
treatment data to predict pre-treatment outcomes as a type of validity check
on the primary synthetic controls estimate. These tests produce error
estimates that are essentially equivalent to the bias estimates that we
propose. The key difference is that these tests are not used to correct the
synthetic controls estimate – instead, practitioners simply examine the pre-
treatment errors to heuristically evaluate the validity of the synthetic
controls estimate. There are at least two problems with this procedure
generalized DiD solves: first, the pre-treatment errors may be large,
suggesting that the primary estimate is invalid. Second, reporting the primary
estimate conditional on small pre-treatment errors may lead to inflated Type I
error rates, akin to issues identified with testing parallel trends in DiD
[22]. Generalized DiD solves this problem by essentially building the placebo
tests into the estimation procedure, specifically by subtracting off some
function of these “placebo” errors from the primary estimate. One can
therefore think of generalized DiD as a more principled version of these “in-
time” placebo tests.

On the other hand, the price of the stable bias assumption is a stronger
positivity assumption (Assumption 3), which must be carefully considered when
applying these methods. In settings where the positivity may be nearly
violated, one may remove $\bar{{\bf Y}}_{i0}$ from ${\bf W}_{i0}$ and only
condition on ${\bf X}_{i}$ (as in the gDiD-0 estimator in Section 4), or one
may modify ${\bf W}_{i0}$ to include functions of the pre-period outcomes for
which there is better overlap, like the trends in the pre-period outcomes
$Y_{it}-Y_{it-1}$ rather than the outcomes themselves (often done in the
synthetic control literature[1]).

The connection between our estimator and typical ignorability estimators opens
up a number of potentially interesting avenues for future research into
potential violations of gDiD’s required assumptions. Another potential way of
addressing positivity violations is to use targeted maximum likelihood [28,
26], which has been shown previously to handle positivity violations well
[20], instead of our approach based on augmented inverse probability
weighting. Also, violations of the stable bias assumption could potentially be
evaluated by adapting methods for omitted variable bias of causal effects
under ignorability [11].

## References

  * [1] Alberto Abadie.  Using synthetic controls: Feasibility, data requirements, and methodological aspects.  Journal of Economic Literature, 59(2):391–425, 2021.
  * [2] Alberto Abadie, Alexis Diamond, and Jens Hainmueller.  Comparative politics and the synthetic control method.  American Journal of Political Science, 59(2):495–510, 2015.
  * [3] Heejung Bang and James M Robins.  Doubly robust estimation in missing data and causal inference models.  Biometrics, 61(4):962–973, 2005.
  * [4] Eli Ben-Michael.  augsynth: The Augmented Synthetic Control Method, 2023.  R package version 0.2.0.
  * [5] Eli Ben-Michael, Avi Feller, and Jesse Rothstein.  The augmented synthetic control method.  Journal of the American Statistical Association, 116(536):1789–1803, 2021.
  * [6] Kirill Borusyak, Xavier Jaravel, and Jann Spiess.  Revisiting event study designs: Robust and efficient estimation.  arXiv preprint arXiv:2108.12419, 2021.
  * [7] David Bruns-Smith, Oliver Dukes, Avi Feller, and Elizabeth L Ogburn.  Augmented balancing weights as undersmoothed regressions.  Cited on, page 3, 2023.
  * [8] Brantly Callaway and Pedro HC Sant’Anna.  Difference-in-differences with multiple time periods.  Journal of Econometrics, 225(2):200–230, 2021.
  * [9] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, and Whitney Newey.  Double/debiased/neyman machine learning of treatment effects.  American Economic Review, 107(5):261–265, 2017.
  * [10] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins.  Double/debiased machine learning for treatment and structural parameters, 2018.
  * [11] Victor Chernozhukov, Carlos Cinelli, Whitney Newey, Amit Sharma, and Vasilis Syrgkanis.  Long story short: Omitted variable bias in causal machine learning.  Technical report, National Bureau of Economic Research, 2022.
  * [12] Victor Chernozhukov, Kaspar Wüthrich, and Yinchu Zhu.  An exact and robust conformal inference method for counterfactual and synthetic controls.  Journal of the American Statistical Association, 116(536):1849–1864, 2021.
  * [13] Jamie R Daw and Laura A Hatfield.  Matching and regression to the mean in difference-in-differences analysis.  Health services research, 53(6):4138–4156, 2018.
  * [14] Max H Farrell, Tengyuan Liang, and Sanjog Misra.  Deep neural networks for estimation and inference.  Econometrica, 89(1):181–213, 2021.
  * [15] Michael F Furukawa, Laura Kimmey, David J Jones, Rachel M Machta, Jing Guo, and Eugene C Rich.  Consolidation of providers into health systems increased substantially, 2016–18: Study examines provider consolidation into vertically-integrated health systems.  Health Affairs, 39(8):1321–1325, 2020.
  * [16] Dae Woong Ham and Luke Miratrix.  Benefits and costs of matching prior to a difference in differences analysis when parallel trends does not hold.  arXiv preprint arXiv:2205.08644, 2022.
  * [17] Joseph DY Kang and Joseph L Schafer.  Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data.  2007\.
  * [18] Edward H Kennedy.  Semiparametric doubly robust targeted double machine learning: a review.  arXiv preprint arXiv:2203.06469, 2022.
  * [19] Eric Polley, Erin LeDell, Chris Kennedy, Sam Lendle, and Mark van der Laan.  Package ‘superlearner’.  CRAN, 2019.
  * [20] Kristin E Porter, Susan Gruber, Mark J Van Der Laan, and Jasjeet S Sekhon.  The relative performance of targeted maximum likelihood estimators.  The international journal of biostatistics, 7(1):0000102202155746791308, 2011.
  * [21] Ashesh Rambachan and Jonathan Roth.  Design-based uncertainty for quasi-experiments.  arXiv preprint arXiv:2008.00602, 2020.
  * [22] Jonathan Roth.  Pretest with caution: Event-study estimates after testing for parallel trends.  American Economic Review: Insights, 4(3):305–322, 2022.
  * [23] Jonathan Roth, Pedro HC Sant’Anna, Alyssa Bilinski, and John Poe.  What’s trending in difference-in-differences? a synthesis of the recent econometrics literature.  Journal of Econometrics, 2023.
  * [24] Donald B Rubin.  Causal inference using potential outcomes: Design, modeling, decisions.  Journal of the American Statistical Association, 100(469):322–331, 2005.
  * [25] Pedro HC Sant’Anna and Jun Zhao.  Doubly robust difference-in-differences estimators.  Journal of Econometrics, 219(1):101–122, 2020.
  * [26] Megan S Schuler and Sherri Rose.  Targeted maximum likelihood estimation for causal inference in observational studies.  American journal of epidemiology, 185(1):65–73, 2017.
  * [27] Mark J Van der Laan, Eric C Polley, and Alan E Hubbard.  Super learner.  Statistical applications in genetics and molecular biology, 6(1), 2007.
  * [28] Mark J Van Der Laan and Daniel Rubin.  Targeted maximum likelihood learning.  The international journal of biostatistics, 2(1), 2006.
  * [29] Aad W Van der Vaart.  Asymptotic statistics, volume 3.  Cambridge university press, 2000.
  * [30] Stefan Wager and Guenther Walther.  Adaptive concentration of regression trees, with application to random forests.  arXiv preprint arXiv:1503.06388, 2015.

##  Appendix A Proof of Proposition 1 (main result)

The proof proceeds by showing that the proposed estimator falls into the class
of double machine learning estimators of [10]. First, we require a set of
regularity conditions on the data-generating process. Define positive
constants $C^{*},c,\epsilon>0$ and $\iota\in(0,1)$.

######  Assumption 15.

We have the following for $t=0,1$:

  1. 1.

$\left(E|Y_{it}|^{q}\right)^{1/q}\leq C^{*}$.

  2. 2.

$\left\\{E|Y_{it}-\mu_{t}({\bf W}{{}_{i}t-1})|^{2}\right\\}^{1/2}>c$.

  3. 3.

$\max E\left[\\{Y_{it}-\mu_{t}({\bf W}_{it-1})\\}^{2}|{\bf
W}_{it-1}\right]\leq C$.

  4. 4.

$n^{-1}\sum_{i=1}^{n}A_{i}\longrightarrow\iota$.

And we require a few further conditions on the nuisance function estimation

######  Assumption 16.

For nuisance functions $\eta({\bf W}_{i0})=\\{\mu_{1}({\bf
W}_{i0}),\mu_{0}({\bf W}_{i-1}),\pi({\bf W}_{i0}),\pi_{0}({\bf W}_{i-1})\\}$
and estimated versions ${\widehat{\eta}}({\bf
W}_{i0})=\\{{\widehat{\mu}}_{1}({\bf W}_{i0}),{\widehat{\mu}}_{0}({\bf
W}_{i-1}),{\widehat{\pi}}({\bf W}_{i0}),{\widehat{\pi}}_{0}({\bf
W}_{i-1})\\}$, we have $\left\\{E|{\widehat{\eta}}({\bf W}_{i0})-\eta({\bf
W}_{i0})|^{q}\right\\}^{1/q}<C^{*},\left\\{E|{\widehat{\eta}}({\bf
W}_{i0})-\eta({\bf W}_{i0})|^{2}\right\\}^{1/2}=o_{p}(1),$ and $\max
E\\{{\widehat{\pi}}({\bf W}_{i0})-1/2\\}\leq 1/2-\epsilon,\max
E\\{{\widehat{\pi}}_{0}({\bf W}_{i-1})-1/2\\}\leq 1/2-\epsilon$.

Under these conditions, we have that the two components of ${\widehat{\tau}}$
are estimators of the form of [10], which implies that ${\widehat{\tau}}$ also
falls into the same class of estimators:

$\displaystyle{\widehat{\tau}}$
$\displaystyle={\widehat{\tau}}^{\text{ign.}}_{1}-{\widehat{\tau}}^{\text{ign.}}_{0}$
$\displaystyle{\widehat{\tau}}^{\text{ign.}}_{1}$
$\displaystyle=n_{1}^{-1}\sum_{i=1}^{n}\phi^{\text{ign.}}_{1}({\bf
O}_{i},\tau^{\text{ign.}},{\widehat{\eta}})+\tau^{\text{ign.}}=n_{1}^{-1}\sum_{i=1}^{n}\left[Y_{i1}A_{i}-\frac{(1-A_{i}){\widehat{\pi}}({\bf
W}_{i0})Y_{i1}+\\{A_{i}-{\widehat{\pi}}({\bf
W}_{i0})\\}{\widehat{\mu}}_{1}({\bf W}_{i0})}{1-{\widehat{\pi}}({\bf
W}_{i0})}\right]$ $\displaystyle{\widehat{\tau}}^{\text{ign.}}_{0}$
$\displaystyle=n_{1}^{-1}\sum_{i=1}^{n}\phi^{\text{ign.}}_{0}({\bf
O}_{i},\tau^{\text{ign.}}_{0},{\widehat{\eta}})+\tau^{\text{ign.}}_{0}=n_{1}^{-1}\sum_{i=1}^{n}\left[Y_{i0}A_{i}+\frac{(1-A_{i}){\widehat{\pi}}_{0}({\bf
W}_{i-1})Y_{i0}+\\{A_{i}-{\widehat{\pi}}_{0}({\bf
W}_{i-1})\\}{\widehat{\mu}}_{0}({\bf W}_{i-1})}{1-{\widehat{\pi}}_{0}({\bf
W}_{i-1})}\right]$

where $\tau^{\text{ign.}}=E(Y_{i1}|A_{i}=1)-E\left\\{E(Y_{i1}|{\bf
W}_{i0},A_{i}=0)|A_{i}=1\right\\}$, and
$\tau^{\text{ign.}}_{0}=E(Y_{i0}|A_{i}=1)-E\left\\{E(Y_{i0}|{\bf
W}_{i-1},A_{i}=0)|A_{i}=1\right\\}$. These estimators satisfy Neyman
orthogonality (see Definition 2.1 in [10]) such that

$\displaystyle E\left\\{\phi^{\text{ign.}}_{1}({\bf
O},\tau^{\text{ign.}},\eta)\right\\}=0$ $\displaystyle
E\left\\{\phi^{\text{ign.}}_{0}({\bf
O},\tau^{\text{ign.}}_{0},\eta)\right\\}=0$
$\displaystyle\partial_{\widetilde{\eta}}E\left\\{\phi^{\text{ign.}}_{1}({\bf
O},\tau^{\text{ign.}},\eta)\right\\}[{\widetilde{\eta}}-\eta]=0$
$\displaystyle\partial_{\widetilde{\eta}}E\left\\{\phi^{\text{ign.}}_{0}({\bf
O},\tau^{\text{ign.}}_{0},\eta)\right\\}[{\widetilde{\eta}}-\eta]=0\qquad\text{for
all }{\widetilde{\eta}}\in{\cal T}_{{\widetilde{\eta}}}$

where ${\cal T}_{{\widetilde{\eta}}}$ is a properly shrinking neighborhood of
the true nuisance functions $\eta$ and
$\partial_{\eta}E\left\\{\phi^{\text{ign.}}_{1}({\bf
O},\tau^{\text{ign.}},\eta)\right\\}[{\widetilde{\eta}}-\eta]$ is the limit of
the pathwise derivative

$\displaystyle\partial_{\widetilde{\eta}}E\left\\{\phi({\bf
O},\tau,\eta)\right\\}[{\widetilde{\eta}}-\eta]$
$\displaystyle=\lim_{r\rightarrow 0}\partial_{r}\left[E\left\\{\phi({\bf
O},\tau,\eta+r({\widetilde{\eta}}-\eta))\right\\}\right].$

Thus, because
${\widehat{\tau}}={\widehat{\tau}}^{\text{ign.}}_{1}-{\widehat{\tau}}^{\text{ign.}}_{0}$
and the influence function for $\tau$, defined in (3) is the difference
between these two $\phi({\bf O})=\phi({\bf
O},\tau,\eta)=\phi^{\text{ign.}}_{1}({\bf
O},\tau^{\text{ign.}},\eta)-\phi^{\text{ign.}}_{0}({\bf
O},\tau^{\text{ign.}}_{0},\eta)$, we obtain Neyman orthogonality for the gDiD
estimator

$\displaystyle E\left\\{\phi({\bf O},\tau,\eta)\right\\}=0$
$\displaystyle\partial_{\eta}E\left\\{\phi({\bf
O},\tau,\eta)\right\\}[{\widetilde{\eta}}-\eta]=0.$

This, along with the regularity conditions in Assumptions 15 and 16, ensures
the results in the main text.

##  Appendix B Multiplier bootstrap variance estimation

In the text, we give plug-in estimators for the variance of
${\widehat{\tau}}$. In this section, we review the use of multiplier bootstrap
estimators for the variance.

###  B.1 Standard iid setting

The multiplier bootstrap proceeds by generating $n$ independent random
variables ${\bf G}=(G_{i})_{i=1,...,n}$ such that $E(G_{i})=1$ and
$\mbox{var}(G_{i})=1$. Then one bootstrap replicate can be computed as

$\displaystyle{\widehat{\tau}}^{*}_{\bf
G}=n_{1}^{-1}\sum_{i=1}^{n}G_{i}{\widehat{\phi}}({\bf O}_{i}).$

It is easy to see that $E({\widehat{\tau}}^{*}_{\bf G}|{\bf
O})=n_{1}^{-1}\sum_{i=1}^{n}{\widehat{\phi}}({\bf O}_{i})={\widehat{\tau}}$,
and $\mbox{var}({\widehat{\tau}}^{*}_{\bf G}|{\bf
O})=n_{1}^{-2}\sum_{i=1}^{n}{\widehat{\phi}}({\bf O}_{i})^{2}$. So,
$\mbox{var}\\{n^{\frac{1}{2}}({\widehat{\tau}}_{\bf
G}^{*}-{\widehat{\tau}})\\}$ converges to $\sigma^{2}$ conditional on ${\bf
O}$. Therefore, the multiplier bootstrap may be used as an alternative to the
plug-in estimator (8). Operationally, this means generating a large number of
replications $B$ of ${\bf G}:\\{{\bf G}_{b}\\}_{b=1,..,B}$, yielding a large
set of estimates $\\{{\widehat{\tau}}^{*}_{{\bf G}_{b}}\\}_{1,...,B}$. One may
then construct a confidence interval for $\tau$ based on the quantiles of
$\\{{\widehat{\tau}}^{*}_{{\bf G}_{b}}\\}_{1,...,B}$ or by estimating
$\sigma^{2}$ as $\mbox{var}\\{{\widehat{\tau}}^{*}_{{\bf
G}_{b}}\\}_{1,...,B}$.

###  B.2 Clustered data setting

The multiplier bootstrap in the clustered data setting is very similar as in
the iid setting, except that we now construct $C$ independent random variables
${\bf G}=(G_{j})_{j=1,...,C}$ and construct the bootstrap estimators as

$\displaystyle{\widehat{\tau}}_{\bf
G}=C^{-1}\sum_{j=1}^{C}G_{j}{\widetilde{\phi}}({\bf O}_{j}).$

It is again easy to see that $E({\widehat{\tau}}^{*}_{\bf G}|{\bf
O})=C^{-1}\sum_{i=1}^{n}{\widehat{\phi}}({\bf O}_{j})={\widehat{\tau}}$, and
$\mbox{var}({\widehat{\tau}}^{*}_{\bf G}|{\bf
O})=C^{-2}\sum_{j=1}^{C}{\widetilde{\phi}}({\bf O}_{j})^{2}$. It is again
straightforward to show that $C^{-1/2}({\widehat{\tau}}^{*}_{\bf
G}-{\widehat{\tau}})|{\bf O}\longrightarrow N(0,\zeta^{2}).$

##  Appendix C Parametric nuisance function estimation

As noted in the text, in some cases we may wish to estimate all nuisance
functions using parametric models (e.g. generalized linear models) and avoid
data splitting altogether (e.g., if sample size is too small to feasibly allow
sample splitting). In this case the asymptotic variance of the estimator
depends on whether the models are correctly specified or not. However,
regardless of whether the models were correctly specified, the asymptotic
variance has a simple form assuming that the nuisance parameters and estimator
were the solutions to estimating equations taking the form:

$\displaystyle\mathbb{P}_{n}\begin{bmatrix}m({\bf
O};{\widehat{\boldsymbol{\beta}}})\\\ \phi({\bf
O},{\widehat{\tau}},{\widehat{\boldsymbol{\beta}}}),\end{bmatrix}=\mathbb{P}_{n}\left\\{\varphi({\bf
O};{\widehat{\boldsymbol{\theta}}})\right\\}=0$

where $\mathbb{P}_{n}\\{f({\bf O})\\}=n^{-1}\sum_{i=1}^{n}f({\bf O}_{i})$,
${\widehat{\boldsymbol{\beta}}}$ is a finite-dimensional parameter, $m({\bf
O},\boldsymbol{\beta})$ is the estimating equation for the nuisance functions,
$\varphi=(m,\phi)$, and
${\widehat{\boldsymbol{\theta}}}=({\widehat{\tau}},{\widehat{\boldsymbol{\beta}}})$.
For example, in the two time-period case, $m$ may be an estimating function,
such as a score equation, for a series of models for
$\eta=(\mu_{1},\mu_{0},\pi,\pi_{0})$ parameterized by the finite-dimensional
vector $\boldsymbol{\beta}$. And $\phi$ is the estimated efficient influence
function (3) that plugs in nuisance function estimates parameterized by
${\widehat{\boldsymbol{\beta}}}$. Then, by Theorem 5.41 in van der Vaart
([29]):

$\displaystyle{\widehat{\boldsymbol{\theta}}}-\boldsymbol{\theta}$
$\displaystyle=\begin{bmatrix}{\widehat{\boldsymbol{\beta}}}\\\
{\widehat{\tau}}\end{bmatrix}-\begin{bmatrix}\boldsymbol{\beta}\\\
\tau\end{bmatrix}$
$\displaystyle=\mathbb{P}_{n}\left\\{-\left(\frac{\partial\mathbb{E}[\varphi({\bf
O};\boldsymbol{\theta})]}{\partial\boldsymbol{\theta}}\right)^{-1}\varphi({\bf
O};\boldsymbol{\theta})\right\\}+o_{p}(1/\sqrt{n})$

and, therefore, the asymptotic variance of ${\widehat{\boldsymbol{\theta}}}$
is given by:

$\displaystyle\left(\frac{\partial\mathbb{E}[\varphi({\bf
O};\boldsymbol{\theta})]}{\partial\boldsymbol{\theta}}\right)^{-1}\mathbb{E}[\varphi({\bf
O};\boldsymbol{\theta})\varphi({\bf
O};\boldsymbol{\theta})^{\top}]\left(\frac{\partial\mathbb{E}[\varphi({\bf
O};\boldsymbol{\theta})]}{\partial\boldsymbol{\theta}}\right)^{-1}$

The variance can be estimated by plugging-in the empirical counterparts (i.e.
the sandwich estimator), and the variance estimate of $\hat{\tau}$ is given by
the final diagonal element of this estimated matrix. We refer to [29] for more
technical details.

##  Appendix D Proof of Proposition 2 (clustered treatment assignment)

The primary difference between the proof of Proposition 2 and Proposition 1 is
the presence of clustering. We thus first establish that the number of
clusters grows with sample size:

######  Assumption 17.

$C\rightarrow\infty$.

Next, we restate the identifiability conditions (Assumptions 1, 2, 3, and 4)
in the clustered data setting.

######  Assumption 18 (Consistency, clustered data setting).

$Y_{ji1}(a)=Y_{ji1}\text{ when }A_{j}=a,a=0,1.$

######  Assumption 19 (No anticipation, clustered data setting).

$Y_{jit}(a)=Y_{jit},a=0,1;t<1.$

######  Assumption 20 (Positivity, clustered data setting).

$P\left\\{p({\widetilde{\bf W}}_{j0})<1-\epsilon\right\\}=1$ for some
$\epsilon>0$.

######  Assumption 21 (Stable bias, clustered data setting).

$\displaystyle E\left[E\left\\{Y_{ji1}(0)|{\bf
V}_{ji0},A_{i}=1\right\\}-E\left\\{Y_{ji1}|{\bf
V}_{ji0},A_{i}=0\right\\}|A_{i}=1\right]=E\left[E\left\\{Y_{ji0}|{\bf
V}_{ji-1},A_{i}=1\right\\}-E\left\\{Y_{ji0}|{\bf
V}_{ji-1},A_{i}=0\right\\}|A_{i}=1\right].$

These are mere extensions to the clustered data setting and do not introduce
any notable additional complexity. Next, define cluster-level versions of the
outcomes and mean functions

$\displaystyle Z_{jt}=\frac{C}{N_{1}}\sum_{i=1}^{n_{j}}Y_{jit}$
$\displaystyle\bar{m}_{jt}=\frac{C}{N_{1}}\sum_{i=1}^{n_{j}}m_{t}({\bf
V}_{jit-1})$

where $N_{1}=\sum_{j=1}^{C}n_{j}A_{j}$ is the total number of treatment
observations. Then we can rewrite ${\widehat{\tau}}$ in terms of these
aggregate values

$\displaystyle{\widehat{\tau}}$
$\displaystyle=C^{-1}\sum_{j=1}^{C}\left[Z_{j1}A_{j}-\frac{(1-A_{j}){\widehat{p}}({\widetilde{\bf
W}}_{j0})Z_{j1}+\\{A_{j}-{\widehat{p}}({\widetilde{\bf
W}}_{j0})\\}\widehat{\bar{m}}_{1}}{1-{\widehat{p}}({\widetilde{\bf
W}}_{j0})}-\right.$
$\displaystyle\qquad\qquad\left.Z_{j0}A_{j}+\frac{(1-A_{j}){\widehat{p}}_{0}({\widetilde{\bf
W}}_{j-1})Z_{j0}+\\{A_{j}-{\widehat{p}}_{0}({\widetilde{\bf
W}}_{j-1})\\}\widehat{\bar{m}}_{0}}{1-{\widehat{p}}_{0}({\widetilde{\bf
W}}_{j-1})}\right]$ $\displaystyle=C^{-1}\sum_{j=1}^{C}{\widetilde{\phi}}({\bf
O}_{j}).$

We may now analyze ${\widehat{\tau}}$ in the same manner as for Proposition 1
in Appendix A. We first establish the primitive regularity conditions

######  Assumption 22.

We have the following for $t=0,1$:

  1. 1.

$\left(E|Z_{jt}|^{q}\right)^{1/q}\leq C^{*}$.

  2. 2.

$\left\\{E|Z_{jt}-\bar{m}_{t}|^{2}\right\\}^{1/2}>c$.

  3. 3.

$\max E\left[\\{Z_{it}-\bar{m}_{t}\\}^{2}|{\bf V}_{jit-1}\right]\leq C$.

  4. 4.

$C^{-1}\sum_{i=1}^{n}A_{j}\longrightarrow\iota$.

And we adapt conditions on the nuisance function estimation

######  Assumption 23 (Nuisance function convergence, clustered data setting).

$\displaystyle
E\left[\left\\{\widehat{\bar{m}}_{1}-\bar{m}_{1}\right\\}^{2}\right]^{1/2}\times
E\left[\left\\{{\widehat{p}}({\widetilde{\bf W}}_{j0})-\pi({\widetilde{\bf
W}}_{j0}\right\\}^{2}\right]^{1/2}=o_{p}(C^{-1/2})$ $\displaystyle
E\left[\left\\{\widehat{\bar{m}}_{0}-\bar{m}_{0}\right\\}^{2}\right]^{1/2}\times
E\left[\left\\{{\widehat{\pi}}_{0}({\bf W}_{i-1})-\pi_{0}({\bf
W}_{i-1}\right\\}^{2}\right]^{1/2}=o_{p}(C^{-1/2})$

######  Assumption 24.

For nuisance functions $\eta({\bf
V}_{ji0})=\\{\bar{m}_{1},\bar{m}_{0},p({\widetilde{\bf
W}}_{j0}),p_{0}({\widetilde{\bf W}}_{j-1})\\}$ and estimated versions
${\widehat{\eta}}({\bf
V}_{ji0})=\\{\widehat{\bar{m}}_{1},\widehat{\bar{m}}_{0},{\widehat{p}}({\widetilde{\bf
W}}_{j0}),{\widehat{p}}_{0}({\widetilde{\bf W}}_{j-1})\\}$, we have

  * •

$\left\\{E|{\widehat{\eta}}({\bf V}_{ji0})-\eta({\bf
V}_{ji0})|^{q}\right\\}^{1/q}<C^{*}$

  * •

$\left\\{E|{\widehat{\eta}}({\bf V}_{ji0})-\eta({\bf
V}_{ji0})|^{2}\right\\}^{1/2}=o_{p}(1)$

  * •

$\max E\\{{\widehat{p}}({\widetilde{\bf W}}_{j0})-1/2\\}\leq 1/2-\epsilon,\max
E\\{{\widehat{p}}_{0}({\widetilde{\bf W}}_{j-1})-1/2\\}\leq 1/2-\epsilon$

,

The argument proceeds just as in Appendix A.

##  Appendix E Complete simulation results

###  E.1 Simulation one

We present additional results when setting $\zeta=0.3,0.5$. In these settings
all ignorability-based estimators and gDiD perform poorly due to limited
overlap between the pre-treatment outcome levels in the treatment and control
groups.

Table 5: Simulation study one results: $\zeta=0.3$ | Linear | Nonlinear
---|---|---
Estimator | Bias | RMSE | CIL | Coverage | Bias | RMSE | CIL | Coverage
cDiD | -0.176 | 0.273 | 0.784 | 86.2 | -0.228 | 0.329 | 0.889 | 83
gDiD-0 | -0.179 | 0.272 | 0.791 | 85.7 | -0.218 | 0.316 | 0.872 | 83.3
DiD | -0.197 | 0.304 | 0.963 | 88.7 | -0.209 | 0.319 | 0.964 | 86.7
gDiD-2 | 0.4 | 3.095 | 2.354 | 46.3 | 0.542 | 1.063 | 1.494 | 39.5
gDiD-1 | 0.406 | 0.57 | 1.233 | 67 | 0.39 | 0.502 | 1.003 | 59.2
Ign-1 | -4.616 | 4.622 | 1.09 | 0 | -4.629 | 4.632 | 0.997 | 0
Ign-2 | -7.245 | 7.71 | 2.373 | 0.7 | -7.134 | 7.166 | 1.752 | 0.3
AugSynth | -9.854 | 9.896 | - | - | -9.772 | 9.819 | - | -
Table 6: Simulation study one results: $\zeta=0.5$ | Linear | Nonlinear
---|---|---
Estimator | Bias | RMSE | CIL | Coverage | Bias | RMSE | CIL | Coverage
cDiD | -0.189 | 0.28 | 0.786 | 84.2 | -0.231 | 0.337 | 0.896 | 80.3
gDiD-0 | -0.19 | 0.281 | 0.792 | 84.2 | -0.226 | 0.324 | 0.874 | 78.7
DiD | -0.211 | 0.329 | 0.965 | 84.7 | -0.214 | 0.328 | 0.964 | 85.7
gDiD-1 | 0.632 | 0.911 | 1.48 | 50 | 0.606 | 0.833 | 1.359 | 46.2
Ign-1 | -7.644 | 7.653 | 1.661 | 0 | -7.686 | 7.696 | 1.639 | 0
AugSynth | -16.301 | 16.328 | - | - | -16.26 | 16.29 | - | -
Ign-2 | -11310059838.405 | 800877093417.37 | 427364180154.257 | 3 | 303.724 | 6421.427 | 1264.803 | 2.2
gDiD-2 | -38559885534.194 | 2106347444114.77 | 1105410193452.54 | 39 | -70686491808.493 | 2145878599062.18 | 1308869531810.81 | 34.9

###  E.2 Simulation two

The table below presents all simulation results from all specifications
explored.

| Linear | Additive | Nonlinear
---|---|---|---
Estimator | Bias | RMSE | CIL | Coverage | Bias | RMSE | CIL | Coverage | Bias | RMSE | CIL | Coverage
$\gamma=0$: $\beta_{t}=\beta$
AugSynth | 0.24 | 0.255 | 0.963 | 100 | - | - | - | - | 0.228 | 0.241 | 0.911 | 100
DID | -0.001 | 0.087 | 0.351 | 94.8 | - | - | - | - | 0.019 | 0.093 | 0.357 | 94.6
Ign-1 | 0.456 | 0.464 | 0.349 | 0 | - | - | - | - | 0.401 | 0.409 | 0.328 | 0
Ign-2 | 0.315 | 0.327 | 0.344 | 5.8 | - | - | - | - | 0.285 | 0.297 | 0.321 | 6.4
gDID-1 | 0.011 | 0.13 | 0.496 | 92.8 | - | - | - | - | 0.023 | 0.126 | 0.489 | 94.8
gDID-2 | 0.005 | 0.126 | 0.486 | 94.6 | - | - | - | - | 0.023 | 0.121 | 0.467 | 94.8
gDiD-0 | -0.001 | 0.105 | 0.401 | 94 | - | - | - | - | 0.017 | 0.101 | 0.337 | 90.8
sDID | -0.002 | 0.159 | 0.605 | 94.2 | - | - | - | - | 0.02 | 0.157 | 0.618 | 94.6
$\gamma=0$: $\beta_{t}-\beta_{t-1}=d$
AugSynth | 0.425 | 0.439 | 1.701 | 100 | - | - | - | - | 0.413 | 0.424 | 1.651 | 100
DID | 0.422 | 0.432 | 0.368 | 0.8 | - | - | - | - | 0.387 | 0.398 | 0.375 | 1
Ign-1 | 0.593 | 0.603 | 0.445 | 0 | - | - | - | - | 0.493 | 0.503 | 0.406 | 0
Ign-2 | 0.46 | 0.473 | 0.425 | 1.4 | - | - | - | - | 0.385 | 0.398 | 0.388 | 3.4
gDID-1 | -0.066 | 0.168 | 0.616 | 93.4 | - | - | - | - | -0.083 | 0.169 | 0.593 | 93.8
gDID-2 | -0.132 | 0.197 | 0.577 | 83.6 | - | - | - | - | -0.133 | 0.193 | 0.545 | 81.6
gDiD-0 | 0.244 | 0.266 | 0.418 | 38 | - | - | - | - | 0.188 | 0.212 | 0.344 | 45
sDID | 0.019 | 0.148 | 0.606 | 96 | - | - | - | - | 0.001 | 0.156 | 0.618 | 96
$\gamma_{t}=\gamma$: $\beta_{t}=\beta$
AugSynth | 0.273 | 0.285 | 1.091 | 100 | 0.122 | 0.169 | 0.558 | 100 | 0.49 | 0.562 | 1.982 | 100
DID | 0.833 | 0.837 | 0.336 | 0 | 0.103 | 0.148 | 0.427 | 84 | 0.398 | 0.422 | 0.547 | 17.4
Ign-1 | 0.287 | 0.3 | 0.345 | 10.8 | 0.368 | 0.378 | 0.372 | 1.8 | 0.132 | 0.733 | 1.093 | 59.8
Ign-2 | 0.281 | 0.296 | 0.344 | 12.2 | 0.337 | 0.35 | 0.366 | 4.8 | -1.766 | 40.604 | 8.424 | 66
gDID-1 | -0.096 | 0.158 | 0.467 | 84.4 | -0.034 | 0.134 | 0.528 | 96 | -0.131 | 0.709 | 1.066 | 95
gDID-2 | -0.081 | 0.157 | 0.483 | 89 | -0.055 | 0.138 | 0.497 | 93.2 | -1.969 | 40.105 | 8.262 | 94.4
gDiD-0 | 0.39 | 0.4 | 0.347 | 1.2 | 0.218 | 0.242 | 0.376 | 38.4 | 0.344 | 0.507 | 0.684 | 32.4
sDID | 0.01 | 0.096 | 0.351 | 93.4 | -0.122 | 0.205 | 0.647 | 88 | 0.038 | 0.169 | 0.66 | 95
$\gamma_{t}=\gamma$: $\beta_{t}-\beta_{t-1}=d$
AugSynth | 0.434 | 0.448 | 1.736 | 100 | 0.437 | 0.456 | 1.75 | 100 | 0.609 | 0.692 | 2.452 | 100
DID | 1.625 | 1.63 | 0.512 | 0 | 0.623 | 0.634 | 0.474 | 0 | 0.716 | 0.733 | 0.638 | 0.4
Ign-1 | 0.561 | 0.572 | 0.446 | 0.6 | 0.631 | 0.641 | 0.474 | 0.2 | 0.088 | 4.598 | 2.06 | 36.6
Ign-2 | 0.455 | 0.468 | 0.427 | 2.2 | 0.572 | 0.585 | 0.464 | 0.4 | 0.163 | 1.591 | 1.624 | 46.6
gDID-1 | -0.202 | 0.242 | 0.507 | 64.4 | -0.098 | 0.178 | 0.589 | 89 | -0.314 | 4.506 | 1.999 | 96.2
gDID-2 | -0.19 | 0.24 | 0.58 | 76.4 | -0.077 | 0.168 | 0.577 | 91 | -0.224 | 1.562 | 1.527 | 95.6
gDiD-0 | 1.01 | 1.018 | 0.489 | 0 | 0.573 | 0.584 | 0.414 | 0.4 | 0.682 | 0.736 | 0.652 | 7.1
sDID | 0.396 | 0.406 | 0.366 | 1.2 | -0.081 | 0.179 | 0.656 | 93 | 0.195 | 0.252 | 0.688 | 83.4
$\gamma_{t}-\gamma_{t-1}=c$: $\beta_{t}=\beta$
AugSynth | 0.252 | 0.265 | 1.008 | 100 | 0.05 | 0.145 | 0.464 | 100 | 0.522 | 0.614 | 2.116 | 100
DID | 1.38 | 1.385 | 0.468 | 0 | 0.283 | 0.307 | 0.451 | 32.6 | 0.654 | 0.677 | 0.697 | 3.2
Ign-1 | 0.253 | 0.268 | 0.342 | 18.8 | 0.333 | 0.348 | 0.391 | 7.6 | 0.027 | 2.392 | 1.52 | 71.2
Ign-2 | 0.255 | 0.27 | 0.343 | 17 | 0.303 | 0.319 | 0.385 | 12 | 0.027 | 1.793 | 1.433 | 76.8
gDID-1 | -0.072 | 0.141 | 0.474 | 89.6 | -0.02 | 0.145 | 0.539 | 93.8 | -0.236 | 2.366 | 1.491 | 93.8
gDID-2 | -0.069 | 0.14 | 0.476 | 89.8 | -0.033 | 0.139 | 0.505 | 94.2 | -0.177 | 1.77 | 1.347 | 93
gDiD-0 | 0.414 | 0.425 | 0.378 | 1.2 | 0.192 | 0.227 | 0.405 | 52.6 | 0.447 | 0.536 | 0.735 | 35.4
sDID | 0.57 | 0.577 | 0.359 | 0 | -0.113 | 0.201 | 0.649 | 90 | 0.231 | 0.291 | 0.695 | 74.8
$\gamma_{t}-\gamma_{t-1}=c$: $\beta_{t}-\beta_{t-1}=d$
AugSynth | 0.44 | 0.452 | 1.759 | 100 | 0.381 | 0.41 | 1.525 | 100 | 0.584 | 0.696 | 2.37 | 100
DID | 2.276 | 2.282 | 0.669 | 0 | 0.771 | 0.781 | 0.51 | 0 | 0.957 | 0.979 | 0.8 | 0.2
Ign-1 | 0.47 | 0.483 | 0.434 | 2.4 | 0.527 | 0.539 | 0.488 | 0.8 | 0.292 | 0.563 | 1.065 | 52.4
Ign-2 | 0.437 | 0.452 | 0.43 | 3.4 | 0.498 | 0.511 | 0.477 | 0.6 | 0.231 | 0.442 | 1.113 | 66.4
gDID-1 | -0.181 | 0.23 | 0.535 | 72.4 | -0.089 | 0.177 | 0.591 | 91.4 | -0.108 | 0.472 | 1.045 | 96
gDID-2 | -0.192 | 0.246 | 0.575 | 74.2 | -0.072 | 0.174 | 0.569 | 90.4 | -0.149 | 0.387 | 1.027 | 95
gDiD-0 | 1.029 | 1.039 | 0.544 | 0 | 0.5 | 0.517 | 0.441 | 2 | 0.737 | 0.805 | 0.771 | 12.3
sDID | 1.016 | 1.022 | 0.423 | 0 | -0.071 | 0.183 | 0.659 | 94.2 | 0.375 | 0.412 | 0.745 | 48.8
$\gamma_{t}-\gamma_{t-1}=c_{t}$: $\beta_{t}-\beta_{t-1}=d_{t}$
AugSynth | 0.529 | 0.538 | 2.114 | 100 | 0.57 | 0.578 | 2.279 | 100 | 0.421 | 0.439 | 1.686 | 100
DID | 0.793 | 0.807 | 0.574 | 0.2 | 0.475 | 0.492 | 0.492 | 3.4 | 0.078 | 0.174 | 0.629 | 92.2
Ign-1 | 1.168 | 1.173 | 0.468 | 0 | 0.849 | 0.855 | 0.434 | 0 | 0.457 | 0.47 | 0.428 | 0.8
Ign-2 | 0.54 | 0.55 | 0.395 | 0.2 | 0.535 | 0.543 | 0.391 | 0 | 0.21 | 0.242 | 0.398 | 39.2
gDID-1 | 0.629 | 0.64 | 0.504 | 0 | 0.262 | 0.296 | 0.537 | 52.4 | -0.061 | 0.173 | 0.63 | 95
gDID-2 | 0.001 | 0.136 | 0.511 | 93.2 | -0.03 | 0.131 | 0.522 | 95.4 | -0.271 | 0.327 | 0.644 | 62.8
gDiD-0 | 0.007 | 0.13 | 0.502 | 93 | 0.149 | 0.184 | 0.391 | 65.4 | -0.212 | 0.259 | 0.515 | 61.6
sDID | -0.557 | 0.585 | 0.694 | 11.6 | -0.846 | 0.864 | 0.738 | 0.4 | -0.702 | 0.733 | 0.872 | 12