# Efficient estimation of longitudinal treatment effects using difference-in-
differences and machine learning

Nicholas Illenberger, Iván Díaz, Audrey Renson
Department of Population Health, New York University Grossman School of
Medicine, New York, United States
Contact<EMAIL_ADDRESS>
###### Abstract

Difference-in-differences is based on a parallel trends assumption, which
states that changes over time in average potential outcomes are independent of
treatment assignment, possibly conditional on covariates. With time-varying
treatments, parallel trends assumptions can identify many types of parameters,
but most work has focused on group-time average treatment effects and similar
parameters conditional on the treatment trajectory. This paper focuses instead
on identification and estimation of the intervention-specific mean – the mean
potential outcome had everyone been exposed to a proposed intervention – which
may be directly policy-relevant in some settings. Previous estimators for this
parameter under parallel trends have relied on correctly-specified parametric
models, which may be difficult to guarantee in applications. We develop
multiply-robust and efficient estimators of the intervention-specific mean
based on the efficient influence function, and derive conditions under which
data-adaptive machine learning methods can be used to relax modelling
assumptions. Our approach allows the parallel trends assumption to be
conditional on the history of time-varying covariates, thus allowing for
adjustment for time-varying covariates possibly impacted by prior treatments.
Simulation results support the use of the proposed methods at modest sample
sizes. As an example, we estimate the effect of a hypothetical federal minimum
wage increase on self-rated health in the US.

##  1 Introduction

Difference-in-differences (DID) methods are popular in policy analysis because
they can identify effects in the presence of unmeasured confounding, which is
often expected in comparisons across jurisdictions. DID relies on parallel
trends assumptions, which state that time trends in average potential outcomes
are independent of treatment assignment (possibly conditional on covariates).
This is in contrast to much of the causal inference literature, where
identification is based on exchangeability assumptions; these state that
levels (rather than trends) of potential outcomes are independent of treatment
assignment (possibly conditional on covariates). Parallel trends can (but does
not necessarily) hold in the presence of unmeasured common causes of the
outcome and treatment (i.e., unmeasured confounding) whereas exchangeability
cannot. For example, when treatment only changes once (i.e. is time-
invariant), parallel trends holds if and only if there is a constant
unmeasured confounding effect pre- and post-treatment [29]. The simluation
data-generating mechanism in Section 3 of this paper provides an example of
how parallel trends may hold under unmeasured confounding in a more
complicated time-varying scenario.

Until recently, most methodological analyses of DID have considered only point
treatments [1, 15], but realistic settings often have time-varying treatments,
making identification and estimation more complex. Time-varying treatments are
usually accomodated in the applied literature via two-way fixed effects (TWFE)
models, which rely on strong parametric assumptions and/or may make illogical
comparisons, resulting in potentially large bias [14]. More recently, methods
have been developed that relax these assumptions by considering DID for time-
varying treatments in a non-parametric or semi-parametric framework [5, 8, 22,
26]. These methods are increasingly used in practice and can produce markedly
different results than TWFE [13].

In DID with time-invariant binary treatments, there is typically only one
parameter of interest, the average treatment effect in the treated (ATT).
However, with time-varying and/or many-valued treatments, there is not
necessarily an obvious generalization of the ATT since there may be many
observed treatment trajectories; therefore, many different causal parameters
can potentially be of interest. The recent DID literature has focused
primarily on parameters built on the ATT concept, such as the group-time
average treatment effects [5] or cohort average treatment effects [30]; like
the ATT, these parameters condition on an observed treatment trajectory and
apply mainly to binary treatments. A different parameter which may also be of
interest and can be identified in time-varying DID settings is the
intervention-specific mean, or the the average potential outcome for the whole
population under a specified treatment plan [22]. Intervention-specific means
have been the major focus in the literature on longitudinal treatment effects
under sequential exchangeability assumptions, which propose that potential
outcomes are independent of a time-varying treatment given the history of
covariates and treatments. Intervention-specific means have received less
attention in the DID literature, but may be more policy-relevant in some
settings. For example, [22] used a DID framework to estimate counterfactual US
mortality rates had all US states remained longer under pandemic stay-at-home
orders, such a parameter may be directly applicable to pandemic planning
decisions at the federal level.

One challenge with time-varying treatments is that there are likely to be
time-varying covariates, which may be necessary to condition on in order for
the parallel trends assumption to hold. Time-varying covariates in the
parallel trends assumption have received limited attention in the DID
literature [8, 3, 22, 26]. When the parallel trends assumption depends on
time-varying covariates, the standard approach in applied literature has been
TWFE models with covariates entered as regressors. When analyzing these
models, it is generally assumed (in addition to the other strong parametric
assumptions) that covariates are not affected by past treatment [34, 14,
e.g.]; such covariates are termed “bad controls” in this literature [2].
Analyses of DID models that relax the TWFE assumptions have generally also
ruled out feedback between treatment and covariates [3, 12]. However, such
treatment-covariate feedback is likely to be common in many DID applications;
[4] and [22] both provide examples.

Previously developed estimators in time-varying DID settings have generally
required that at least one of either the propensity score or outcome
regression are correctly-specified parametric models [5, 22]. Such models may
be high-dimensional and difficult to specify correctly in practice, especially
when estimation depends on the history of time-varying covariates. Estimators
that allow for appropriate inference when data-adaptive models are used to
estimate nuisance parameters have been developed primarily for settings where
identification is based on exchangeability [10, e.g.]. Such an approach has
been extended to estimating the ATT in a time-invariant DID setting [6]. The
extension of this approach to a time-varying DID setting has, to our
knowledge, only been considered in the context of g-estimation of structural
nested models [26]. Though this approach allows estimation of the
intervention-specific mean and other marginal parameters, the procedure is
quite complex and its main strength is in studying time-varying effect
heterogeneity.

Building on the identification results of [22], we develop estimators for the
intervention-specific mean under parallel trends based on a nonparametric
statistical model. Our estimators are based on the efficient influence
function for the identified parameter, and thus are nonparametrically
asymptotically efficient, multiply robust, and can accomodate data-adaptive
machine learning to estimate the propensity scores and/or outcome regression
functions. Our approach depends on a parallel trends assumption that
conditions on the history of time-varying covariates up to each treatment time
point, and does not require an assumption that covariates are unaffected by
prior treatment. Though our focus is on longitudinal settings with repeated
observations among units, our results could be extended to repeated cross-
sectional data as well. The remainder of this paper is organized as follows.
Section 2 presents identification results and estimators based on the
efficient influence function, Section 3 presents a simulation study exploring
finite-sample properties of the proposed estimators, Section 4 presents a
worked example studying the effect of minimum wages on self-rated health, and
Section 5 concludes.

##  2 Methods

###  2.1 Data

Suppose we collect data on $i=1,...,n$ participants at $t=0,...,\tau$ discrete
timepoints. For participant $i$ at timepoint $t$, we observe a discrete
exposure variable $A_{it}$, a set of time-varying covariates $X_{it}$, and the
outcome of interest $Y_{it}$. The full collection of observed data takes the
form of $n$ i.i.d. copies of
$O_{i}=(X_{i0},A_{i0},Y_{i0},X_{i1},...,X_{i\tau},A_{i\tau},Y_{i\tau})$. Let
$\overline{A}_{t}=(A_{0},A_{1},...,A_{t})$ denote the treatment history up to
and including time $t$, and let $\overline{A}\equiv\overline{A}_{\tau}$ denote
the entire history. We define $\overline{X}_{t}$ and $\overline{Y}_{t}$
similarly. For convenience, we define $Z_{t}=\\{\emptyset\\}$ when $t<0$ for
any random variable $Z$. We assume that data come from a staggered
discontinuation design [22] such that $A_{0}=a_{0}^{*}$ for all units, where
$a_{0}^{*}$ is the value that $A_{0}$ takes under the counterfactual treatment
trajectory to which the target causal effect refers (defined in the next
section). In other words, we assume all units begin follow-up with treatment
value consistent with the counterfactual regime of interest, but may later
discontinue. Note that this design includes as a special case staggered
adoption designs, often the focus of time-varying DID studies [5], where
$A_{t}\in\\{0,1\\}$, $A_{0}=0$ for everyone, and $A_{t}=1$ implies $A_{t+1}=1$
with probability 1, but also includes settings where treatments switch on and
off and/or are non-binary.

###  2.2 Target causal effects

We will work within the potential outcomes framework, wherein
$Y_{t}(\overline{a}_{t})$ represents the value of $Y_{t}$ that would be
observed if an individual’s exposure history had been set by intervention to
$\overline{a}_{t}$. Note that this potential outcome is only observed for (at
most) those participants whose observed treatment pattern follows
$\overline{A}_{t}=\overline{a}_{t}$.

