Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
86,004 characters · 21 sections · 46 citation commands
Identification of dynamic treatment effects when treatment histories are partially observed
\begingroup \footnote{$^\dag$Department of Econometrics and Business Statistics, Monash University. Emails: [email removed], [email removed]. We would like to thank Christian Cox, Chris Muris, Hidehiko Ichimura, Irene Botosaru, Hashem Pesaran, James Powell, Otavio Bartalotti, Tiemen Woutersen, and Wooyong Lee for comments.} \addtocounter{footnote}{-1} \endgroup
JEL Classification Codes: C14, C21, C23
Keywords: Parallel trends, Missing treatments, Panel data, Dynamic treatment effects, Robust, Difference-in-differences
Estimating dynamic treatment effects with panel data is often a central goal in applied research. Many empirical settings involve binary time-varying treatments (such as health shocks, medicare churning, or union membership) whose effects may persist well beyond the initial intervention period and depend on the full history of prior exposures. Difference-in-differences (DID) provides a compelling framework for studying such dynamics by leveraging repeated observations to control for unobserved time-invariant heterogeneity. However, identification of such path-dependent effects is challenging if a complete history of treatment decisions is not observed, whether that is due to i) survey non-response pepper2001response, ii) attrition in repeated surveys ghanem2024correcting, or iii) not observing some individuals in certain time periods, as in the case of rotating panels bellego2024chained. Standard approaches such as complete case analysis or imputation methods are valid only under relatively restrictive assumptions about the missingness mechanism or correct model specification. These conditions can often be violated in practice and lead to biased or inefficient estimates.
To illustrate this challenge, consider a stylized example of a balanced panel constructed from two non-consecutive survey waves. In this setting, a standard pre-post DID analysis that ignores the intermediate history generally fails to identify a causal parameter. We formally show that this estimand (which ignores persistence) identifies a non-convex weighted average of different path-dependent average treatment effects (PDATTs), which may not correspond to a causally meaningful quantity. This failure in identification is related to problems in the literature that discuss incorrect aggregation of heterogeneous causal effects.\footnote{See goodman2021difference, ishimaru2021we, callaway2021difference, sun2021estimating, imai2021use, de2020two, and others.} Moreover, the common alternative of relying on complete cases (CC), which essentially excludes observations with missing histories, only recovers specific PDATTs under strong assumptions that limit the heterogeneity of treatment paths by excluding adoption in period one, excluding dropouts or late-adopters, ruling out persistence, or imposing staggered adoption.\footnote{We show that a specific convex weighted average of PDATTs is partially identified under a monotone treatment response condition molinari2010missing.}
In this paper, we introduce a general framework for identifying path-dependent treatment effects in short panels when treatment histories are partially observed. Our setting allows for binary time-varying treatments which can switch on or off in each period, with no adoption in the initial period, and nests staggered adoption as a special case.\footnote{See roth2022s for a recent synthesis of the current DID literature with discussions on staggered adoption, violation of parallel trends, and design-based inference.} We develop a novel identification result for PDATTs under a missing-at-random selection mechanism and straightforward extensions of the DID identifying assumptions. The estimand adjusts for missing treatment histories by re-weighting observations based on the probability of experiencing a particular treatment path (propensity score) and the probability of it being observed in the population (missing data probability), combined with the conditional mean of outcomes (outcome regression). These elements are combined into a new augmented inverse probability weighted (AIPW) type estimand. Importantly, identification holds if any two of the three models involved - the outcome model, the propensity score model, or the missing treatment model - are correctly specified. Our identification result also nests missingness-adjusted versions of (i) outcome regression (OR), (ii) inverse probability weighting (IPW), and (iii) doubly robust (DR) estimands. However, each of these alternatives requires a correct missing data model along with at least one additional correctly specified model. In contrast, our proposed estimand identifies the target parameter even if the missing data model is misspecified.
Based on this identification result, we construct a two-step robust (R) estimator. The first step estimates the true nuisance functions (probability weights and outcome regression)\footnote{Our identification results are agnostic about the nature of the first-step estimators, and hence machine learning methods may be employed for estimating the three nuisance functions, especially in situations where large administrative datasets are available. Note that our theoretical results for inference are developed for parametric first-step estimators, and therefore do not cover the choice of machine learning methods or cross-fitting procedures.} and the second step plugs-in the estimated first-stage parameters into the sample analogue of the proposed estimand. We establish formal results on identification, estimation, and inference with the robust procedure. When all three models are correctly specified, the robust estimator attains the semiparametric efficiency bound for the PDATT parameter, making it efficient within the class of missingness-adjusted estimators. This result follows from our derivation of the associated efficiency bound. Moreover, since OR, IPW, and DR are nested within the robust proposal, inference with these methods is also made available.
Although the proposed method identifies the target parameter under misspecification of at-most one model, inference may still be affected. Specifically, the effect of estimating first-stage parameters indexing a misspecified model may propagate into the second step, altering the form of the asymptotic variance of the robust estimator. To reduce the effect of misspecification on inference, we further refine the asymptotic properties of the robust estimator by proposing two alternatives. Our approach builds on the recommendations in vermeulen2015bias, who propose minimizing the first-order effect of nuisance parameter estimation on a second-stage parameter of interest, and generalizes sant2020doubly's improved estimation results to settings with persistent treatment effects and/or missing treatment histories.
Numerical experiments help to evaluate the performance of the robust estimator against other missingness-adjusted estimators and CC-DID methods. First, we show that the robust estimator remains unbiased if either the missing data model, propensity score model, or outcome regression model is misspecified. In contrast, DR, IPW, and OR are biased whenever the missing data model is misspecified, irrespective of the specification of the other nuisance functions. Second, inference based on the robust estimator has accurate test size across all experiments, whereas the DR, IPW, and OR estimators show considerable size distortions when the missing data model is misspecified. Third, experiments varying the extent of missingness and the degree of misspecification in the missing data model reveal that the bias in CC-DR and DR estimators grows as missingness rates or misspecification severity increases. In contrast, the robust estimator continues to perform well.
We conclude by demonstrating the practical relevance of the robust estimator with two empirical applications. The first investigates the persistent effects of COVID-19 cases on county-level voter turnout in the 2022 U.S. general elections, where case histories are missing for 51% of the counties. The robust estimator suggests a statistically significant reduction in voter turnout of 0.18% points for counties that experienced above-average number of cases in 2020 and 2021, while standard CC-DID estimates suggest a negligible and statistically insignificant reduction in voter turnout between 0.01% and 0.03% points. The second application uses individual-level data from the Current Population Survey (CPS) to study the effects of worker disability, job certification, and work absence on family income and hours worked. While missingness in treatment histories is modest (under 5%), a meta-analysis across the three treatments reveals that the robust estimator can yield estimates that differ substantially from those obtained using CC-DR methods.
\paragraph{Relation to the literature:} We contribute to a growing body of literature which allows outcomes to be affected by the entire treatment path. An early example is hull2018estimating, who studies two-way-fixed-effects regressions for mover panels and imposes some version of conditional mean impersistence. strezhnev2018semiparametric develops inverse propensity score weighted DID estimators for estimating persistent treatment effects with multiple time periods. de2022two and de2024difference extend their earlier work on interpretations of two-way-fixed-effects regressions to allow for several treatments and treatment lags, respectively. viviano2021dynamic propose a dynamic covariate balancing method for estimating the effects of different treatment trajectories. None of these papers address the challenge of missing treatment histories.
In the panel data literature, our paper is broadly related to the strand studying missing covariates or treatments. abrevaya2017gmm present a generalized method of moments (GMM) estimator for dealing with missing regressors. muris2020efficient provides a GMM framework for efficient parameter estimation with incomplete data and botosaru2018difference propose a proxy-variable solution to address the problem of a missing treatment variable within a standard DID analysis with repeated cross sections. Finally, coe2019estimation proposes an inverse probability weighted solution for pooled ordinary least squares and first-differenced moments.
There is also a rich literature on robust estimation of treatment effects.\footnote{See robins1994estimation, scharfstein1999adjusting, graham2012inverse, bang2005doubly, sloczynski2018general, lewbel2023over, negi2024doubly for doubly robust estimators in cross-sectional settings.} In the context of panel data, \citet*{arkhangelsky2021double} develop an augmented doubly robust two-way-fixed effects estimator and arkhangelsky2022doubly integrate design-based and model-based identification strategies to construct a doubly robust alternative. Our paper is closely related to the papers by sant2020doubly (SZ) and callaway2021difference (CS), who propose doubly robust estimators for ATTs in simple and staggered adoption settings, respectively. Our PDATT estimator equals the proposal in SZ and CS in special cases. yanagi2022doubly generalizes SZ and CS to allow for general treatment patterns across multiple time periods. Implicitly, these papers assume that treatments are fully observed. Recently, bellego2024chained propose a chained DID method that combines short-term treatment effects from many incomplete unbalanced panels to estimate long-run effects in staggered settings.\footnote{In the high-dimensional literature, farrell2015robust introduces a doubly robust estimator for constructing confidence intervals for the ATE after model selection and chernozhukov2022locally propose locally robust orthogonal moment conditions that also exhibit doubly robust properties.}
The statistics literature explores multiply robust estimation of average treatment effects. This spans mediation analysis \citep*{xia2023identification, tchetgen2014estimation, jiang2022multiply}, missing outcomes han2014multiply, han2013estimation, and missing treatment information \citep*{zhang2016causal}. \citet*{shi2020multiply} propose a multiply robust ATE estimator in the presence of categorical unmeasured confounding and negative controls, while wang2018bounded use instrumental variables, and \citet*{wei2023multiply} study nonrandom assignment and missing outcomes. zhang2016causal examine robust estimation of ATEs in a cross-sectional setting with missing treatment data with a binary outcome, and propose an estimator which exhibits properties similar to ours under model misspecification.
The rest of the paper is organized as follows. Section (ref) introduces our framework, the parameters of interest, and the identifying assumptions. Section (ref) presents the identification, estimation, and inferential results with the proposed approach. Section (ref) discusses other missingness-adjusted estimands that, while less robust, are nested within our theoretical results. Section (ref) presents numerical experiments comparing the different estimators and Section (ref) illustrates the estimators in two empirical applications. Section (ref) concludes.
Let $Y_t$ be the observed outcome at time period $t$, and $D_t$ be a binary treatment which is equal to one if an individual is treated in period $t$ or zero otherwise. Assume that there is no treatment in the baseline with $D_0=0$. Additionally, we observe a $k$-dimensional vector of pre-treatment characteristics $\mathbf{X}$.\footnote{It is standard in the DID literature to only consider time-invariant covariates or to condition only on the pre-treatment values for any time-varying covariates (see callaway2019quantile and references therein).} For ease of exposition, consider a setting with three time periods denoted by $t=0,1,2$. The treatment history is then denoted by $\mathbf{D}=(D_1,D_2)$ and we define $\Delta Y= Y_2-Y_0$. Empirically, treatment histories may be partially observed. To formalize this, let $S$ be a binary indicator which is equal to one if $D_1$ is observed and zero otherwise. Extension of this framework to a general short panel setting with an arbitrary number of time periods and general missing treatment history patterns is discussed in Section (ref).
A parameter that is of interest to policy makers is the effect of a particular treatment history on final period ($t=2$) outcomes. We define the average effect of experiencing treatment path $\mathbf{D} = \mathbf{d}$ compared to $\mathbf{d}^\prime$ for individuals who experienced path $\mathbf{d}$ in period $t=2$ as
where $Y_t(\mathbf{d})$ denotes the potential outcome in period $t$ if the treatment history $\mathbf{D}$ takes the value $\mathbf{d}=(d_1, d_2)\in \{0,1\}^2$. Our parameters of interest always consider $\mathbf{d^\prime} = (0,0)$.\footnote{Supplementary Appendix \ref*{sec:cpt_alt} discusses identification of $\tau_{\mathbf{dd^\prime}}$ which involves comparisons with $\mathbf{d^\prime} = (1,1)$. Such comparisons would necessitate assuming parallel trends in the treated counterfactual distribution - an assumption that is practically never invoked in the literature.} The observed outcome is $Y_t = Y_t(\mathbf{D})$. With two treatments, the definition in (ref) allows for three PDATTs. Summaries of these effects may also be of interest. For instance, the average effect of receiving the second treatment ($D_2=1$) compared to not receiving any treatment equals $\tau_{(11)(00)}\mathbb{P}(D_1=1|D_2=1)+\tau_{(01)(00)}\mathbb{P}(D_1=0|D_2=1)$.
Our setup covers a wide range of policy-relevant treatment settings, including both sequential and simultaneous interventions. In an educational context, $D_1$ and $D_2$ could represent college enrollment and graduation, respectively. Here, $\tau_{(11)(00)}$ measures the effect of the whole college program, $\tau_{(10)(00)}$ measures the impact of enrollment for dropouts, and $\tau_{(01)(00)}$ captures the effect of graduation for late-adopters. nibbering2024instrument show that in general, treatment programs with dropouts and late enrollment are ubiquitous. Similarly, in a workforce development program katz2022sectoral, $D_1$ and $D_2$ might represent an initial job training program followed by an internship, with PDATTs defined analogously.
With staggered treatment adoption, such as an irreversible implementation of state-level minimum wage laws callaway2021difference, PDATTs will capture differential effects of early ($\tau_{(11)(00)}$) versus late adoption ($\tau_{(01)(00)}$). One can also use this framework to study simultaneous treatments that may be correlated in time de2023two. For instance, with labor market policies, minimum wage regulations and working hours restrictions may be implemented together, making it important to jointly account for them when studying outcomes of interest. The PDATTs remain relevant here for capturing heterogeneous treatment effects corresponding to simultaneous, and potentially, interactive policies.
In order to identify the PDATTs in (ref), we extend the difference-in-differences assumptions to allow potential outcomes to depend on the full history of treatment decisions.
Assumption (ref).1 rules out any anticipatory effects of future treatment on outcomes at $t=0$. Violations arise if individuals change their behavior in anticipation of the treatment. Assumption (ref).2 imposes that the average trend in the untreated potential outcome of the treatment and comparison groups would have evolved in parallel between $t=0$ and $t=2$, conditional on $\mathbf{X}$. This is commonly referred to as conditional parallel trends. Assumption (ref).3 is an overlap or common support condition which bounds the propensity score $p_\mathbf{d}(\mathbf{X})$ away from one.
A natural starting point when only partial information on $D_1$ is available, is to completely ignore $D_1$ and conduct a conditional pre/post DID analysis using the first and final time period, while varying $D_2$. As we show below, this strategy fails to identify a causal parameter.
Proof is deferred to Appendix (ref). The conditional DID estimand which ignores $D_1$ identifies a non-convex weighted average of three PDATTs, with weights given by different conditional treatment probabilities. This estimand does not have a causal interpretation unless one is willing to impose additional assumptions on treatment adoption. For instance, if dropouts and late-adopters are ruled out, and $D_1=D_2$, the estimand identifies $\tau_{(11)(00)}$. When treatment in period $t=1$ can be ruled out i.e.\ $D_1=0$, we recover $\tau_{(01)(00)}$. Under the assumption of treatment impersistence ($Y_2(0,d_2)=Y_2(1,d_2)$ for each $d_2$) or the assumption of staggered adoption ($D_2\geq D_1$), we identify the weighted average $\tau_{(11)(00)}\cdot \mathbb{P}(D_1=1|D_2=1)+\tau_{(01)(00)}\cdot \mathbb{P}(D_1=0|D_2=1)$. Hence, this DID approach only recovers specific causal parameters under strong assumptions, but cannot identify individual PDATTs that are permissible in the general case.\footnote{Supplementary Appendix \ref*{sec:partial_DID} shows that a convex weighted average of specific PDATTs is partially identified under a monotone treatment response assumption.}
Participation in treatments may be missing for a variety of reasons. First, missingness may arise due to item non-response where survey questions elicit information about sensitive behaviors such as drug use and alcohol consumption pepper2001response. Second, treatment participation may only be reported partially, thereby obscuring whether individuals adhered to the full treatment program, dropped-out, or adopted late silliman2022labor, zimmerman2014returns. For example, information about when treatment was initiated may be available, but actual adoption timing might be unknown. Additionally, treatment data may also be missing for individuals due to noncompliance with the assigned treatment. Finally, when panel data are constructed from repeated surveys, attrition can pose a significant problem ghanem2024correcting. If the attrited sample is systematically different from the observed sample, attrition bias can distort treatment effect estimates. We impose the following assumptions on the missing treatment mechanism.
Assumption (ref).1 is a novel missing-at-random (MAR) assumption tailored to our DID setting with missing treatments. It permits missingness in $D_1$ to be correlated with the fully-observed treatment and covariates, and subsumes both the stronger version of MAR, $S \perp (\Delta Y, \mathbf{D})|\mathbf{X}$, and the missing completely at random (MCAR) version, $S \perp (\Delta Y, \mathbf{D}, \mathbf{X})$.\footnote{Supplementary Appendix \ref*{sec:weakmar} discusses an alternative MAR assumption that allows missingness to be correlated with $\Delta Y$, explains why existing estimands and those proposed in this paper are biased in this case, and proposes a novel estimand that remains unbiased under certain conditions.} Assumption 2.2 ensures that for each group defined by $(D_2,\mathbf{X})$, there is a positive probability of observing $D_1$.
From Assumption (ref).1, it follows that $\mathbbm{E}[\Delta Y|\mathbf{D}=\mathbf{d}, \mathbf{X}, S=s] = \mathbbm{E}[\Delta Y|\mathbf{D}=\mathbf{d}, \mathbf{X}]$, which has implications for the potential outcomes. For the comparison group with $\mathbf{D}=\mathbf{0}$, this imposes conditional parallel trends between the observed and unobserved comparison groups: $\mathbbm{E}[Y_2(\mathbf{0})-Y_0(\mathbf{0})|\mathbf{D}=\mathbf{0}, \mathbf{X}, S=s] = \mathbbm{E}[Y_2(\mathbf{0})-Y_0(\mathbf{0})|\mathbf{D}=\mathbf{0}, \mathbf{X}]$. This assumption is made on the outcome trends instead of levels, and therefore the latter can still depend on the missingness mechanism. For the treatment groups, we have $\mathbbm{E}[Y_2(\mathbf{d})-Y_2(\mathbf{0})+Y_2(\mathbf{0})-Y_0(\mathbf{0})|\mathbf{D}=\mathbf{d}, \mathbf{X}, S=s] = \mathbbm{E}[Y_2(\mathbf{d})-Y_2(\mathbf{0})+Y_2(\mathbf{0})-Y_0(\mathbf{0})|\mathbf{D}=\mathbf{d}, \mathbf{X}]$ which requires (i) conditional parallel trends between the observed and unobserved treated groups and (ii) conditional independence between the treatment effects and the missingness mechanism. In particular, Assumption (ref) can be violated if treatment effects vary with the missingness mechanism even after conditioning on observables, despite conditional trends being the same across observed and unobserved groups. shin2024difference discusses similar assumptions in a DID setting with missing outcomes.
A common empirical strategy to deal with missing data is to restrict the analysis to the subset of observations for whom treatment histories are fully observed, also known as complete-case (CC) analysis. In the current setting, this entails conducting the DID analysis on the set of observations for whom $S=1$. While this approach is simple and avoids the need for imputation or weighting adjustments, it will produce inconsistent estimates of PDATT unless the missingness mechanism satisfies the stricter MCAR assumption. The following proposition provides an explicit expression of selection bias introduced by this method within our setup.
Proof is in Appendix (ref). Proposition (ref) shows that in the presence of missing treatments ($\mathbb{P}(S=0|\mathbf{D}=\mathbf{d})\neq 0$), the CC estimand is biased if the covariate distributions for the groups experiencing treatment path $\mathbf{d}$ are different between the observed and unobserved subpopulations: $\mathbb{P}(\mathbf{X}|\mathbf{D}=\mathbf{d},S=1)\neq\mathbb{P}(\mathbf{X}|\mathbf{D}=\mathbf{d},S=0)$. Such selection bias arises due to the fact that the PDATTs require the integration over $\mathbf{X}$ under Assumption (ref).2, despite the fact that $\Delta Y$ does not depend on $S$ given $D_2$ and $\mathbf{X}$ under Assumption (ref).1. This bias disappears if $D_1$ is MCAR, but the efficiency loss from discarding incomplete observations may be substantial.
For ease of exposition, we consider a setting with three time periods and a partially missing $D_1$ throughout this paper. Our main results hold in a setting with a general number of time periods in which we allow for general missing treatment history patterns. We briefly discuss this setup here, and defer the details to Appendix (ref).
First, we extend our framework to $1<T<<n$ time periods with $n$ the number of individuals or units. The treatment history is denoted by $\mathbf{D} = (D_1, \ldots, D_T)$, and we maintain that no unit receives treatment in $t=0$. The PDATT in (ref) now generalises to $\tau_{\mathbf{dd^\prime}} = \mathbbm{E}[Y_T(\mathbf{d}) - Y_T(\mathbf{d}^\prime)|\mathbf{D} = \mathbf{d}]$, with $Y_T=Y_T(\mathbf{d})$ and $\mathbf{d}^\prime = (0,\ldots,0) \equiv \mathbf{0}_T$.
Second, we generalize our missing treatment mechanism as follows. Partition $\mathbf{D}$ into two vectors $\mathbf{D}_{-h}$ and $\mathbf{D}_{h}$, such that the elements in $\mathbf{D}_{-h}$ are observed with probability one and the elements in $\mathbf{D}_{h}$ may be missing. With $T=2$ and $D_1$ partially missing, this boils down to $\mathbf{D}_{-h}=D_2$ and $\mathbf{D}_{h}=D_1$. However, we can also capture $D_2$ missing by setting $\mathbf{D}_{-h}=D_1$ and $\mathbf{D}_{h}=D_2$. With $T>2$, multiple time periods of the treatment history may be missing, from which follows that $\mathbf{D}_{h}$ may include multiple time periods.
The binary indicator $S$ now indicates whether all elements in $\mathbf{D}_h$ are observed, and the missing at random assumption imposes that this indicator is independent of $\mathbf{D}_{h}$ and $\Delta Y=Y_T-Y_0$, given $\mathbf{D}_{-h}$ and $\mathbf{X}$. When $T$ is large, a stronger assumption may be invoked to make estimation of a missing data model feasible. For instance, one that only requires conditioning on time periods adjacent to the ones in $\mathbf{D}_{h}$ instead of all time periods in $\mathbf{D}_{-h}$. Alternatively, more flexible missingness patterns can be allowed by defining separate missingness indicators for each time period in $\mathbf{D}_{h}$.
We consider three models corresponding to the true unknown outcome means, propensity scores, and missing treatment probabilities, respectively. More precisely, $\mu_{\mathbf{d}}(\mathbf{X})$ represents a model for the outcome mean $m_{\mathbf{d}}(\mathbf{X}) \equiv \mathbb{E}[\Delta Y|\mathbf{D}=\mathbf{d},\mathbf{X}]$. The models $\pi_{d_1|d_2}(\mathbf{X})$ and $\pi_{d_2}(\mathbf{X})$ represent the propensity scores $p_{d_1|d_2}(\mathbf{X})\equiv \mathbb{P}({D_1}={d_1}|D_2=d_2,\mathbf{X})$ and $p_{d_2}(\mathbf{X})\equiv \mathbb{P}({D_2}={d_2}|\mathbf{X})$, respectively. Finally, $\phi_{d_2}(\mathbf{X})$ is a model for the missing treatment probability, $q_{d_2}(\mathbf{X})= \mathbb{P}(S=1|D_2=d_2,\mathbf{X})$. Our first result shows how correctly specified $\mu_{\mathbf{d}}(\mathbf{X})$ and $\pi_{d_1|d_2}(\mathbf{X})$ can be identified under the MAR assumption, even when $D_1$ is missing.
The proof is deferred to Appendix (ref). The expressions on the right-hand sides of the equations in Lemma (ref) only depend on observables which implies that outcome models and propensity score models are identified. We combine these models with the missing treatment probability models into a single estimand. This estimand identifies the PDATT in (ref) even if one of the three models is misspecified. This robust estimand is given as
where the hajek1971discussion-type weights are defined as {
} The following result states the robustness property of the estimand:
Proof is deferred to Appendix (ref). The intuition behind this result follows from the fact that when any two of the three models are replaced by their true counterparts, certain components of the estimand —best described as adjustment terms— vanish in expectation, and the remaining term equals the true parameter. Consider the following decomposition of the estimand:
The proof of Theorem (ref) shows that with correct propensity score and outcome models, (I)-(II)=(VI) and (III)-(IV)=0. The term (V) is only a function of propensity score and outcome models, and provided that these models are correctly specified, identifies the target parameter:
In case one of the correctly specified models is the missing treatment model $\phi_{d_2}(\mathbf{X})$, the proof of Theorem (ref) shows that (V)=(VI). If, in addition, the outcome model is correct we have (III)-(IV)=0 and (I)-(II) identifies the PDATT, or the propensity score model is correct and (II)-(IV)=0 and (I)-(III) identifies the PDATT:
Corollaries (ref) and (ref) directly follow from the proof of Theorem (ref). Note that all results in this section are derived in Appendix (ref) for the general case discussed in Section (ref).
The main novelty of the estimand in Theorem (ref) is that it enables identification of PDATTs with partially observed treatment histories, even with misspecification in the missing treatment model. We highlight this contribution with three examples. First, consider a setting in which the missing treatment model depends on the observed covariates $\mathbf{X}$ in an unknown way. If Assumption (ref) holds, and the propensity score and outcome models are correctly specified in $\mathbf{X}$, the robust estimand would identify the PDATTs. Second, suppose that the missingness mechanism depends on a different set of covariates than those required for the conditional parallel trends. Let $\mathbf{X}=(\mathbf{X}_1,\mathbf{X}_2)$, where Assumption (ref), the propensity score, and the outcome models depend on $\mathbf{X}_1$, but the missing treatment model depends on $\mathbf{X}_2$. In such case, the researcher only requires knowledge on how the propensity scores and the outcome models vary with $\mathbf{X}_1$ to identify the PDATTs. Third, consider a situation where missingness is driven by unobserved factors. If Assumption (ref), the propensity score, and outcome models depend on $\mathbf{X}$, while Assumption (ref) depends on unobservables, we still achieve identification.
To investigate the conditions under which the robust estimand is efficient, we first derive the semiparametric efficiency bound for the PDATT parameter in (ref) in the presence of missing treatments. The semiparametric efficiency bound serves as a benchmark for the asymptotic variance of any $\sqrt{n}$-consistent estimator of $\tau_{\mathbf{dd^\prime}}$. In spirit, one can think of this as the semiparametric analogue of the Cramer-Rao lower bound for parametric models.
Proof is deferred to Appendix (ref). The derivation of the bound for the data $(Y_2, Y_0, \mathbf{D}, \mathbf{X})$ is self-contained and can be seen to follow previous results in the literature (see for example, hahn1998role and sant2020doubly). From there on, we employ the result in Theorem 7.2 in tsiatis2006semiparametric to derive the bound under our MAR assumption.
The expression in (ref) suggests that the robust estimand can be estimated with a two-step procedure, provided that a random sample is available.
Assumption (ref) covers a setting in which panel data are available.\footnote{The $i.i.d$ assumption can be relaxed to allow for intra-cluster correlations in cases where data have a clustering dimension. Our identification and estimation results will continue to hold in such case, while inference will have to be adjusted to account for such correlation structure.} Estimation can then proceed as follows. First, the models for the true unknown outcome means, propensity scores, and missing data probabilities are estimated. Second, the predicted values for these estimated models are plugged into the sample analogue of $\tau_{\mathbf{dd^\prime}}^{\textup{R}}$.
The first step requires a choice of models and estimators for the outcome means, propensity scores, and missing data probabilities. So far, we have simply postulated the existence of models for each of these functions but have not committed to it either being parametric or non-parametric in nature. We derive the asymptotic behavior of the estimator for $\tau_{\mathbf{dd^\prime}}^{\textup{R}}$ assuming parametric first-stage estimators, which allows us to derive asymptotic theory for general parametric estimators. These estimators are often preferred in applied work due to their simplicity, and due to the fact that nonparametric estimators may suffer from challenges such as the curse of dimensionality or tuning parameter selection.
Let $\mu(\bm{\beta}_\mathbf{d})$, $\pi(\bm{\gamma}_{\mathbf{d}})$, and $\phi(\bm{\delta}_{d_2})$ be parametric models for $m_{\mathbf{d}}(\mathbf{X})$, $p_{\mathbf{d}}(\mathbf{X})$, and $q_{d_2}(\mathbf{X})$, respectively, where we suppress the dependence of these models on data for notational convenience. Define the pseudo-true parameter values as $\bm{\beta^\ast}_{\mathbf{d}}$, $\bm{\gamma^\ast}_{\mathbf{d}} = (\bm{\gamma^\ast}_{d_1|d_2}, \bm{\gamma^\ast}_{d_2})$, and $\bm{\delta^\ast}_{d_2}$. Let $\bm{\widehat{\beta}}_{\mathbf{d}}$, $\bm{\widehat{\gamma}}_{\mathbf{d}}$, $\bm{\widehat{\delta}}_{d_2}$ denote $\sqrt{n}$-consistent estimators of these pseudo-true values. The estimator of the robust estimand $\widehat{\tau}^\textup{R}_{\mathbf{dd^\prime}}$ is given by
where $\mathbb{E}_n(\cdot)$ denotes the empirical mean and the weights are estimated as {
} with $\bm{\widehat{\gamma}}=(\bm{\widehat{\gamma}}_{\mathbf{d}},\bm{\widehat{\gamma}}_{\mathbf{d^\prime}})$, and the dependence of these weights on the data is suppressed. Define $\bm{\beta^\ast}=(\bm{\beta^\ast}_{\mathbf{d}},\bm{\beta^\ast}_{\mathbf{d^\prime}})$, $\bm{\gamma^\ast}=(\bm{\gamma^\ast}_{\mathbf{d}},\bm{\gamma^\ast}_{\mathbf{d^\prime}})$, and $\bm{\delta^\ast}=(\bm{\delta^\ast}_{d_2}$, $\bm{\delta^\ast}_{d_2^\prime})$. Theorem (ref) derives the asymptotic properties of $\widehat{\tau}^\textup{R}_{\mathbf{dd^\prime}}$ using some weak high-level conditions on the estimators for the generic parametric models, which are outlined in Appendix (ref):
Proof is deferred to Appendix (ref). Theorem (ref) shows that $\widehat{\tau}^\textup{R}_{\mathbf{dd^\prime}}$ is $\sqrt{n}$-consistent and asymptotically normal provided that at least two of the three models are correct. This result suggests that we can use the sample analogue to $\Omega$ to conduct asymptotically valid inference. Estimation of the nuisance parameters affects the asymptotic variance of the robust estimator. This effect is proportionate to the average change in the influence function of the robust estimator from locally perturbing the first-stage parameters around their probability limits. When these probability limits index a correctly specified population model, small changes in $(\bm{\beta}, \bm{\gamma}, \bm{\delta})$ have no effect on the influence function of the robust estimator, causing the estimation effect from the first stage to disappear. In the special case when all three models are correctly specified, we show that $\widehat{\tau}^\textup{R}_{\mathbf{dd^\prime}}$ achieves the semiparametric efficiency bound.
Proof is deferred to Appendix (ref). A practical implication of Corollary (ref) is that when all three models are correct, the choice of first-step estimators does not influence the asymptotic variance of the robust estimator. However, this property is lost as soon as one of the working models is misspecified. In this case, the expression for $\Omega$ in Theorem (ref) includes terms that depend on the first-stage estimators, making inference sensitive to such choice. Supplementary Appendix \ref*{sec:inferencerobust} explores two inference-robust alternatives whose asymptotic variance remains unaffected under misspecification of any one model. The search for such an alternative is inspired from sant2020doubly who propose improved DID estimators in the standard DID setup with similar robustness properties. Building on the insights from vermeulen2015bias, our inference-robust proposals use first-stage estimators that are specifically designed to minimize the effect of the nuisance parameters on the robust estimator.
Corollary (ref) presents estimands based on a correct missing data model along with either a correct propensity score or mean outcome model thereby giving us missingness-adjusted outcome regression (OR) and inverse probability weighting (IPW) estimands. These are given by
and
It is important to note that unlike the standard OR method, which only depends on a correct outcome model, the missingness-adjusted OR estimand given in (ref) depends on both a correct missing treatment and outcome model. In a similar spirit, the adjusted IPW estimand in (ref) depends not only on a correct propensity score but also a correct missing treatment model, thereby requiring both probability weights to be correct to identify the target parameter.
For our setting, we can also combine the two results in Corollary (ref) to give us the missingness-adjusted DR estimand, which is given by
It follows from Theorem (ref), and the discussion around the decomposition in (ref), that this estimand identifies $\tau_{\mathbf{dd^\prime}}$ when either the missing data model and the outcome model are correct, or the missing data model and the propensity score model are correct. While our proposed approach (R) is robust to misspecification in the missing data model, the missingness-adjusted OR, IPW, and DR strategies presented above will not identify the target parameter if the missing data model is incorrect, making it a less preferred alternative compared to R.
The proof follows directly from Corollary (ref) and Theorem (ref) for OR, IPW, and DR estimands, respectively. Replacing the population models in (ref)-(ref) with their estimated counterparts allows us to propose estimators which are given by
Since these are incomplete versions of the robust estimator, their asymptotic variances can be easily and directly obtained from the variance of $\widehat{\tau}^\textup{R}_{\mathbf{dd^\prime}}$. Supplementary Appendix \ref*{sec:simsadd} presents the asymptotic influence function representations for all three alternatives.
In this section, we first conduct a Monte Carlo study which analyzes the finite sample performance of different estimators for PDATTs. We then study how varying amounts of missingness and degrees of misspecification in the missing data model affect their performance.
The data generating process is defined as
where $U_1$, $U_2$, and $U_3$ are three independently distributed random variables with a standard uniform distribution, and $\varepsilon$ has a standard normal distribution.
The covariates in the propensity scores, missing treatment probabilities, and outcome means are denoted by $\mathbf{X}_p$, $\mathbf{X}_m$, and $\mathbf{X}_o$, respectively. We set $\mathbf{X}_g=\eta_g\mathbf{X}+(1-\eta_g)\mathbf{{Z}}$ with $\eta_g=0,1$ and $g=p,m,o$. The vector $\mathbf{X}$ includes an intercept and four independently distributed standard normal random variables $X_1, \ldots, X_4$. We then use the transformations defined in kang2007demystifying: $\tilde{Z}_1= \text{exp}(0.5 X_1)$, $\tilde{Z}_2= 10+X_2/(1+\text{exp}(X_1))$, $\tilde{Z}_3 = \left(0.6+X_1 X_3/25\right)^3$ and $\tilde{Z}_4 = \left(20+X_2+X_4\right)^2$. This gives us the vector $\mathbf{Z}$ which includes an intercept and $\tilde{Z}_1, \ldots, \tilde{Z}_4$ that are standardized to have mean 0 and variance 1. Since $\mathbf{X}$ is treated as the vector of observed covariates, setting $\eta_g=0$ results in a misspecified working model. The values for the parameters in (ref) are provided in Table (ref). The percentage of missing values for $D_1$ is governed by $c$, where $c=0$ corresponds to approximately 50% missingness.
We estimate the PDATTs $\tau_{(11)(00)}$, $\tau_{(10)(00)}$, and $\tau_{(01)(00)}$ using the estimators $\widehat{\tau}_{\mathbf{dd^\prime}}^{\textup{R}} $, $\widehat{\tau}^\textup{OR}_{\mathbf{dd^\prime}}$, $\widehat{\tau}^\textup{IPW}_{\mathbf{dd^\prime}}$, and $\widehat{\tau}^\textup{DR}_{\mathbf{dd^\prime}}$ defined in (ref), (ref), (ref), and (ref), respectively. The asymptotic distributions of the OR, IPW, and DR estimators are derived in Supplementary Appendix \ref*{sec:simsdet}. Additionally, we estimate the PDATTs using the CC-OR, CC-IPW, and CC-DR estimators which rely on the observed sample, without adjusting for missing histories through a missing data model. All estimators use working models as specified in Appendix (ref).
We study the finite sample performance of the estimators with four Monte Carlo experiments. Each experiment consists of 10,000 replications with each a sample of $(\Delta Y_i, S_i, S_iD_{1i},D_{2i},\mathbf{X}_i)$ with 10,000 observations generated from (ref) with $c=0$. In these experiments, either only the missing data model (M), only the propensity score (P), only the outcome regression (OR), or none of the models are misspecified (None). The four scenarios in which more than one working model is misspecified are discussed in Supplementary Appendix \ref*{sec:simsresults}. Since each scenario conditions on different covariates, the values for the PDATTs vary across the experiments; ${\tau}_{(11)(00)}$ varies between 1.68 and 1.71, ${\tau}_{(10)(00)}$ between 0.58 and 0.63, and ${\tau}_{(01)(00)}$ between 1.15 and 1.42.
Figure (ref) shows the bias of the missingness-adjusted estimators for the PDATTs across the four experiments. We find that the proposed robust estimator (R) is unbiased across all settings under consideration. For the other estimators, there is a bias when the missing data model is misspecified. As expected, IPW also shows bias when the propensity score model is misspecified, OR when the outcome regression is misspecified, and all missingness-adjusted estimators are unbiased when all models are correctly specified.
For the CC estimators, the bias becomes substantial compared to the bias caused by misspecified working models in the case of missingness-adjusted estimators. This indicates that selection bias can have a relatively large impact on the accuracy of the PDATT estimates. Detailed results for these estimators are reported in Supplementary Appendix \ref*{sec:simsresults}.
Figure (ref) shows the test size of testing the null-hypothesis that a PDATT equals its true value at a nominal level of 5% across the four experiments. We find that the robust estimator obtains nominal test size for all PDATTs across all experiments. The tests corresponding to the other estimators are oversized when the missing data model is misspecified. In addition, we find major distortions for the OR and IPW estimators if the outcome regression or propensity score models are incorrect, respectively. When all models are correctly specified, all four estimators appropriately control size. In this case, the asymptotic variance of R is close to the semiparametric efficiency bound. Supplementary Appendix \ref*{sec:simsresults} provides the estimates for the asymptotic variances of the estimators together with the semiparametric efficiency bounds.
Figure (ref) shows the statistical power of the test of $H_0:\tau_{(11)(00)}=0$ with nominal test size of 5%. The panels correspond to experiments with different misspecified models. Each panel shows the power curves of the estimators that theoretically should control size under the misspecification at hand by dotted lines, and the power curves in the experiments with none of the models misspecified by solid lines. Note that the power curves of R and DR are identical in the second and third panel, and hence the latter are not displayed. We find that the power curves of the R estimator are close to the curves of the other estimators across all experiments. This indicates that the potential power loss of using the most robust estimator relative to the OR, IPW, and DR estimators is small. Moreover, since the power curves in the experiments with only correctly specified models are close to the power curves with misspecified models, the power loss due to misspecification also seems small.
Second, we illustrate the importance of appropriately accounting for missing treatment histories at varying rates of missingness. The percentage of missingness in the Monte Carlo experiments equals approximately 50%. By varying $c$, the intercept in the missing data models, we explore how the proportion of missingness influences the bias of different estimators. We assume that only the missing data model is misspecified with $\eta_p=\eta_o=1-\eta_m=1$ and generate one million observations from (ref) with $c \in \{-5,-4.5,\dots,4.5,5\}$.
Figure (ref) shows the bias in R, DR, and CC-DR for an increasing percentage of missingness. We find that the estimates from R have negligible bias, even when the amount of missingness is large. However, both DR and CC-DR show a bias that increases in the percentage of missingness. These biases follow from a misspecified missing data model in DR or from sample selection in CC-DR, and may already affect estimates when only a small percentage of treatment histories is missing.
Third, we examine the effect of different degrees of misspecification in the missing data model on the bias of different estimators. Since $\eta_m=1$ corresponds to a correctly specified model, we can explore the effect of an increasing amount of misspecification by decreasing $\eta_m \in \{1,0.9,\dots,0.1,0\}$. We assume that only the missing data model is misspecified with $\eta_p=\eta_o=1$ and $c=0$, and generate one million observations from (ref) for each value of $\eta_m$.
Figure (ref) shows the bias in R, DR, and CC-DR for an increasing degree of misspecification. Again, we find negligible bias in the estimates from R across all degrees. The bias in DR increases in the amount of misspecification, but is small compared to the bias in CC-DR, which does not necessarily depend on the degree of misspecification. Both findings align with our theoretical results, which show that the DR estimator requires the missing data model to be correct and CC-DR requires the sample selection bias to be negligible. The first directly hinges on the degree of misspecification, while the latter also depends on other features of the data generating process. The illustrations in Figures (ref) and (ref) show that minor violations of these assumptions can already cause substantial biases in these PDATT estimators.
In this section, we demonstrate the proposed method by applying it to two distinct empirical settings. The first application investigates the effects of covid case surges on voter turnout in the 2022 U.S. general elections.\footnote{Numerous studies have investigated the effects of COVID cases and COVID-led policies on a range of outcomes callaway2023evaluating, callaway2023policy, kim2021impact, morgenstern2022interaction, badinlou2024trajectories, reuschke2024impacts,herrnson2023impact.} This setting involves aggregated county-level data and has treatment histories missing for around 50% of the sample. The second application conducts a comprehensive meta-analysis utilizing individual-level household data from CPS to examine the effects of worker disability, job certification, and work absenteeism on family income and hours worked. Treatment histories are generally missing here for less than 5% of the sample.
We obtain daily county-level COVID-19 transmission data from the CDC from 2020 to 2022.\footnote{Following callaway2023evaluating, we obtain data on weekly number of covid cases from \url{https://data.cdc.gov/Public-Health-Surveillance/United-States-COVID-19-County-Level-of-Community-T/8396-v7yb/about_data}.} This is combined with county-level voter turnout data between 2004-2022 from the National Neighborhood Data Archive at Inter-university Consortium for Political and Social Research. We supplement these data with county-level covariates from the US Census Bureau and United State's Department of Agriculture Economic Research Service database.
Our binary treatment measures whether, during a given month, a county's weekly number of confirmed cases per 100,000 people ever exceeds the national weekly average in a given year. We refer to our treatment variable as “above-average” cases in a particular year. To illustrate our method, we consider August across three years; 2018 is the pre-treatment period, 2020 is the middle period, and 2021 is the final period.\footnote{While we could have used 2022 to be the last treatment period, COVID cases had declined significantly that year with mass vaccination already underway.} In the data, 51% of the counties were missing confirmed cases data in at least one week of August 2020 (middle period). Even though the last treatment period is 2021, turnout rates are only available in 2022 since there were no general elections in the previous year. Our outcome is defined as changes in county-level voter turnout between the first and final period. The final sample has 3,096 counties.
We condition on county-level covariates\footnote{See Supplementary Appendix \ref*{sec:emp_sc} for additional details.} that can account for differential trends in voter turnout. Turnout patterns are typically stable over time at the county level, and the timing of case surges is plausibly exogenous to underlying electoral dynamics. Since covid cases are likely to be missing due to factors like public health infrastructure and reporting practices, which could be correlated with the county characteristics in the observed covariates, our MAR assumption is also plausible in this setting.
Table (ref) reports the estimated effects of having above-average number of cases on turnout rates in the 2022 general elections along with standard errors and estimated confidence intervals. Based on the robust method, we find that having above-average cases in 2020 and 2021 reduces turnout rates for counties that experienced it by 0.18% points, on average. Unlike the estimates obtained using adjusted-DR or CC methods, this estimate is statistically significant. In general, the CC-OR, CC-IPW, and CC-DR estimates are smaller than their adjusted counterparts with narrower confidence intervals.
We use publicly available household survey data from the CPS which is accessed through the Integrated Public Use Microdata Series (IPUMS). The CPS is a nationally representative monthly survey conducted jointly by the U.S. Census Bureau and the Bureau of Labor Statistics, and serves as the official source of labor force statistics for the U.S. population.
We construct a three-period panel by using the monthly observations within a given year with the household head (HH) as the unit of analysis. We define the initial and final periods as the first and final months a household is observed. Treatment in the middle period is defined as whether the HH receives the treatment during the intermediate month(s). Disability, job certification, and work absence treatment status may be missing in the middle period for various reasons: the household may not have been surveyed during those months due to CPS rotation design, responses may be missing due to item or unit non-response, or responses may be unknown (which CPS codes as NIU). Overall, missingness in the middle period remains low, affecting fewer than 5% of the samples.\footnote{See Supplementary Appendix \ref*{sec:emp_sc} for additional details and sample construction.} The CPS also collects extensive demographic information on the household members which includes region, race, sex, marital status, level of education, nativity, which are standardized before before being used in estimation. When missingness is uncorrelated with treatment status in the middle time period and income or hours worked, conditional on the covariates, our MAR assumption holds in this setting.
Based on the robust estimates, we find PDATTs align with economic intuition and vary across time. For example, in 2009, having a disability in both periods reduced hours worked by 3.559 hours (with a standard error of 1.149), while being disabled in only one period in 2009 does not have a statistically significant effect. In contrast, having job certification in both periods or the second period in 2017 has a statistically significant positive effect on family income: the effect of job certification in both periods is \$668.652, only in the first period is \$148.938, and only in the second period is \$391.694, with standard errors equal to 297.751, 225.406, and 73.537, respectively. For work absence, we find estimates indicating that only the effect in the final period is significant: being absent from work in both periods increases family income in 2009 by \$20.421 (13.576), while the effects of absence in the first and second period equal a reduction in income of \$10.949 (7.441) and \$78.985 (17.976).
The differences between the robust estimates and the estimates from the DR and CC-DR estimators can be substantial. Table (ref) reports the mean, median, and maximum values of the absolute percent differences in the PDATT estimates of these methods across all outcome-by-year combinations for each treatment variable. Consider $\tau_{(01)(00)}$ for the disability treatment. On average, the R estimate is 16% larger than the DR estimate with a maximum percent difference of 57%, and 36% larger than the CC-DR estimate with a maximum percent difference of 97%. For job certification, the mean differences between R and DR and CC-DR for $\tau_{(01)(00)}$ are 1.5% and 18%. Similarly, the mean differences for absence are smaller compared to disability. The maximum differences show that the estimates from R can be very different from CC-DR, with percentage differences reaching up to 157%.
In this paper, we consider a difference-in-differences framework with a binary time-varying treatment, no treated units in the pre-treatment period, but otherwise no restrictions on treatment-path heterogeneity. We identify and estimate the effect of each treatment history on final period outcomes with treatment histories partially observed. We propose a novel AIPW estimand which identifies the target parameter as long as any two of the outcome, propensity score, or missing treatment models are correctly specified. This method generalizes and improves upon other missingness-adjusted alternatives (such as IPW, OR, and DR) which require the missing treatment model to be correctly specified alongside another correctly specified component.
We present numerical experiments which compare the performance of the missingness-adjusted estimators and their complete-case counterparts. We find that the robust estimand remains unbiased and controls size across the three cases of model misspecification whereas the other adjusted estimators exhibit bias and size distortions once the missing treatment model is misspecified. By varying the degree of missingness and misspecification, we show that the bias in R remains negligible compared to the bias in DR and CC-DR estimators. We further demonstrate the applicability of the missingness-adjusted methods compared to the practice of dropping data through two empirical applications. First, we find an economically and statistically significant treatment effect of covid cases on voter turnout across U.S.\ counties in the presence of 51% missingness using the proposed estimator. Second, a meta-analysis on CPS household data shows that the proposed method can produce estimates very different from existing methods in a wide range of settings even with missingness below 5%.
\singlespacing