We will focus on causal quantities of the form
$\mu_{t}=\mathbb{E}[Y_{t}(\overline{a}^{*})]$, where $\overline{a}^{*}$
represents a counterfactual trajectory of interest. Thus, the quantity
$\mu_{t}$ is an “intervention-specific mean”; i.e., an average potential
outcome, had everyone’s exposure trajectory $\overline{A}$ been set by
intervention to $\overline{a}^{*}$. Typically, counterfactual outcomes are
compared to observed outcomes in some way, so that the target causal quantity
is actually $\mathbb{E}[Y_{t}-Y_{t}(\overline{a}^{*})]$. However, for such
comparisons, estimation and inference are typically straightforward once these
are achieved for the counterfactual component $\mu_{t}$; hence, the focus of
this article is the latter. Throughout, we assume there is only one
counterfactual trajectory $\overline{a}^{*}$ of interest; the methods
presented in this article also apply for multiple trajectories, but assuming
parallel trends for multiple trajectories places restrictions on treatment
effect homogeneity which may not be desirable [22, callaway2021continuous].

###  2.3 Identification

Because potential outcomes under a given regime are unobserved, certain
assumptions are required to connect the distribution of these potential
outcomes to that of the observable data. In our context, identification of the
causal effect of $A$ on $Y$ will rely on a generalization of the parallel
trends assumption to a longitudinal setting. Specifically, we require the
following identification assumptions:

######  Assumption 1 (No interference).

$Y_{it}(\overline{a}_{i}^{*},\overline{a}_{i^{\prime}}^{*})=Y_{it}(\overline{a}_{i})^{*}$
for $i\neq i^{\prime}$

######  Assumption 2 (Treatment version irrelevance).

If $\overline{A}_{it}=\overline{a}^{*}_{t}$, then
$Y_{it}=Y_{it}(\overline{a}^{*}_{t})$

######  Assumption 3 (No anticipation).

$Y_{it}(\overline{a}^{*}_{s})=Y_{it}(\overline{a}_{t}^{*})$ for $s>t$

######  Assumption 4 (Parallel trends).

For $t\in\\{1,2,...,\tau\\}$ and $k\in\\{0,1,...,t-1\\}$, assume there is a
$\overline{W}_{k}\subseteq(\overline{X}_{k},\overline{Y}_{k-2})$ such that

$E[Y_{t}(\overline{a}^{*})-Y_{t-1}(\overline{a}^{*})|\overline{W}_{k},\overline{A}_{k-1}=\overline{a}_{k-1}^{*}]=E[Y_{t}(\overline{a}^{*})-Y_{t-1}(\overline{a}^{*})|\overline{W}_{k},\overline{A}_{k}=\overline{a}_{k}^{*}]$

where $\overline{W}_{k-1}\subseteq\overline{W}_{k}$.

######  Assumption 5 (Positivity).

If $f(w_{t}|\overline{w}_{t-1},\overline{a}_{t-1}^{*})>0$, then
$f(\overline{a}_{t}^{*}|\overline{w}_{t},\overline{a}_{t-1}^{*})>0$ for
$\overline{w}_{t}\in\overline{\mathcal{W}}_{t}$, $t\in\\{0,1,...,\tau\\}$

Assumptions 1-3 are standard in the longitudinal DID literature; we briefly
summarize these here. No interference requires that outcomes of one unit are
not affected by the treatment a separate unit receives (i.e. the potential
outcome only depends on that unit’s observed treatment). Treatment version
irrelevance requires that exposures are defined and measured with sufficient
granularity that, if a given unit is observed under the treatment pattern
$\overline{a}^{*}$, then their observed outcome is equal to their potential
outcome under an intervention to set $\overline{A}=\overline{a}^{*}$.
Assumptions 1 and 2 are sometimes combined as the Stable Unit Treatment Value
Assumption. No anticipation requires that treatments observed at a time $s>t$
after outcomes $Y_{it}$ did not have the ability to impact that prior outcome;
for example, this may be violated if units altered their behavior in
anticipation of planned policy changes.

The parallel trends assumption (Assumption 4) requires that there is some
subset of the covariate history (up to time $k$) and the outcome history (up
to time $k-2$), conditional on which the average changes in potential outcomes
on the additive scale over successive periods are equal across values of the
current treatment. (Note that if $Y_{k-1}$ were included, then parallel trends
would reduce to a sequential exchangeability assumption and DID methods would
not be necessary.) Assumption 4 is identical to the parallel trends assumption
previously considered in the context of the intervention-specific mean [22].
It differs from parallel trends assumptions considered in other popular
longitudinal DID methods in that it conditions on the history of time-varying
covariates; for example, the methods of [5] condition only on baseline
covariates, those of [8] condition on contemporaneous covariates, and those of
[3] condition on the entire covariate history (i.e. including covariates
measured later than the outcome). Importantly, we require no assumption that
covariates are unaffected by treatment; thus, our methods can account for
treatment-covariate feedback. Finally, positivity requires that there be
positive probability of receiving the proposed treatment pattern $a_{t}^{*}$,
for individuals with treatment history $\overline{a}_{t-1}^{*}$ and at all
levels of the covariate/outcome history needed for parallel trends.

Following [22], under Assumptions 1-4, we have:

$\displaystyle\mu_{t}$
$\displaystyle=\mathbb{E}[Y_{0}(\overline{a}^{*})]+\sum_{k=1}^{t}\mathbb{E}[Y_{k}(\overline{a}^{*})-Y_{k-1}(\overline{a}^{*})]$
$\displaystyle=\mathbb{E}[Y_{0}(\overline{a}^{*})]+\sum_{k=1}^{t}\mathbb{E}(\mathbb{E}[Y_{k}(\overline{a}^{*})-Y_{k-1}(\overline{a}^{*})|\overline{W}_{1},A_{1}])$
$\displaystyle=\mathbb{E}[Y_{0}(\overline{a}^{*})]+\sum_{k=1}^{t}\mathbb{E}(\mathbb{E}[Y_{k}(\overline{a}^{*})-Y_{k-1}(\overline{a}^{*})|\overline{W}_{1},A_{1}=a_{1}^{*}])$
$\displaystyle\vdots$
$\displaystyle=\mathbb{E}[Y_{0}(\overline{a}^{*})]+\sum_{k=1}^{t}\mathbb{E}(\mathbb{E}[\cdots\mathbb{E}\\{Y_{k}(\overline{a}^{*})-Y_{k-1}(\overline{a}^{*})|\overline{W}_{k},\overline{A}_{k}\\}\cdots|\overline{W}_{1},A_{1}=a_{1}^{*}])$
$\displaystyle=\mathbb{E}[Y_{0}(\overline{a}^{*})]+\sum_{k=1}^{t}\mathbb{E}(\mathbb{E}[\cdots\mathbb{E}\\{Y_{k}(\overline{a}^{*})-Y_{k-1}(\overline{a}^{*})|\overline{W}_{k},\overline{A}_{k}=\overline{a}_{k}^{*}\\}\cdots|\overline{W}_{1},A_{1}=a_{1}^{*}])$
$\displaystyle=\mathbb{E}[Y_{0}(a_{0}^{*})]+\sum_{k=1}^{t}\mathbb{E}(\mathbb{E}[\cdots\mathbb{E}\\{Y_{k}(\overline{a}^{*}_{k})-Y_{k-1}(\overline{a}^{*}_{k})|\overline{W}_{k},\overline{A}_{k}=\overline{a}_{k}^{*}\\}\cdots|\overline{W}_{1},A_{1}=a_{1}^{*}])$
$\displaystyle=\mathbb{E}[Y_{0}]+\sum_{k=1}^{t}\mathbb{E}(\mathbb{E}[\cdots\mathbb{E}\\{Y_{k}-Y_{k-1}|\overline{W}_{k},\overline{A}_{k}=\overline{a}_{k}^{*}\\}\cdots|\overline{W}_{1},A_{1}=a_{1}^{*}])$

The first equality of the above proof holds by telescoping, the second and
fourth by law of iterated expectation, the third and fifth by parallel trends,
the sixth by no anticipation, and the seventh by no interference and treatment
version irrelevance because
$Y_{k}(\overline{a})-Y_{k-1}(\overline{a})=Y_{k}-Y_{k-1}$ in the event
$\overline{A}_{k}=\overline{a}_{k}$. Note that in the last line above,
covariates $\overline{W}_{k}$ and treatments $\overline{A}_{k}$ may be
measured later than outcomes $Y_{k-1}$; thus, in estimators we consider below,
outcomes may be regressed on future variables, which may appear counter-
intuitive. This concern may be addressed by noting that differences
$Y_{k}-Y_{k-1}$ are only determined in the future of $\overline{W}_{k}$ and
$\overline{A}_{k}$.

In order to simplify notation in our definitions of estimators, we
additionally define:

$\displaystyle\mathbb{E}[Y_{0}]$
$\displaystyle+\sum_{k=1}^{t}\mathbb{E}(\mathbb{E}[\cdots\mathbb{E}\\{Y_{k}-Y_{k-1}|\overline{W}_{k},\overline{A}_{k}=\overline{a}_{k}^{*}\\}\cdots|\overline{W}_{1},A_{1}=a_{1}^{*}])$
$\displaystyle=\mathbb{E}[Y_{0}]+\sum_{k=1}^{t}(\phi_{k,k}-\phi_{k-1,k})$
$\displaystyle=\psi_{t}$

with $\psi_{t}=\mathbb{E}[Y_{0}]+\sum_{k=1}^{t}(\phi_{k,k}-\phi_{k-1,k})$ and
$\phi_{j,k}=\mathbb{E}\\{\mathbb{E}[\cdots\mathbb{E}(Y_{j}\mid\overline{W}_{k},\overline{A}_{k}=\overline{a}_{k}^{*})\cdots|\overline{W}_{1},A_{1}=a_{1}^{*}]\\}$.
Note that $\phi_{j,k}$ is the time-varying g-formula [23], which identifies
causal effects of the type $\mu_{t}$ under sequential exchangeability. Here,
$\phi_{j,k}$ has no causal interpretation in itself, but instead forms the
building blocks of the parameter $\psi_{t}$ (which equals $\mu_{t}$ under our
assumptions). However, established methods to estimate $\phi_{j,k}$ based on
its efficient influence function form the basis of our proposed estimators.

###  2.4 Estimators

For the rest of this paper we assume that Assumptions 1-5 hold, so that the
intervention-specific mean is identified as $\mu_{t}=\psi_{t}$, and focus
attention on estimating the quantity $\psi_{t}.$ Before defining the efficient
influence function and resulting estimator, we will introduce some additional
notation. Let
$g_{m}(\overline{W}_{m})=\prod_{k=1}^{m}f(A_{k}=a_{k}^{*}|\overline{W}_{k},\overline{A}_{k-1}=\overline{a}_{k-1}^{*})$
denote the cumulative propensity score up to time $m$. Define
$Q^{j,k,k+1}\equiv Y_{j}$, and let $Q^{j,k,m}$ be iteratively defined as

$Q^{j,k,m}(\overline{W}_{m})=E[Q^{j,k,m+1}(\overline{W}_{m+1})|\overline{A}_{m}=\overline{a}_{m}^{*},\overline{W}_{m}].$

It can be seen that $\phi_{j,k}=E[Q^{j,k,1}(\overline{W}_{1})].$ In what
follows, we use the notation $P_{n}X$ to denote the empirical average of $X$,
$n^{-1}\sum_{i=1}^{n}X_{i}$. Whenever there is no ambiguity, we drop the
argument $\overline{W}_{m}$ from $g_{m}(\overline{W}_{m})$ and
$Q^{j,k,m}(\overline{W}_{m})$ for conciseness. For example, we let
$\widehat{Q}_{i}^{j,k,m}$ denote an estimate of $Q^{j,k,m}(\overline{W}_{m})$
for observation $i$. For a function $f$ of $\overline{W}_{m}$, we let
$||f||^{2}=\int[f(\overline{w}_{m})]^{2}dP(\overline{w}_{m})$.

The following theorem defines the efficient influence function (EIF) for the
parameter $\psi_{t}$, upon which multiply robust and efficient data-adaptive
estimators will be based.

######  Theorem 1.

The Efficient Influence Function for $\psi_{t}$ is given by

$\displaystyle IF(\psi_{t})$
$\displaystyle=IF(\mathbb{E}[Y_{0}])+\sum_{k=1}^{t}\big{\\{}IF(\phi_{k,k})-IF(\phi_{k-1,k})\big{\\}}$

where $IF(\mathbb{E}[Y_{0}])=Y_{0}-\mathbb{E}[Y_{0}]$ and
$IF(\phi_{j,k})=\sum_{m=1}^{k}\frac{I(\overline{A}_{m}=\overline{a}_{m}^{*})}{g_{m}(\overline{W}_{m})}(Q^{j,k,m+1}-Q^{j,k,m})+Q^{j,k,1}-\phi_{j,k}.$

Given the EIF a “one-step” estimator for $\psi_{t}$ is given by

$\widehat{\psi}_{t}=P_{n}\Bigg{\\{}Y_{0}+\sum_{k=1}^{t}\bigg{(}\widetilde{\phi}_{k,k}(\widehat{\Theta})-\widetilde{\phi}_{k-1,k}(\widehat{\Theta})\bigg{)}\Bigg{\\}}$
(1)

where:

$\widetilde{\phi}_{jki}(\widehat{\Theta})=\sum_{m=1}^{k}\frac{I(\overline{A}_{mi}=\overline{a}^{*}_{m})}{\widehat{g}_{mi}}(\widehat{Q}_{i}^{j,k,m+1}-\widehat{Q}_{i}^{j,k,m})+\widehat{Q}_{i}^{j,k,1}$

where $\widehat{Q}_{i}^{j,k,m}$ and
$\widehat{g}_{mi}=\widehat{g}_{m}(\overline{W}_{mi})$ denote estimators of the
referenced quantities, and
$\widehat{\Theta}=\\{\widehat{Q}_{i}^{j,k,m},\widehat{g}_{m}:m=1,...,k;k=1,...,t\\}$
denotes the full collection of nuisance function estimators. In the next
section we will describe asymptotic behavior of this estimator, including
multiple robustness, efficiency, and convergence rates under difference
choices of nuisance function estimators.

###  2.5 Asymptotic Results

This section describes the asymptotic distribution of the one-step estimator
$\widehat{\psi}_{t}$ of the intervention-specific mean. This discussion
includes (1) inferential techniques based on the asymptotic normality of the
estimator, (2) conditions under which the estimator is unbiased, and (3)
description of a sample splitting procedure that allows for efficient
estimation under flexible data-adaptive estimators of the nuisance functions.

######  Theorem 2.

Assume the following conditions:

  * (C1)

The estimated nuisance functions, $\hat{\Theta}$, are in a Donsker class, and
are such that
$||IF[\phi_{j,k}(\hat{\Theta})]-IF[\phi_{j,k}]||\xrightarrow{p}0$, and

  * (C2)

For some constant $\epsilon>0$,
$P[\epsilon<\widehat{g}_{k}(\overline{W}_{m})]=P[\epsilon<g_{k}(\overline{W}_{m})]=1$,
for all $k$ and $\overline{a}^{*}$,

where $\xrightarrow{p}$ denotes convergence in probability. Let

$R_{2,n}=\sum_{m=1}^{t}\left\lVert
g_{m}-\widehat{g}_{m}\right\rVert\sum_{k=m}^{t}\bigg{(}\left\lVert
Q^{k,k,m}-\widehat{Q}^{k,k,m}\right\rVert+\left\lVert
Q^{k-1,k,m}-\widehat{Q}^{k-1,k,m}\right\rVert\bigg{)}$

Under (C1) and (C2), we have:

$\widehat{\psi}_{t}-\psi_{t}=(P_{n}-P)IF(\psi_{t})+o_{p}(n^{-1/2})+O_{p}(R_{2,n})$

###### Proof.

By Lemma 1 of [19] we have

$\widetilde{\phi}_{jk}-\phi_{jk}=(P_{n}-P)IF(\phi_{jk})+(P_{n}-P)(\widehat{IF}(\phi_{jk})-IF(\phi_{jk}))+Rem_{jk},$

where $\widehat{IF}(\phi_{jk})$ denotes the IF of Theorem 1 with $\Theta$
replaced by $\widehat{\Theta}$, and where

$Rem_{jk}=\sum_{m=1}^{k}\int
c_{m}(\overline{w}_{m})\\{g_{m}(\overline{w}_{m})-\widehat{g}_{m}(\overline{w}_{m})\\}\\{Q^{j,k,m}(\overline{w}_{m})-\widehat{Q}^{j,k,m}(\overline{w}_{m})\\}dP(\overline{w}_{m}),$

where $c_{m}(\overline{W}_{m})$ is a constant given in the original Lemma,
which is bounded in probability according to C2. Under C1, Lemma 19.24 of [32]
shows that the term $(P_{n}-P)(\widehat{IF}(\phi_{jk})-IF(\phi_{jk}))$ is
$o_{P}(n^{-1/2})$. Likewise, using the Cauchy-Schwarz inequality we have

$|Rem_{jk}|\leq\left\lVert g_{m}-\widehat{g}_{m}\right\rVert\left\lVert
Q^{j,k,m}-\widehat{Q}^{j,k,m}\right\rVert.$

Applying the above to

$\psi_{t}=\mathbb{E}[Y_{0}]+\sum_{k=1}^{t}(\phi_{k,k}-\phi_{k-1,k});\quad\widetilde{\psi}_{t}=P_{n}Y_{0}+\sum_{k=1}^{t}(\widetilde{\phi}_{k,k}-\widetilde{\phi}_{k-1,k})$

together with Theorem 1 yields the result. ∎

Theorem 2 describes the asymptotic behavior of our proposed estimator. The
first error component illustrates the asymptotic bias of our estimator, as
well as the importance of using doubly robust estimating equations such as the
EIF when using machine learning. First, note that $R_{2,n}$ is defined
essentially by sums of product of errors in the estimation of the nuisance
parameters $g_{m}$ and $Q^{k,k,m}$ as well as $Q^{k-1,k,m}$. Machine learning
regression algorithms are specifically designed to make each of these errors
small. Therefore, the product of each pair of errors is expected to be smaller
than each of its components. In comparison, if these nuisance regressions are
estimated with parametric models, large regression errors from model
misspecificaiton are more likely. Therefore, one of these errors could be
large, which would make the product larger than each of its components [24],
leading to large bias.

Estimators constructed using the efficient influence function as above are
often called doubly-robust because they yield consistent estimators even in
the case of misspecification of one of the models under consideration.
According to the above theorem, the estimator $\widehat{\psi}_{t}$ is
consistent if $R_{2.n}=o_{p}(1)$. This holds if, for each time point $m\leq
t$, we have either $g_{m}$ is consistently estimated, or the sequence
$S_{m}=\\{(Q^{k,k,m},Q^{k-1,k,m}):m\leq k\leq t\\}$ is consistenly estimated,
where we note that consistent estimation of $S_{m}$ requires consistent
estimation of $S_{m+1}$. Assume that we start with a consistent outcome
regression estimator at time point $t$, but then at time $m\leq t$ we have an
inconsistent outcome regression. Then we would require that the propensity
score is estimated consistently at all times prior to and including $m$, since
in this case $S_{j}:j\leq m$ cannot be consistently estimated anymore. Since
we have $t+1$ possible values for $m$, this means that we have $t+1$ chances
of consistency. Importantly, this also implies that using an adaptation of the
sequential doubly robust estimators of [19] and [25] would not yield
estimators that are more robust than a simple one-step estimator.

The result of the theorem also allows us to study the weak convergence of the
estimator. In particular, if $R_{2,n}=o_{p}(n^{-1/2})$, then it will be the
case that

$\sqrt{n}(\widehat{\psi}_{t}-\psi_{t})\xrightarrow{d}N[0,\left\lVert
IF(\psi_{t})\right\rVert^{2}]$

where $\xrightarrow{d}$ denotes convergence in distribution. In this case, an
estimator of the asymptotic variance of $\widehat{\psi}_{t}$ is given by:

$\widehat{V}(\widehat{\psi}_{t})=\frac{1}{n}\sum_{i=1}^{n}\left\\{Y_{0i}+\sum_{k=1}^{t}\bigg{(}\widetilde{\phi}_{k,k,i}(\widehat{\Theta})-\widetilde{\phi}_{k-1,k,i}(\widehat{\Theta})\bigg{)}-\widehat{\psi}_{t}\right\\}^{2}.$

A sufficient condition for the second order term of the theorem to be
$o_{p}(n^{-1/2})$ is that all nuisance parameters are estimated consistently
at rates $n^{-1/4}$. These rates are achievable by some machine learning
algorithms, under assumptions.

At this point it is important to clarify a semantic note on the use of the
term “double robustness”. In this paper, as in most manuscripts that use
machine learning for estimation of causal effects, we use the term double
robustness to refer to the property that the estimator is consistent under
inconsistent estimation of some of the nuisance estimators, as stated above.
Some of the literature on causal inference, especially the literature focused
on parametric models, uses the term double robustness to refer to the property
that the estimators are $\sqrt{n}$-consistent and asymptotically normal under
misspecification of some of the nuisance models. While our proposed estimator
would also enjoy asymptotic normality under misspecification of some models if
$\hat{\Theta}$ is estimated within pre-specified parametric models, in this
paper we are focused on addressing model misspecification bias and therefore
focus on results for machine learning nuisance estimation, using the term
double-robustness to refer to consistency of the estimator under model
misspecification.

###  2.6 Sample-Split Estimator

The asymptotic results for the proposed estimator relies on the condition that
the estimated nuisance functions, $\hat{\Theta}$, are within the Donsker
class. This condition controls how close functions of the estimated nuisance
functions can get to their limiting versions (e.g. the difference between
$(P_{n}-P)\tilde{\phi}_{jk}(\widehat{\Theta})$ and
$(P_{n}-P)\tilde{\phi}_{jk}(\Theta)$) [16]. This condition can be restrictive
if estimators with unbounded variation norm are used, as these are generally
not Donsker [9, e.g.,]. This condition can be removed by implementing sample-
split versions of the estimators. Partition the index set $\\{1,...,n\\}$ into
$M$ mutually exclusive groups $V_{1},...,V_{M}$ such that
$\bigcup_{m=1}^{M}V_{m}=\\{1,...,n\\}$ and $V_{m}\bigcap
V_{m^{{}^{\prime}}}=\emptyset$ for $m\neq m^{{}^{\prime}}$. For each
validation set, $m$, let
$\mathcal{T}_{m}=\\{1,...,n\\}\setminus\mathcal{V}_{m}$ denote the
corresponding training set. Let $m(i)$ denote the index of the validation set
containing unit $i$ (e.g. $m(i)=s$ if $i\in\mathcal{V}_{s}$). Finally, let
$\widehat{\Theta}_{m}=(\widehat{Q}_{m},\widehat{g}_{m})$ denote the estimated
nuisance functions obtained using training data $\mathcal{T}_{m}$.

The sample-split estimator of the desired causal estimand is obtained by
replacing the full-sample estimated nuisance functions in Equation 1 with
their training-sample estimated counterparts for each unit.

$\widehat{\psi}_{t,ss}=P_{n}\Bigg{\\{}Y_{0i}+\sum_{k=1}^{t}\bigg{(}\widetilde{\phi}_{k,k}(\widehat{\Theta}_{m(i)})-\widetilde{\phi}_{k-1,k}(\widehat{\Theta}_{m(i)})\bigg{)}\Bigg{\\}}$

Note that within each validation set, $V_{m}$, the estimated nuisance
functions $\widehat{\Theta}$ are fixed. This fact allows us to apply the
empirical process properties that our asymptotic results are based upon
without invoking Donsker class conditions.

Because estimated treatment effects depend on the random partitioning of the
data, we may repeat the procedure over different partitions and combine
estimates of the treatment effect and variance. Let $\widehat{\psi}_{t,ss,k}$
for $k=1,...,K$ and $\widehat{V}(\widehat{\psi}_{t,ss,k})$ denote the
estimated treatment effect and variance for the $k$-th repetition. The
combined treatment effect can be estimated as the median estimate across the
$K$ partitions of the data,
$\widetilde{\psi}_{t,ss}=\text{median}_{k}(\widehat{\psi}_{t,ss,k})$ [35]. As
opposed to the mean, the median is more robust to outliers across the
partitioned estimates [7]. The estimated variance of this estimator is given
by
$\widehat{V}(\widetilde{\psi}_{t,ss})=\text{median}_{k}(\widehat{V}(\widehat{\psi}_{t,ss,k})+(\widehat{\psi}_{t,ss,k}-\widetilde{\psi}_{t,ss})^{2})$.
This estimator incorporates the uncertainty arising from the random
partitioning of the data.

##  3 Simulation study

We generated 500 datasets each at sample sizes $n=1,000$, $5,000$, and
$20,000$ according to the following structural model for $t=0,1,2$:

$\displaystyle U$ $\displaystyle\sim N(0,1)$ $\displaystyle W_{1t}$
$\displaystyle\sim
N(\gamma_{01t}+\gamma_{11t}W_{1,t-1}+\gamma_{21t}A_{t-1},0.1)$ $\displaystyle
W_{2t}$ $\displaystyle\sim
N(\gamma_{02t}+\gamma_{12t}W_{1t}+\gamma_{22t}W_{2,t-1}+\gamma_{32t}A_{t-1},0.1)$
$\displaystyle W_{3t}$ $\displaystyle\sim
N(\gamma_{03t}+\gamma_{13t}W_{1t}+\gamma_{23t}W_{2t}+\gamma_{33t}W_{3,t-1}+\gamma_{43t}A_{t-1},0.1)$
$\displaystyle A_{t}$
$\displaystyle\sim\text{Bernoulli}\\{(1-A_{t-1})\text{logit}(\alpha_{0t}+\alpha_{1t}W_{1t}+\alpha_{2t}\cos
W_{2t}+\alpha_{3t}W_{3t}^{2}+\alpha_{4t}U)\\}$ $\displaystyle Y_{t}$
$\displaystyle\sim N(\beta_{0t}+\beta_{1t}\sin
W_{1t}+\beta_{2t}W_{2t}W_{3t}+\beta_{3t}A_{t}+U,0.1)$

All true parameter values were generated once from a standard normal
distribution. To see that parallel trends holds in this simulation, note that
if $A_{t}$ is set to zero for all $t$, differences $Y_{t}-Y_{t-1}$ will not
depend on $U$ because $U$ enters the model for $Y_{t}$ additively and does not
enter the model for $W_{t}$. In each simulated dataset, we estimated the
parameter $\mu_{t}=\mathbb{E}[Y_{t}(\overline{A}=\overline{0})]$, the
intervention-specific mean under the policy to set everyone to be untreated.
Nuissance functions were estimated using parametric models with treatment and
outcome correctly specified (true), and misspecifying the treatment (gfal),
outcome (qfal), or both (bfal); all without sample splitting. For each set of
parametric nuissance functions, the estimator $\widehat{\psi}_{t}$ and
variance estimator $\widehat{V}(\widehat{\psi}_{t})$ were calculated. Finally,
nuissance functions were estimated using Super Learner with a library
including the mean, lasso, multivariate adaptive regression splines, and
kernal support vector machines, using sample splitting with $M=10$ random
partitions. The estimator $\widehat{\psi}_{t,ss}$ and variance estimator
$\widehat{V}(\widehat{\psi}_{t,ss})$ were then calculated.

Results from the simulation study are shown in Table 1 and Figure 1. At all
sample sizes considered (ranging from $n=1,000$ to $50,000$) and whenever
either set of models is correctly specified, estimators appear to be
converging at root-n rate to a normal distribution centered at zero and with
variance close to the EIF-based variance estimate. This is true regardless of
whether data-adaptive or parametric nuisance estimators are used, suggesting
there is little loss of efficiency in using data-adaptive estimators. Though
the consistency of the proposed variance estimator depends theoretically on
all models being correctly specified, in this particular example the results
do not appear very sensitive to this requirement. That is, the mean of
estimated variances appears very close to the empirical variance of estimates
over repeated simulations, even when only either the propensity score or
outcome models are correct.

$n_{obs}$ | method | bias $(\widehat{\psi}_{1})^{2}$ | $\widehat{V}_{sim}(\widehat{\psi}_{1})$ | $\widehat{V}_{eif}(\widehat{\psi}_{1})$ | bias $(\widehat{\psi}_{2})^{2}$ | $\widehat{V}_{sim}(\widehat{\psi}_{2})$ | $\widehat{V}_{eif}(\widehat{\psi}_{2})$
---|---|---|---|---|---|---|---
1000 | bfal | 0.1269 | 0.1952 | 0.2024 | 0.2119 | 0.3131 | 0.3790
| gfal |  $<$0.0001 | 0.1889 | 0.1882 | 0.0001 | 0.3122 | 0.3229
| qfal |  $<$0.0001 | 0.1915 | 0.2138 |  $<$0.0001 | 0.3233 | 0.4417
| super | 0.0490 | 0.2224 | 0.1670 | 0.0029 | 0.3424 | 0.2847
| true |  $<$0.0001 | 0.1911 | 0.1918 |  $<$0.0001 | 0.3181 | 0.3333
5000 | bfal | 0.1280 | 0.0349 | 0.0404 | 0.2219 | 0.0676 | 0.0759
| gfal |  $<$0.0001 | 0.0349 | 0.0376 | 0.0004 | 0.0629 | 0.0648
| qfal |  $<$0.0001 | 0.0350 | 0.0427 | 0.0002 | 0.0682 | 0.0881
| super | 0.0017 | 0.0353 | 0.0347 | 0.0003 | 0.0681 | 0.0601
| true |  $<$0.0001 | 0.0351 | 0.0383 | 0.0004 | 0.0631 | 0.0669
20000 | bfal | 0.1328 | 0.0094 | 0.0101 | 0.2184 | 0.0189 | 0.0190
| gfal |  $<$0.0001 | 0.0094 | 0.0094 | 0.0002 | 0.0185 | 0.0162
| qfal |  $<$0.0001 | 0.0095 | 0.0107 | 0.0001 | 0.0190 | 0.0220
| super | 0.0002 | 0.0085 | 0.0093 |  $<$0.0001 | 0.0154 | 0.0162
| true |  $<$0.0001 | 0.0096 | 0.0096 | 0.0002 | 0.0187 | 0.0167
50000 | bfal | 0.1310 | 0.0046 | 0.0040 | 0.2158 | 0.0061 | 0.0076
| gfal |  $<$0.0001 | 0.0047 | 0.0038 | 0.0001 | 0.0059 | 0.0065
| qfal |  $<$0.0001 | 0.0046 | 0.0043 |  $<$0.0001 | 0.0060 | 0.0088
| super | 0.0001 | 0.0042 | 0.0038 | 0.0002 | 0.0060 | 0.0066
| true |  $<$0.0001 | 0.0046 | 0.0038 | 0.0001 | 0.0058 | 0.0067
Table 1: Results of simulation study, with $n_{sims}=500$, and $M=2$ (bias and
variance numbers multiplied by 100). $\widehat{V}_{sim}$ is the empirical
variance across simulation runs, and $\widehat{V}_{eif}$ is the mean of
influence curve-based variance estimates.  Figure 1: Simulation results.
Histograms are error estimates scaled by $\sqrt{n}$. Curves are normal density
functions with mean 0 and variance equal to the mean of estimated efficient
influence function-based variances across simulation runs.

##  4 Minimum wages and self-rated health

In this section, we apply our proposed approach to estimate the effect of
recent minimum wage increases on self-reported health using data from the
Behavioral Risk Factor Surveillence Survey (BRFSS). BRFSS is a cross-sectional
telephone survey focused conducted among more than 400,000 adults annually
[28]. BRFSS data are supplemented using annual state-level minimum wage data
and control variables related to state economic conditions from the U.S.
Department of Labor, Bureau of Labor Statistics, and Bureau of Economic
Analysis [31]. These data were originally compiled by [27] for use in a
similar difference-in-difference based analysis.

Our target estimand is the expected self-reported health for the US
population, had the US federal government counterfactually implemented a
minimum wage on par with that of New York State (NYS). Specifically, consider
a Likert scale measure of self-reported health ranging from 1-5 corresponding
to to self-rated “excellent”, “very good”, “good”, “fair”, and “poor” general
health. Let $Y_{st}$ denote the proportion of residents in a given state, $s$,
who report having “excellent” or “very good” health during year $t$. Let
$A_{st}$ be an indicator of whether state $s$ has minimum wage greater than or
equal to that of NYS in year $t$. Then the target estimand is the
intervention-specific mean $\mu_{t}=E[Y_{st}(\bar{a}^{*}_{t}=\bar{1})]$: the
expected proportion of individuals who would report excellent or very good
health if the US had a federal minimum wage law consistent with that of NYS.
If a given state has implemented a minimum wage that is greater than that of
NYS, then this state is consistent the hypothetical federal minimum wage and
is not intervened upon. Importantly, because the exposure for each state in a
given year explicitly depends upon the minimum wage in NYS, data can not be
considered as independent and identically distributed, which is an issue
because the proposed inferential techniques rely on this assumption. Instead,
we draw all conclusions conditional on knowledge of NYS’ minimum wage
trajectory over the period of study. Conditional on this information, data may
be considered IID and inference can proceed normally.

Our analysis is similar in scope to that of [27]. In their analysis, the
authors estimate the following two-way fixed effects model:

$E[Y_{istm}]=\beta_{0}+\beta_{1}MW_{s,t-1}+\beta_{2}Z_{itm}+\beta_{3}X_{st}+State_{s}+Year_{t}+Month_{m}+\phi_{st}+\epsilon_{istm}$

Where $MW_{st}$ denotes the minimum wage in state $s$ during year $t-1$,
$Z_{itm}$ denotes individual-level control variables, and $X_{st}$ denotes
state-level controls. Note that our approach differs from [27] in three
important ways: (1) the estimand of interest, (2) the identification and
modeling assumptions underlying our approach, and (3) the framing of the
minimum wage exposure. While our estimator targets the intervention-specific
mean, the target estimand of the TWFE model is a weighted-average of all
possible 2$\times$2 DID estimators [14], and is thus less immediately
interpretable. Moreoever, the TWFE approach relies on additional assumptions
for identification; in particular, that adjustment covariates $X_{st}$ are not
affected by previous treatment. This may not be plausible: state-level control
covariates such as unemployment rate and per capita income are very likely
affected by previous minimum wage changes.

We incorporate data from the 50 U.S. states and Washington, D.C.. Between 2013
and 2019, $n=3,111,090$ participants were surveyed. State-level confounding
variables considered in our analysis include economic indices such as the
maximum Temporary Assistance for Needy Families (TANF) benefit for a family of
four, the maximum Supplemental Nutrition Assistance Program (SNAP) benefit for
a family of four, the state Earned Income Tax Credit (EITC), the unemployment
rate, and the per capita personal income as well as state-level demographic
summaries of sex (proportion male), educational attainment (proportion
graduated high school, attended college or technical school, and graduated
from college or technical school), race (proportion non-hispanic white, non-
hispanic black, and other), and age (proportion older than 64). State-level
summaries were obtained by averaging the individual-level BRFSS data using
survey design weighting. Figure 3 presents the minimum wages for the 50 states
and Washington DC over the period from 2013 to 2019. In 2013, when the NYS
minimum wage was set to $7.25, all states had minimum wages at least as high
as NYS. By 2019, only 5 states (CA, DC, MA, VT, and WA) had maintained minimum
wages at least as high as NYS throughout the study period.

Estimating the target parameter required first estimating nuisance functions:
(1) the outcome models $Q^{j,k,m}$ for $j=2019,\ldots,2014$, $k=j\text{ or
}j+1$, and $m=2013,\ldots,j$ and (2) the exposure propensity model $g_{m}$ for
$m=2013,\ldots,2019$. Outcome models were estimated using SuperLearner, an
ensemble learning approach which models the outcome as an optimally weighted
combination of predictions from a set of candidate learners. Our candidate
learners included variations of random forests and multivariate adaptive
regression. For each type of learner for the outcome in year $j$, we
incorporated a version which used all observed covariates up to year $j$, as
well as versions which screened covariates based on having a significant
correlation with the selected outcome or based on being selected from a LASSO
regression. To avoid Donsker class restrictions, we used sample splitting
based on $K=2$ splits of the data. The cumulative propensity model, $g_{m}$,
was constructed using the product of time-specific treatment propensities,
$f(A_{k}=1|\bar{W}_{k},\bar{A}_{k-1}=\bar{1})$. These propensity models were
also fit using SuperLearner, with candidate learners including random forests
and elastic net regression. Note that treatment propensities only needed to be
modelled at times when the number of “consistent” states changed. For example,
from 2013 to 2014, all 51 states were consistent and no estimation of
propensity was required ($f(A_{2014}=1|\overline{W}_{2014},A_{2013}=1)=1$).
Models for treatment propensity had to be fit for the years $2015,\ 2016,\
2018$, and $2019$. Because the propensity score model conditions on having
remained consistent with NYS up to year $t-1$, small sample sizes in the final
years led to difficulties in estimating these models. For example, 8 states
were consistent with NYS minimum wage in 2018 and, of these states, 6 remained
consistent in 2019. In particular, when taking into account the sample
splitting required for both avoiding Donsker class conditions and inherent in
the optimization of SuperLearner model parameters were occasionally
unidentifiable. To address this concern, we assume the same propensity score
model for all years (i.e. $g_{m}=g$). In this simplified model, exposure in
the following year is predicted based on observed covariates over the previous
two years.

|  | Observed |  | Counterfactual |  | Difference
---|---|---|---|---|---|---
Year |  | Est. | SE | Lower | Upper |  | Est. | SE | Lower | Upper |  | Est. | SE | Lower | Upper
2015 |  | 52.0% | 0.66 | 50.7% | 53.3% |  | 51.6% | 0.78 | 50.1% | 53.2% |  | -0.40 | 0.35 | -1.09 | 0.28
2016 |  | 51.5% | 0.64 | 50.2% | 52.8% |  | 51.7% | 0.82 | 50.0% | 53.3% |  | 0.16 | 0.44 | -0.70 | 1.01
2017 |  | 50.3% | 0.67 | 49.0% | 51.6% |  | 50.7% | 0.84 | 49.1% | 52.4% |  | 0.42 | 0.50 | -0.55 | 1.40
2018 |  | 50.5% | 0.62 | 49.3% | 51.7% |  | 50.5% | 0.90 | 48.8% | 52.3% |  | 0.01 | 0.56 | -1.08 | 1.10
2019 |  | 49.7% | 0.62 | 48.5% | 50.9% |  | 50.4% | 0.86 | 48.8% | 52.1% |  | 0.75 | 0.64 | -0.51 | 2.00
Table 2: Estimated prevalence of “excellent” or “very good” self-rated health
in the US, based on data from the 2013-2019 Behavioral Risk Factor
Surveillance Study. Counterfactual estimates are based on our proposed
estimators implemented with Super Learner. Figure 2: Estimated proportion of
participants reporting ”Excellent” or ”Very Good” Health between 2015 and
2019. Proportions provided for the observed population and for a hypothetical
population of states with minimum wage laws consistent with NYS. 95% point-
wise confidence intervals included.

The primary results of this analysis can be found in Table 2 and Figure 2.
Table 2 presents the estimated proportion (with associated standard errors) of
individuals reporting “excellent” or “very good” among (1) the observed
population, and (2) a hypothetical population in which the federal government
had a minimum wage law consistent with that of NYS for each year between 2015
and 2019. Also provided are the estimated differences in proportions between
the two populations under each year of study. Figure 2 presents a step
function visualizing these results. Importantly, we do not observe sufficient
evidence to claim that setting the federal minimum wage to that of NYS would
have impacted the self-rated health of individuals during the years under
consideration. The largest estimated difference in proportion of individuals
with high-levels of self-reported health occured in the year 2019. In this
year, the observed proportion of participants who reporting high-levels of
health was 49.7% (95% CI: 48.5%, 50.9 %).

Figure 3: State-level minimum wage levels ($USD) between 2013 and 2019. New
York State minimum wage highlighted in bold.

##  5 Discussion

This paper proposed novel machine learning-based estimators for longitudinal
DID settings with time-varying treatments and/or covariates, focused on the
intervention-specific mean. By being based on the efficient influence
function, the proposed estimators are efficient and allow appropriate
inference even when the nuisance functions are estimated using complex machine
learning algorithms. In our simulations, the estimators achieved similar
finite-sample precision when using machine learning as when correctly
specified parametric models were used. Investigators rarely have knowledge
about the functional form of nuisance models, especially when these are higher
dimensional; thus interest in data adaptive estimators has been growing [10,
7, 6, 26]. Our focus on the intervention-specific mean stems from the fact
that these parameters can sometimes be highly policy relevant. For example,
[18] estimated counterfactual COVID-19 incidence rates if US-state-level
eviction moratoria had not expired; these results were then cited in the CDC
order calling for a federal eviction moratoria [11]. We also are among the
first DID papers to accomodate time-varying covariates; particularly in the
likely scenario that these may be affected by prior treatment. For example,
[4] discuss economic outcomes of pandemic-related policies. These policy
decisions are based on the current infection rates, which are in turn both
influenced by prior policy and influence future economic outcomes. Such
treatment-confounder, which we conjecture is operating in many contexts where
DID is used, would invalidate most existing DID methods but is adequately
addressed in our approach, so long as it is consistent with the parallel
trends assumption (the simulation setup in Section 3 provides such an
examples).

In contrast to exchangeability assumptions, the parallel trends assumption is
semiparametric and cannot be fully justified using causal diagrams; thus it
can be difficult to justify a priori. An active area of research aims to
describe, outside of any parametric model, what conditional independencies
and/or functional form restrictions in a structural model imply and/or are
implied by parallel trends [12, 20]. In parametric models, parallel trends
generally requires that any unmeasured confounders are time-invariant in both
their distribution and their association with outcomes over time. Thus, we may
reasonably expect that a semiparametric structural model consistent with
parallel trends would need to be partially linear with respect to the
unmeasured covariates, though there is more to be explored here. It is
important to note the scale-dependence of parallel trends; the version
presented in this paper is linear (as typical in DID), but the assumption
could be expressed on a nonlinear scale as well. A linear scale may not be
reasonable for many kinds of outcomes, such as survival times or other
distributions expected to be skewed.

Regarding estimation in practice, choosing the machine learning library and/or
the number of partitions in cross-fitting may be challenging. In general,
selecting a larger number of partitions reduces bias by ensuring nuisance
models are fit on a larger portion of the data, but can increase computation
time. Though our approach using a one-step estimator and (single-) cross-
fitting approach guarantees root-n consistency under the stated assumptions,
finite sample bias may be further reduced by applying other approaches. Recent
work by [35] describe a “double cross-fitting” approach that eliminates
second-order bias [21], wherein the data are repeatedly partitioned into three
components: (1) a training set for the outcome regression models, (2) a
training set for the exposure propensity model, and (3) an treatment effect
evaluation set. To our knowledge, this approach has not been extended to time-
varying settings such as ours. It may also be useful to develop target maximum
likelihood estimators [33, 17] for this setting, as the proposed approach does
not yield a substitution estimator and so could produce estimates outside the
known parameter range (e.g. risks outside of $[0,1]$), possibly contributing
to finite sample bias.

## References

  * [1] Alberto Abadie  “Semiparametric difference-in-differences estimators”  In _The review of economic studies_ 72.1 Wiley-Blackwell, 2005, pp. 1–19
  * [2] Joshua D Angrist and Jörn-Steffen Pischke  “Mostly harmless econometrics: An empiricist’s companion”  Princeton university press, 2009
  * [3] Carolina Caetano, Brantly Callaway, Stroud Payne and Hugo Sant’Anna Rodrigues  “Difference in differences with time-varying covariates”  In _arXiv preprint arXiv:2202.02903_ , 2022
  * [4] Brantly Callaway and Tong Li  “Policy evaluation during a pandemic”  In _Journal of Econometrics_ 236.1 Elsevier, 2023, pp. 105454
  * [5] Brantly Callaway and Pedro HC Sant’Anna  “Difference-in-differences with multiple time periods”  In _Journal of Econometrics_ 225.2 Elsevier, 2021, pp. 200–230
  * [6] Neng-Chieh Chang  “Double/debiased machine learning for difference-in-differences models”  In _The Econometrics Journal_ 23.2 Oxford University Press, 2020, pp. 177–191
  * [7] Victor Chernozhukov et al.  “Double/debiased machine learning for treatment and structural parameters”  Oxford University Press Oxford, UK, 2018
  * [8] Clément De Chaisemartin and Xavier d’Haultfoeuille  “Difference-in-differences estimators of intertemporal treatment effects”, 2022
  * [9] Iván Díaz  “Statistical inference for data-adaptive doubly robust estimators with survival outcomes”  In _Statistics in medicine_ 38.15 Wiley Online Library, 2019, pp. 2735–2748
  * [10] Iván Díaz, Nicholas Williams, Katherine L Hoffman and Edward J Schenck  “Nonparametric causal effects based on longitudinal modified treatment policies”  In _Journal of the American Statistical Association_ 118.542 Taylor & Francis, 2023, pp. 846–857
  * [11] Centers Disease Control and Prevention  “Temporary halt in residential evictions to prevent the further spread of COVID-19”  In _” ”_ , 2021
  * [12] Dalia Ghanem, Pedro HC Sant’Anna and Kaspar Wüthrich  “Selection and parallel trends”  In _arXiv preprint arXiv:2203.09001_ , 2022
  * [13] Dana E Goin and Corinne A Riddell  “Comparing Two-way Fixed Effects and New Estimators for Difference-in-Differences: A Simulation Study and Empirical Example”  In _Epidemiology_ 34.4 Wolters Kluwer, 2023, pp. 535–543
  * [14] Andrew Goodman-Bacon  “Difference-in-differences with variation in treatment timing”  In _Journal of econometrics_ 225.2 Elsevier, 2021, pp. 254–277
  * [15] James J Heckman, Hidehiko Ichimura and Petra E Todd  “Matching as an econometric evaluation estimator: Evidence from evaluating a job training programme”  In _The review of economic studies_ 64.4 Wiley-Blackwell, 1997, pp. 605–654
  * [16] Edward H Kennedy  “Semiparametric theory and empirical processes in causal inference”  In _Statistical causal inferences and their applications in public health research_ Springer, 2016, pp. 141–167
  * [17] Mark J Laan and Susan Gruber  “Targeted minimum loss based estimation of causal effects of multiple time point interventions”  In _The international journal of biostatistics_ 8.1 De Gruyter, 2012
  * [18] Kathryn M Leifheit et al.  “Expiring eviction moratoriums and COVID-19 incidence and mortality”  In _American journal of epidemiology_ 190.12 Oxford University Press, 2021, pp. 2503
  * [19] Alexander R Luedtke, Oleg Sofrygin, Mark J Laan and Marco Carone  “Sequential double robustness in right-censored longitudinal models”  In _arXiv preprint arXiv:1705.02459_ , 2017
  * [20] Philip Marx, Elie Tamer and Xun Tang  “Parallel trends and dynamic choices”  In _Journal of Political Economy Microeconomics_ 2.1 The University of Chicago Press Chicago, IL, 2024, pp. 129–171
  * [21] Whitney K Newey and James R Robins  “Cross-fitting and fast remainder rates for semiparametric estimation”  In _arXiv preprint arXiv:1801.09138_ , 2018
  * [22] Audrey Renson et al.  “Identifying and estimating effects of sustained interventions under parallel trends assumptions”  In _arXiv preprint arXiv:2206.05788_ , 2022
  * [23] James Robins  “A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect”  In _Mathematical modelling_ 7.9-12 Elsevier, 1986, pp. 1393–1512
  * [24] James Robins, Mariela Sued, Quanhong Lei-Gomez and Andrea Rotnitzky  “Comment: Performance of double-robust estimators when” inverse probability” weights are highly variable”  In _Statistical Science_ 22.4 JSTOR, 2007, pp. 544–559
  * [25] Andrea Rotnitzky, James Robins and Lucia Babino  “On the multiply robust estimation of the mean of the g-functional”  In _arXiv preprint arXiv:1705.08582_ , 2017
  * [26] Zach Shahn et al.  “Structural nested mean models under parallel trends assumptions”  In _arXiv preprint arXiv:2204.10291_ , 2022
  * [27] Liam Sigaud, Angela Daley, Jonathan Rubin and Caroline Noblet  “The effects of recent minimum wage increases on self-reported health in the United States”  In _Social Science & Medicine_ 305 Elsevier, 2022, pp. 115110
  * [28] Nicole M Silva  “The behavioral risk factor surveillance system”  In _The International Journal of Aging and Human Development_ 79.4 SAGE Publications Sage CA: Los Angeles, CA, 2014, pp. 336–338
  * [29] Tamar Sofer et al.  “On negative outcome control of unobserved confounding as a generalization of difference-in-differences”  In _Statistical science: a review journal of the Institute of Mathematical Statistics_ 31.3 NIH Public Access, 2016, pp. 348
  * [30] Liyang Sun and Sarah Abraham  “Estimating dynamic treatment effects in event studies with heterogeneous treatment effects”  In _Journal of econometrics_ 225.2 Elsevier, 2021, pp. 175–199
  * [31] U.S. Bureau of Labor Statistics  “Changes in basic minimum wages in non-farm employment under state law: selected years 1968 to 2023”  In _” ”_ Office of Communications, WageHour Division Washington, DC, 2024
  * [32] Aad W Vaart  “Asymptotic statistics”  Cambridge university press, 2000
  * [33] Mark J Van Der Laan and Daniel Rubin  “Targeted maximum likelihood learning”  In _The international journal of biostatistics_ 2.1 De Gruyter, 2006
  * [34] Bret Zeldow and Laura A Hatfield  “Confounding and regression adjustment in difference-in-differences studies”  In _Health services research_ 56.5 Wiley Online Library, 2021, pp. 932–941
  * [35] Paul N Zivich and Alexander Breskin  “Machine learning for causal inference: on the use of cross-fit estimators”  In _Epidemiology (Cambridge, Mass.)_ 32.3 NIH Public Access, 2021, pp. 393

##  Appendix A Proof of Theorem 1

By the linearity property of influence functions, it can be seen that:

$IF(\psi_{t})=IF(\mathbb{E}[Y_{0}])+\sum_{k=1}^{t}(IF(\phi_{k,k})-IF(\phi_{k-1,k})).$

Finding the influence function for $\psi_{t}$ simplifies to identifying the
influence function for $\phi_{j,k}$. By [17] this is given by:

$IF(\phi_{j,k})=\sum_{m=1}^{k}\frac{I(\overline{A}_{m}=\overline{a}_{m}^{*})}{g_{m}(\overline{W}_{m})}(Q^{j,k,m+1}-Q^{j,k,m})+Q^{j,k,1}-\phi_{j,k}$

To see why this is true, consider the following:

$\displaystyle IF(\phi_{j,k})$ $\displaystyle=IF(E[Q^{j,k,1}])$
$\displaystyle=IF\bigg{(}\int Q^{j,k,1}dF(\overline{w}_{1})\bigg{)}$
$\displaystyle=\int IF(Q^{j,k,1})dF(\overline{w}_{1})+\int
Q^{j,k,1}IF(dF(\overline{w}_{1}))$

Focus attention on the second term:

$\displaystyle\int Q^{j,k,1}IF(dF(\overline{w}_{1}))$ $\displaystyle=\int
Q^{j,k,1}IF(f(\overline{w}_{1}))d\overline{w}_{1}$ $\displaystyle=\int
Q^{j,k,1}(I(\overline{W}_{1}=\overline{w}_{1})-f(\overline{w}_{1}))d\overline{w}_{1}$
$\displaystyle=\int
Q^{j,k,1}I(\overline{W}_{1}=\overline{w}_{1})d\overline{w}_{1}-\int
Q^{j,k,1}f(\overline{w}_{1})d\overline{w}_{1}$
$\displaystyle=Q^{j,k,1}-E[Q^{j,k,1}]$ $\displaystyle=Q^{j,k,1}-\phi_{j,k}.$

Next, we will show that for $m=0,...,k$ the following relationship holds:

$\displaystyle\int
IF(Q^{j,k,m})\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=\overline{w}_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle=\int
IF(Q^{j,k,m+1})\prod_{s=0}^{m+1}dF(w_{s}|\overline{W}_{s-1}=\overline{w}_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle\qquad+\frac{I(\overline{A}_{m}=\overline{a}_{m}^{*})}{g_{m}(\overline{w}_{m})}(Q^{j,k,m+1}-Q^{j,k,m})$
(2)

Consider the following:

$\displaystyle\int IF(Q^{j,k,m})$
$\displaystyle\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle=\int IF\bigg{\\{}\int
Q^{j,k,m+1}dF(w_{m+1}|\overline{w}_{m},\overline{a}^{*}_{m})\bigg{\\}}\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle=\int\int
IF\bigg{\\{}Q^{j,k,m+1}\bigg{\\}}\prod_{s=0}^{m+1}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle\qquad+\int\int
Q^{j,k,m+1}IF\bigg{\\{}dF(w_{m+1}|\overline{w}_{m},\overline{a}^{*}_{m})\bigg{\\}}\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$

Focus on the second term:

$\displaystyle\int\int Q^{j,k,m+1}$ $\displaystyle
IF\bigg{\\{}dF(w_{m+1}|\overline{w}_{m},\overline{a}^{*}_{m})\bigg{\\}}\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle=\int\int
Q^{j,k,m+1}\bigg{\\{}\frac{I{(\overline{A}_{m}=\overline{a}_{m}^{*},\overline{W}_{m}=\overline{w}_{m})}}{g_{m}(\overline{w}_{m})\prod_{l=0}^{m}f(w_{l}|\overline{w}_{l-1},\overline{a}_{l-1}^{*})}(I(W_{m+1}=w_{m+1})-f(w_{m+1}|\overline{w}_{m},\overline{a}_{m}^{*}))dw_{m+1}\bigg{\\}}$
$\displaystyle\qquad\times\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle=\int\int
Q^{j,k,m+1}\bigg{\\{}\frac{I{(\overline{A}_{m}=\overline{a}_{m}^{*},\overline{W}_{m}=\overline{w}_{m})}}{g_{m}(\overline{w}_{m})\prod_{l=0}^{m}f(w_{l}|\overline{w}_{l-1},\overline{a}_{l-1}^{*})}I(W_{m+1}=w_{m+1})dw_{m+1}\bigg{\\}}\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle\qquad-\int\int
Q^{j,k,m+1}\bigg{\\{}\frac{I{(\overline{A}_{m}=\overline{a}_{m}^{*},\overline{W}_{m}=\overline{w}_{m})}}{g_{m}(\overline{w}_{m})\prod_{l=0}^{m}f(w_{l}|\overline{w}_{l-1},\overline{a}_{l-1}^{*})}f(w_{m+1}|\overline{w}_{m},\overline{a}_{m}^{*})dw_{m+1}\bigg{\\}}\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle=\int\frac{I{(\overline{A}_{m}=\overline{a}_{m}^{*},\overline{W}_{m}=\overline{w}_{m})}}{g_{m}(\overline{w}_{m})\prod_{l=0}^{m}f(w_{l}|\overline{w}_{l-1},\overline{a}_{l-1}^{*})}\bigg{\\{}\int
Q^{j,k,m+1}I(W_{m+1}=w_{m+1})dw_{m+1}\bigg{\\}}\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle\qquad-\int\frac{I{(\overline{A}_{m}=\overline{a}_{m}^{*},\overline{W}_{m}=\overline{w}_{m})}}{g_{m}(\overline{w}_{m})\prod_{l=0}^{m}f(w_{l}|\overline{w}_{l-1},\overline{a}_{l-1}^{*})}\bigg{\\{}\int
Q^{j,k,m+1}f(w_{m+1}|\overline{w}_{m},\overline{a}_{m}^{*})dw_{m+1}\bigg{\\}}\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle=\int\frac{I{(\overline{A}_{m}=\overline{a}_{m}^{*},\overline{W}_{m}=\overline{w}_{m})}}{g_{m}(\overline{w}_{m})\prod_{l=0}^{m}f(w_{l}|\overline{w}_{l-1},\overline{a}_{l-1}^{*})}Q^{j,k,m+1}(\overline{a}^{*},\overline{W}_{m+1}=(W_{m+1},\overline{w}_{m}))\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle\qquad-\int\frac{I{(\overline{A}_{m}=\overline{a}_{m}^{*},\overline{W}_{m}=\overline{w}_{m})}}{g_{m}(\overline{w}_{m})\prod_{l=0}^{m}f(w_{l}|\overline{w}_{l-1},\overline{a}_{l-1}^{*})}Q^{j,k,m}(\overline{a}^{*},\overline{w}_{m})\prod_{s=0}^{m}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle=\frac{I(\overline{A}_{m}=\overline{a}_{m}^{*})}{g_{m}(\overline{W}_{m})}(Q^{j,k,m+1}-Q^{j,k,m})$

This implies the overall influence function because:

$\displaystyle IF(\phi_{j,k})$ $\displaystyle=\int
IF(Q^{j,k,1})dF(\overline{w}_{1})+\int Q^{j,k,1}IF(dF(\overline{w}_{1}))$
$\displaystyle=\int
IF(Q^{j,k,k+1})\prod_{s=0}^{k+1}dF(w_{s}|\overline{W}_{s-1}=w_{s-1},\overline{A}_{s-1}=\overline{a}_{s-1}^{*})$
$\displaystyle\qquad+\sum_{m=1}^{k}\frac{I(\overline{A}_{m}=\overline{a}_{m}^{*})}{g_{m}(\overline{W}_{m})}(Q^{j,k,m+1}-Q^{j,k,m})$
$\displaystyle\qquad+Q^{j,k,1}-\phi_{j,k}$
$\displaystyle=\sum_{m=1}^{k}\frac{I(\overline{A}_{m}=\overline{a}_{m}^{*})}{g_{m}(\overline{W}_{m})}(Q^{j,k,m+1}-Q^{j,k,m})+Q^{j,k,1}-\phi_{j,k}$

Where the last equality follows from the fact that
$IF(Q^{j,k,k+1})=IF(Y_{j})=0$

##  Appendix B Proof of Theorem 2

By Lemma 1 of [19] we have

$\widetilde{\phi}_{jk}-\phi_{jk}=(P_{n}-P)IF(\phi_{jk})+(P_{n}-P)(\widehat{IF}(\phi_{jk})-IF(\phi_{jk}))+Rem_{jk},$

where $\widehat{IF}(\phi_{jk})$ denotes the IF of Theorem 1 with $\Theta$
replaced by $\widehat{\Theta}$, and where

$Rem_{jk}=\sum_{m=1}^{k}\int
c_{m}(\overline{w}_{m})\\{g_{m}(\overline{w}_{m})-\widehat{g}_{m}(\overline{w}_{m})\\}\\{Q^{j,k,m}(\overline{w}_{m})-\widehat{Q}^{j,k,m}(\overline{w}_{m})\\}dP(\overline{w}_{m}),$

where $c_{m}(\overline{W}_{m})$ is a constant given in the original Lemma,
which is bounded in probability according to C2. Under C1, Lemma 19.24 of [32]
shows that the term $(P_{n}-P)(\widehat{IF}(\phi_{jk})-IF(\phi_{jk}))$ is
$o_{P}(n^{-1/2})$. Likewise, using the Cauchy-Schwarz inequality we have

$|Rem_{jk}|\leq\left\lVert g_{m}-\widehat{g}_{m}\right\rVert\left\lVert
Q^{j,k,m}-\widehat{Q}^{j,k,m}\right\rVert.$

Applying the above to

$\psi_{t}=\mathbb{E}[Y_{0}]+\sum_{k=1}^{t}(\phi_{k,k}-\phi_{k-1,k});\quad\widetilde{\psi}_{t}=P_{n}Y_{0}+\sum_{k=1}^{t}(\widetilde{\phi}_{k,k}-\widetilde{\phi}_{k-1,k})$

together with Theorem 1 yields the result.