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.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Efficient Difference-in-Differences Estimation when Outcomes are Missing at Random
abstractThe Difference-in-Differences (DiD) method is a fundamental tool for causal inference, yet its application is often complicated by missing data. Although recent work has developed robust DiD estimators for complex settings like staggered treatment adoption, these methods typically assume complete data and fail to address the critical challenge of outcomes that are missing at random (MAR) -- a common problem that invalidates standard estimators. We develop a rigorous framework, rooted in semiparametric theory, for identifying and efficiently estimating the Average Treatment Effect on the Treated (ATT) when either pre- or post-treatment (or both) outcomes are missing at random. We first establish nonparametric identification of the ATT under two minimal sets of sufficient conditions. For each, we derive the semiparametric efficiency bound, which provides a formal benchmark for asymptotic optimality. We then propose novel estimators that are asymptotically efficient, achieving this theoretical bound. A key feature of our estimators is their multiple robustness, which ensures consistency even if some nuisance function models are misspecified. We validate the properties of our estimators and showcase their broad applicability through an extensive simulation study.
Introduction
The Difference-in-Differences (DiD) method stands as a cornerstone of modern applied research in economics lechner2011estimation, social sciences greene2021review, and public health wing2018designing, prized for its intuitive and powerful approach to estimating the causal effects of policies and interventions from observational panel, or longitudinal, data. In its canonical form, the method is used to estimate the average treatment effect on the treated (ATT), by leveraging a core identifying condition, the parallel trends assumption. This assumption posits that, had the treatment not occurred, the outcomes for the treated and control groups would have followed similar (parallel) trajectories. This core assumption allows researchers to use the observed change in the control group's outcome to construct the counterfactual path of the treated group, thereby isolating the treatment's causal impact by differencing out confounding factors that are constant over time.
While the classic two-group, two-period setup provides a clear theoretical foundation, empirical applications are rarely so simple. The last decade has seen a critical reexamination of how DiD methods are applied in more complex settings -- see callaway2023difference, de2023credible, roth2023s for a review. A major focus has been on scenarios with multiple time periods and variation in when different units receive treatment, known as \enquote{staggered adoption} callaway2021difference, de2020two, goodman2021difference, sun2021estimating. A second parallel line of inquiry has developed methods to handle missingness in post-treatment outcomes, often under stringent assumptions. For example, rathnayake2024difference, shin2024difference, viviens2025difference provide estimation strategies for the ATT in an unconditional framework where available covariates are disregarded; while bellego2025chained, working in a staggered adoption setting, assumes that outcomes are observed at least twice for each unit. However, despite the progress made, none of these papers discusses efficiency in either nonparametric or semiparametric models in two-group, two-period DiD setups with available covariates. A further persistent challenge, commonly encountered in practice, is how to conduct DiD analysis when pre-treatment outcome data are missing for a subset of the sample. This situation is common in many real scenarios: in labor economics, individuals entering a job training program may lack pre-treatment earnings records; in education policy, students may transfer into a school district without their prior test scores; and in health research, electronic health records may not contain a baseline measurement for patients who did not have a clinical visit during the pre-period.
When the pre-treatment outcome is unavailable for some units, the standard DiD estimator cannot be computed. A simple \enquote{complete-case} analysis that discards observations with missing data can be deeply flawed, as it may introduce significant selection bias. In fact, if the mechanism that causes the data to be missing is related to treatment or other characteristics that influence the outcome, the remaining sample is no longer representative of the population of interest and the resulting estimates will be biased.
This paper addresses this critical gap by bridging modern DiD estimation with semiparametric statistical theory bickel1993efficient, kennedy2024semiparametric, tsiatis2006semiparametric. We provide a formal framework for identifying and efficiently estimating the ATT when pre-treatment outcomes are missing at random (MAR). In the Appendix, we also extend our results to the symmetric scenario where all pre-treatment outcomes are observed, but post-treatment outcomes can be missing at random.
Finally, again in the Appendix, we also further extend our framework to accommodate the scenario in which outcomes can be missing both before and after treatment.
Problem setup
We assume a two-group, two-period DiD setup\footnote{The framework can be extended to staggered adoption settings as in callaway2021difference, but we work in a two-group, two-period setting for ease of exposition.}. In particular, we assume to have access to a collection of $n$ i.i.d. observed-data samples $\mathcal{D}_i = (X_i, R_{i0}, R_{i0} Y_{i0}, A_i, Y_{i1})\sim\mathbbmss{P}^\star$, $i=1,\dots,n$, where $X_i\in\mathbbmss{R}^p$ is a $p$-dimensional vector of pre-treatment covariates; $R_{i0}\in\{0,1\}$ is a binary variable indicating whether the baseline, pre-treatment outcome $Y_{i0}$ is observed ($R_{i0}=1$) or not ($R_{i0}=0$); $A_i\in\{0,1\}$ is a binary variable indicating whether observation $i$ has received a treatment ($A_i=1$) or was in the control group ($A_i=0$); finally $Y_{i1} = A_i Y_{i1}^{(1)} + (1-A_i) Y_{i1}^{(0)}$ is the observed outcome after treatment, where $Y_{i1}^{(1)}$ and $Y_{i1}^{(0)}$ denote the potential outcomes for observation $i$ under treatment and control, respectively. We let $\mathcal{D}=(X,R_0,R_0Y_0,A,Y_1)$, with $Y_1= A Y_1^{(1)} + (1-A)Y_1^{(0)}$, denote an independent copy of $\mathcal{D}_i$. From a theoretical viewpoint, it is also useful to construct a prototypical full-data sample as $\mathcal{D}^F = (X,Y_0,Y_1^{(0)}, Y_1^{(1)})$ -- see tsiatis2006semiparametric.
remarkWe want to emphasize that our framework differs from both the balanced panel data and the repeated cross-section data analyzed by sant2020doubly. The former postulates that each sample is observed both before and after treatment: the latter assumes that each sample is observed either before or after treatment. Instead, our setup mirrors the so-called unbalanced, or partially missing, panel data framework, where the outcomes of some samples are observed both before and after the treatment, while some other samples miss pre-treatment outcomes.
Our goal is the estimation of the average treatment effect on the treated (ATT), defined as:
equation[equation omitted — 276 chars of source]
The ATT is a function of the full data $\mathcal{D}^F$, and as such is not directly identifiable using only observed data $\mathcal{D}$. In the next Section, we provide a minimal set of assumptions that make the target $\theta^\star$ a function of the observed data $\mathcal{D}$.
Identification and efficiency bounds
We provide a first set of assumptions for identifiability.
assumption[Identifiability]
Let the following identifiability assumptions hold:
\begin{enumerate}[label=\alph*.]
• Conditional parallel trends. $Y^{(0)}_1 - Y_0 \perp \!\!\!\! \perp A \mid X$.
• Consistency. $Y_1= A Y_1^{(1)} + (1-A)Y_1^{(0)}$.
• Positivity. $\pi^\star(x) = \mathbbmss{P}\left[A=1\mid X=x\right]\in(0,1)$ almost surely for every $x\in\mathbbmss{R}^p$.
\end{enumerate}
remark{\color{black} Assumption a is stronger than the conditional mean assumption standardly employed in the econometrics literature on DiD usually expressed as:
\begin{equation}
\mathbbmss{E}\left[Y^{(0)}_1 - Y_0 \mid X, A=1\right] = \mathbbmss{E}\left[Y^{(0)}_1 - Y_0 \mid X, A=0\right]\,.
\end{equation}
We use the stronger conditional independence notation in Assumption a for \enquote{notational consistency} because it aligns with the standard conditional independence language used in semiparametric theory and in our subsequent MAR assumptions.} Importantly, Assumption a does not imply conditional independence of either $Y^{(0)}_1$ or $Y_0$, but only of their incremental difference. Assumptions b and c are standard in causal inference. In particular, Assumption b guarantees the observed outcome at time $t=1$ corresponds to the potential outcome associated with the treatment status. Assumption \textbf{c} forces the probability of being treated, for each observation, to be strictly positive.
Then, we also provide two different sets of assumptions for the missingness mechanism in the pre-treatment outcome.
assumption[Outcome independent missing at random]
Let the following MAR assumptions hold:
\begin{enumerate}[label=\alph*.]
• No unmeasured confounding. $Y_0 \perp \!\!\!\! \perp R_0 \mid X,A$.
• Weak overlap. $\gamma^\star(x,a) = \mathbbmss{P}\left[R_0 = 1 \mid X=x, A=a\right] \in (0,1)$ almost surely for every $x\in\mathbbmss{R}^p$ and $a\in\{0,1\}$.
\end{enumerate}
assumption[Outcome dependent missing at random]
Let the following MAR assumptions hold:
\begin{enumerate}[label=\alph*.]
• No unmeasured confounding. $Y_0 \perp \!\!\!\! \perp R_0 \mid X,Y_1,A$.
• Weak overlap. $\gamma^\star(x,y_1,a) = \mathbbmss{P}\left[R_0 = 1 \mid X=x, Y_1=y_1, A=a\right] \in (0,1)$ almost surely for every $x\in\mathbbmss{R}^p$, $y_1\in\mathbbmss{R}$, and $a\in\{0,1\}$.
\end{enumerate}
example[Medical records motivation]
To illustrate the practical distinction between these assumptions, consider a longitudinal study on the progression of a chronic disease, where $Y_0$ is the baseline health status and $Y_1$ is the post-treatment status. Under Assumption (ref), missingness in baseline records ($R_0$) is driven solely by covariates ($X$) or treatment assignment ($A$). For instance, missingness might arise because a specific hospital system ($X$) used paper records that were not digitized, regardless of how sick the patients currently are. Under Assumption (ref), missingness may be driven by the future outcome. For example, in retrospective chart reviews, the availability of baseline data ($R_0$) often depends on the patient's current condition ($Y_1$). Clinicians may be required to diligently hunt down and record historical data ($Y_0$) only for patients who currently exhibit severe complications ($Y_1$). Conversely, for patients who have recovered (low $Y_1$), the baseline charts may never be retrieved or digitized ($R_0=0$). In this setting, conditioning on $Y_1$ is necessary to block the dependence between missingness and the unobserved baseline outcome. See Figure (ref) for a graphical representation using DAGs of the causal relationships under the two different assumptions.
figure[figure omitted — 243 chars of source]
remarkThese assumptions are novel to our setting and are equivalent to a missing at random (MAR) missingness design. Notice that we let the missingness pattern depend on both covariates $X$ and the treatment $A$, and eventually on the post-treatment outcome. In other words, we admit the possibility that pre-treatment outcomes can be missing due to covariates, the treatment that will be administered, and the post-treatment outcome values.
Equipped with the previous sets of assumptions, we can now identify the ATT as a function of the observed data $\mathcal{D}$ as shown in the following Lemma.
lemma[Identification of ATT]
Under Assumption (ref) and Assumption (ref), the ATT can be identified as a function of the observed data:
\begin{equation}
\begin{split}
\theta^\star &= \mathbbmss{E}\left[\frac{A}{\mathbbmss{E}\left[A\right]} \left( Y_1 - \mathbbmss{E}\left[Y_0\mid X,A=1, R_0=1\right]\right)\right]}\left[Y_0\mid X,A=1, R_0=1\right]\right)} \\
&- \mathbbmss{E}\left[\frac{A}{\mathbbmss{E}\left[A\right]} \left(\mathbbmss{E}\left[Y_1\mid X, A=0\right] - \mathbbmss{E}\left[Y_0\mid X, A=0, R_0=1\right] \right)\right]ght] - \mathbbmss{E}\left[Y_0\mid X, A=0, R_0=1\right] \right)}\,.
\end{split}
\end{equation}
Under Assumption (ref) and Assumption (ref), the ATT can be identified as a function of the observed data:
\begin{equation}
\begin{split}
\theta^\star &= \mathbbmss{E}\left[\frac{A}{\mathbbmss{E}\left[A\right]} \left( Y_1 - \mathbbmss{E}\left[Y_0\mid X,Y_1,A=1, R_0=1\right]\right)\right]ft[Y_0\mid X,Y_1,A=1, R_0=1\right]\right)} \\
&- \mathbbmss{E}\left[\frac{A}{\mathbbmss{E}\left[A\right]} \left(\mathbbmss{E}\left[Y_1\mid X, A=0\right] - \mathbbmss{E}\left[\mathbbmss{E}\left[Y_0\mid X,Y_1, A=0, R_0=1\right]\mid X, A=0\right]=1\right]\mid X, A=0} \right)\right]\left[\mathbbmss{E}\left[Y_0\mid X,Y_1, A=0, R_0=1\right]\mid X, A=0\right]=1\right]\mid X, A=0} \right)}\,.
\end{split}
\end{equation}
Once identified, the next natural question is how well the ATT can be estimated under each set of assumptions. Before proceeding, we denote the regression functions as $\mu^\star_t(x,a) = \mathbbmss{E}\left[Y_t \mid X=x, A=a\right]$ for $t\in\{0,1\}$, and $\mu^\star_0(x,y_1,a) = \mathbbmss{E}\left[Y_0 \mid X=x, Y_1=y_1, A=a\right]$. Furthermore, denote the nested regression function as $\eta^\star_t(x,a) = \mathbbmss{E}\left[\mu^\star_t(x,Y_1,a)\mid X=x,A=a\right]$. Notice that we cannot use the law of iterated expectations to simplify $\eta^\star_t(x,a)$ into $\mathbbmss{E}\left[Y_0\mid X=x,A=a\right]$. The inner conditional expectation $\mu^\star_t(x,Y_1,a)$ is in fact implicitly conditioning on $R_0=1$, while the external conditional expectation is not. With this in mind, we suppress dependence of $\mu^\star_t(x,Y_1,a)$ on $R_0=1$ for notational simplicity. The following Proposition provides the semiparametric efficiency bounds.
proposition[Semiparametric efficiency bounds]
Under Assumption (ref) and Assumption (ref), the efficient observed-data influence function is given by:
\begin{equation}
\begin{split}
\varphi\left(\mathcal{D}\right) &= \frac{A}{\mathbbmss{E}\left[A\right]}\left(Y_1 - \left(\mu^\star_0(X,1) + \frac{R_0}{\gamma^\star(X,1)}\left(Y_0 - \mu^\star_0(X,1)\right) \right) - \mu^\star_1(X,0) + \mu^\star_0(X,0)\right) \\
&- \frac{(1-A)\pi^\star(X)}{(1-\pi^\star(X))\mathbbmss{E}\left[A\right]} \left(Y_1 - \frac{R_0}{\gamma^\star(X,0)}\left(Y_0 - \mu^\star_0(X,0)\right) - \mu^\star_1(X,0) \right) - \frac{A}{\mathbbmss{E}\left[A\right]}\theta^\star\,.
\end{split}
\end{equation}
Under Assumption (ref) and Assumption (ref), the efficient observed-data influence function is given by:
\begin{equation}
\begin{split}
\varphi\left(\mathcal{D}\right) &= \frac{A}{\mathbbmss{E}\left[A\right]}\left(Y_1 - \left(\mu^\star_0(X,Y_1,1) + \frac{R_0}{\gamma^\star(X,Y_1,1)}\left(Y_0 - \mu^\star_0(X,Y_1,1)\right) \right) - \mu^\star_1(X,0) + \eta^\star_0(X,0)\right) \\
&- \frac{(1-A)\pi^\star(X)}{(1-\pi^\star(X))\mathbbmss{E}\left[A\right]} \left(Y_1 - \left(\mu^\star_0(X,Y_1,0) + \frac{R_0}{\gamma^\star(X,Y_1,0)}\left(Y_0 - \mu^\star_0(X,Y_1,0)\right) \right) - \mu^\star_1(X,0) + \eta^\star_0(X,0) \right) \\
&- \frac{A}{\mathbbmss{E}\left[A\right]}\theta^\star\,.
\end{split}
\end{equation}
The semiparametric efficiency bound under each assumption is given by $\mathbbmss{V}\left[\varphi\left(\mathcal{D}\right)\right]=\mathbbmss{E}\left[\varphi^2\left(\mathcal{D}\right)\right]$.
remarkWhen the pre-treatment outcome is always observed, we recover the results by sant2020doubly (Proposition 1) and hahn1998role (Theorem 1). In fact, the efficient influence function in this case simplifies to
\begin{equation}
\begin{split}
\varphi^F\left(\mathcal{D}\right) &= \frac{A}{\mathbbmss{E}\left[A\right]}\left(Y_1 - Y_0 - \mu^\star_1(X,0) + \mu^\star_0(X,0)\right) \\
&- \frac{(1-A)\pi^\star(X)}{(1-\pi^\star(X))\mathbbmss{E}\left[A\right]} \left(Y_1 - Y_0 - \mu^\star_1(X,0) + \mu^\star_0(X,0) \right) - \frac{A}{\mathbbmss{E}\left[A\right]}\theta^\star\,.
\end{split}
\end{equation}
Notice, however, that by allowing for missingness, our approach provides a better description of the loss in efficiency incurred when pre-treatment outcomes are missing at random. In particular, the efficiency loss incurred under Assumption (ref) is:
\begin{equation}
\begin{split}
\mathbbmss{V}\left[\varphi\left(\mathcal{D}\right)\right] - \mathbbmss{V}\left[\varphi^F\left(\mathcal{D}\right)\right] &= \mathbbmss{V}\left[\varphi\left(\mathcal{D}\right) - \varphi^F\left(\mathcal{D}\right)\right] \\
&= \mathbbmss{E}\left[\frac{\pi^\star(X)\mathbbmss{V}\left[Y_0\mid X,A=1\right]}{\mathbbmss{E}\left[A\right]^2} \frac{1-\gamma^\star(X,1)}{\gamma^\star(X,1)}\right]}{\gamma^\star(X,1)}} \\
&+
\mathbbmss{E}\left[\frac{\pi^\star(X)^2 \mathbbmss{V}\left[Y_0\mid X,A=0\right]}{(1-\pi^\star(X))\mathbbmss{E}\left[A\right]^2} \frac{1-\gamma^\star(X,0)}{\gamma^\star(X,0)} \right]{\gamma^\star(X,0)} }\,,
\end{split}
\end{equation}
where we are exploiting the fact that the influence function without missingness is a projection of the observed-data influence function, and as such the variance of the observed-data influence functions decomposes nicely due to the Pythagorean theorem tsiatis2006semiparametric.
Similarly, the efficiency loss incurred under Assumption (ref) is:
\begin{equation}
\begin{split}
\mathbbmss{V}\left[\varphi\left(\mathcal{D}\right)\right] - \mathbbmss{V}\left[\varphi^F\left(\mathcal{D}\right)\right] &= \mathbbmss{E}\left[\frac{\pi^\star(X) \mathbbmss{V}\left[Y_0\mid X,Y_1,A=1\right]}{\mathbbmss{E}\left[A\right]^2} \frac{1-\gamma^\star(X,Y_1,1)}{\gamma^\star(X,Y_1,1)}\right]ammatarget(X,Y_1,1)}} \\
&+ \mathbbmss{E}\left[\frac{\pi^\star(X)^2 \mathbbmss{V}\left[Y_0\mid X,Y_1,A=0\right]}{(1-\pi^\star(X))\mathbbmss{E}\left[A\right]^2} \frac{1-\gamma^\star(X,Y_1,0)}{\gamma^\star(X,Y_1,0)}\right]ammatarget(X,Y_1,0)}} \\
&+ \mathbbmss{V}\left[ \left(\frac{A}{\mathbbmss{E}\left[A\right]} - \frac{(1-A)\pi^\star(X)}{(1-\pi^\star(X))\mathbbmss{E}\left[A\right]}\right)(\eta^\star_0(X,0) - \mu^\star_0(X,0)) \right]\,.
\end{split}
\end{equation}
We now turn to the construction of estimators that can asymptotically match the semiparametric efficiency bound derived above.
Estimation and inference
The nuisance functions $\mu^\star$, $\pi^\star$, $\gamma^\star$, and $\eta^\star$ are unknown, and must be estimated from the data at hand. We employ cross-fitting to avoid restrictive Donsker conditions and to retain full-sample efficiency bickel1988estimating, chernozhukov2018double, robins2008higher, schick1986asymptotically,zheng2010asymptotic. Cross-fitting works as follows. We first randomly split the observations $\{\mathcal{D}_1,\dots,\mathcal{D}_{n}\}$ into $J$ disjoint folds (without loss of generality, we assume that the number of observations $n$ is divisible by $J$). For each $j=1,\ldots, J$ we form $\hat{\mathbbmss{P}}^{[-j]}$ with all but the $j$-th fold, and $\mathbbmss{P}_{n}^{[j]}$ with the $j$-th fold. Then, we learn $\hat{\mu}^{[-j]}$, $\hat{\pi}^{[-j]}$, $\hat{\gamma}^{[-j]}$, and $\hat{\eta}^{[-j]}$ on $\hat{\mathbbmss{P}}^{[-j]}$, and compute the final estimator $\hat{\theta}$ by solving the estimating equation
equation[equation omitted — 212 chars of source]
For simplicity, we assume that $\hat{\eta}^{[-j]}$ is estimated on an independent subsample of $\hat{\mathbbmss{P}}^{[-j]}$, that is, it is independent from $\hat{\mu}^{[-j]}$, $\hat{\pi}^{[-j]}$, $\hat{\gamma}^{[-j]}$. The solution to the previous estimating equation can be more conveniently expressed as:
equation[equation omitted — 79 chars of source]
where the form of $\hat{\theta}^{[j]}$ depends on the missing at random assumption. Under Assumption (ref) and Assumption (ref), it is equal to:
equation[equation omitted — 653 chars of source]
while under Assumption (ref) and Assumption (ref), it is equal to:
equation[equation omitted — 1,011 chars of source]
remark[Multiple robustness]
The structure of the estimators in Eq. (ref) and Eq. (ref) sheds light on their multiple robustness property. Under Assumption (ref), we need either the model for $\mu^\star$ or both the models for $\pi^\star$ and $\gamma^\star$ to be well-specified in order to achieve consistency. Under Assumption (ref), we also need the nested regression to be consistent if the propensity score is misspecified. See Appendix Table (ref) for a description of the possible combinations of nuisance functions required for consistency.
Estimating the nested regression
The nested regression function $\eta^\star_0(x,0) = \mathbbmss{E}\left[\mu^\star_0(x,Y_1,0)\mid X=x,A=0\right]$ that appears in our estimator under Assumption (ref) is of particular interest. Here, we showcase a regression-based approach to estimate it, and we defer to Appendix Section (ref) a second approach based on conditional densities. For simplicity, we develop the arguments in this Section by assuming that the dataset is splitted in two folds, the former, $\hat{\mathbbmss{P}}$, being employed for nuisance training and the latter, $\mathbbmss{P}_n$, for influence function averaging.
In the regression-based approach, we treat $\mu^\star_0(x,Y_1,0)$ as the response variable with $X$ and as predictor, and we fit a model to learn $\eta^\star_0(x,0)$. Of course, $\mu^\star_0(x,Y_1,0)$ is not known and thus has to be estimated on $\hat{\mathbbmss{P}}$. Interestingly, this nested approach shares many commonalities with the DR-Learner, a method commonly used to estimate heterogeneous treatment effects foster2023orthogonal,kennedy2023towards, and counterfactual regression yang2023forster. First, we introduce some additional notation. Let $\hat{\eta}_0(x,0) = \hat{\mathbbmss{E}}_n\left[\hat\mu_0(x,Y_1,0)\mid X=x,A=0\right]$ be the regression of $\hat\mu_0(x,Y_1,0)$ on the covariates in the averaging sample $\mathbbmss{P}_n$, and $\Tilde{\eta}_0(x,0) = \hat{\mathbbmss{E}}_n\left[\mu^\star_0(x,Y_1,0)\mid X=x,A=0\right]$ be the corresponding oracle estimator regressing the true $\mu^\star_0(x,Y_1,0)$ onto the covariates. Finally, denote the oracle risk $\Tilde{\Delta}_n^2(x) = \mathbbmss{E}\left[(\Tilde{\eta}_0(x,0) - \eta^\star_0(x,0))^2\right]$ and the conditional bias as $\hat b(x,y_1,0) = \mathbbmss{E}\left[\hat\mu_0(X,Y_1,0) -\mu^\star_0(X,Y_1,0)\mid \hat{\mathbbmss{P}}, X=x, Y_1=y_1\right]$. We can now introduce the definition of stable estimator.
definition[Stability of estimator]
The regression estimator $\hat{\mathbbmss{E}}_n\left[\cdot\mid X=x,A=0\right]$ is defined as stable at $X = x$ and $A=0$ (with respect to a distance metric $d$) if
\begin{equation}
\frac{\hat\eta_0(X,0) - \Tilde{\eta}_0(X,0) - \hat{\mathbbmss{E}}_n\left[\hat b(X,Y_1,0) \mid X, A=0\right]}{\Tilde{\Delta}_n(X)} \overset{p}{\rightarrow} 0\,,
\end{equation}
as $d(\hat\mu_0(X,Y_1,0), \mu^\star_0(X,Y_1,0))\overset{p}{\rightarrow} 0$.
remarkThis pointwise (or local) definition, introduced by kennedy2023towards and extended to $\mathbbmss{L}_2$-norm by rambachan2022robust, is inspired by empirical processes and can be rewritten as
\begin{equation}
\hat\eta_0(X,0) - \Tilde{\eta}_0(X,0) = \hat{\mathbbmss{E}}_n\left[\hat b(X,Y_1,0) \mid X, A=0\right] + o_{\mathbbmss{P}^\star}\left(\Tilde{\Delta}_n(X) \right)\,,
\end{equation}
as $d(\hat\mu_0(X,Y_1,0), \mu^\star_0(X,Y_1,0))\overset{p}{\rightarrow} 0$.
When estimating the nested regression function $\eta^\star_0(x,0) = \mathbbmss{E}\left[\mu^\star_0(x,Y_1,0)\mid X=x,A=0\right]$, a misspecified model for $\mu^\star_0(x,y_1,0)$ would imply a misspecified model for $\eta^\star_0(x,0)$. In fact, assuming that the regression estimator $\hat{\mathbbmss{E}}_n\left[\cdot\mid X=x,A=0\right]$ is stable, the bias term shows first-order dependence on the estimation error in the nuisance function. To gain robustness, one can instead learn an estimate of the augmented nested regression $\mathbbmss{E}\left[\mu^\star_0(x,Y_1,0) + R_0 \left(Y_0 - \mu^\star_0(x,Y_1,0)\right) / \gamma^\star(x,Y_1,0) \mid X=x,A=0\right]$. This strategy provides protection against misspecification in $\mu^\star_0(x,y_1,0)$, provided that a good estimate $\gamma^\star(x,y_1,0)$ is available. We can make this property more formal.
theorem[Oracle property]
Assume that the regression estimator $\hat{\mathbbmss{E}}_n\left[\cdot\mid X=x,A=0\right]$ is stable with respect to distance $d$. Assume also that the estimated pseudo-outcomes converge in probability to truth, i.e.
\begin{equation}
d(\hat\mu_0(x,Y_1,0) + R_0 \left(Y_0 - \hat\mu_0(x,Y_1,0)\right) / \hat\gamma(x,Y_1,0), \mu^\star_0(x,Y_1,0) + R_0 \left(Y_0 - \mu^\star_0(x,Y_1,0)\right) / \gamma^\star(x,Y_1,0))\overset{p}{\rightarrow} 0\,.
\end{equation}
Then
\begin{equation}
\hat\eta_0(X,0) - \Tilde{\eta}_0(X,0) = \hat{\mathbbmss{E}}_n\left[\hat b(X,Y_1,0) \mid X=x,A=0\right] + o_{\mathbbmss{P}^\star}\left(\Tilde{\Delta}_n(X) \right)\,,
\end{equation}
where
\begin{equation}
\hat b(x,y_1,0) = \frac{(\hat\mu_0(x,y_1,0) - \mu^\star_0(x,y_1,0) ) (\hat\gamma(x,y_1,0) - \gamma^\star(x,y_1,0))}{\hat\gamma(x,y_1,0)}\,,
\end{equation}
and $\hat\eta_0(X,0)$ is oracle efficient if $\hat{\mathbbmss{E}}_n\left[\hat b(X,Y_1,0) \mid X=x,A=0\right] = o_{\mathbbmss{P}^\star}\left(\Tilde{\Delta}_n(x) \right)$.
remarkThe previous Theorem gives conditions for achieving the oracle rate, which can be phrased in terms of nuisance smoothness and problem dimension. For example, if we assume that $\mu^\star_0(X,Y_1,0)$ is $\alpha$-smooth in $X$ and $\beta$-smooth in $Y_1$, and $\gamma^\star(X,Y_1,0)$ is $\rho$-smooth, then we can build nonparametric estimators $\hat\mu_0(X,Y_1,0)$ and $\hat\gamma(X,Y_1,0)$ whose convergence rates are $n^{-1/(2+p/\alpha + 1/\beta)}$ and $n^{-1/(2+(p+1)/\rho)}$, respectively. By the previous Theorem, we can achieve the oracle rate whenever $\rho \geq (p+1) / [\alpha(2+1/\beta)(2+p/\alpha+1/\beta) / p - 2]$. This is relevant because, without augmentation, we would typically only achieve the oracle rate when $p/\alpha\to0$, which holds only when $\alpha\to\infty$. In contrast, with augmentation, for each finite $\alpha$ we can find a suitable $\rho$ that satisfies the above oracle rate condition.
Inference
We can now introduce a set of high-level assumptions that are useful for inferential purposes.
assumption[Inference]
Let the number of cross-fitting folds be fixed at $J$, and assume that:
\begin{enumerate}[label=\alph*.]
• For each $j\in\{1,\dots,J\}$, one has \begin{equation}
\varphi\left(\mathcal{D};\theta^\star;\hat{\mu}^{[-j]};\hat{\pi}^{[-j]};\hat{\gamma}^{[-j]};\hat{\eta}^{[-j]}\right)\overset{\mathbbmss{L}^{2}}{\rightarrow} \varphi\left(\mathcal{D};\theta^\star;\mu^\star;\pi^\star;\gamma^\star;\eta^\star\right)\,.
\end{equation}
• For each $j\in\{1,\dots,J\}$, one has
\begin{equation}
\sqrt{n} \sum_{j=1}^J \kappa^{[j]} = o_\mathbbmss{P}(1)\,,
\end{equation}
where the remainder term $\kappa^{[j]}$ is defined as $\kappa^{[j]} = \theta\left(\hat{\mathbbmss{P}}^{[-j]}\right) - \theta(\mathbbmss{P}^\star) + \mathbbmss{E}\left[\varphi\left(\mathcal{D};\theta^\star;\hat{\mu}^{[-j]};\hat{\pi}^{[-j]};\hat{\gamma}^{[-j]};\hat{\eta}^{[-j]}\right)\right]$.
• Given $\xi>0$, for each $j\in\{1,\dots,J\}$, $\hat{\pi}^{[-j]}$ and $\hat{\gamma}^{[-j]}$ are bounded away from $\xi$ and $1-\xi$ with probability 1.
\end{enumerate}
remarkThese are all standard assumptions routinely employed in causal inference. See kennedy2024semiparametric for a review. Assumption a is used to control the empirical process term; Assumption b is used to control the remainder term; Assumption c is commonly referred to as strong overlap.
theorem[Asymptotic normality and efficiency]
Under Assumption (ref), the estimator $\hat\theta$ obtained by solving the estimating equation (ref) is asymptotically normal and achieves the semiparametric efficiency bound as $n$ goes to infinity, that is
\begin{equation}
\sqrt{n} \left( \hat\theta - \theta^\star \right) \leadsto \mathcal{N}\left(0,\mathbbmss{V}\left[\varphi\left(\mathcal{D}\right)\right]\right)\,.
\end{equation}
Endowed with a procedure for valid asymptotic inference, we now turn to the empirical validation of our estimators.
Simulation study
To evaluate the finite-sample performance of our proposed estimators and to validate their theoretical properties, we conduct an extensive Monte Carlo simulation study. The primary goal is to assess the estimators' multiple robustness, particularly when some or all of the nuisance function models are misspecified. Our data generating process (DGP), which closely resembles the one pioneered by kang2007demystifying and then revisted by sant2020doubly, is designed to create scenarios where each model can be independently well or misspecified. For each simulation run, we generate a data set with a sample size of $n=2000$ (Appendix Section (ref) provides simulations for additional sample sizes). The true ATT is fixed at $\theta^\star=5$. We repeat our simulations for 500 runs. The DGP is structured as follows:
itemize• Covariates. We first generate a set of four true baseline covariates $Z=(Z_1,Z_2,Z_3,Z_4)$ from a standard normal distribution. From these, we create a corresponding set of four observed covariates $X=(X_1,X_2,X_3,X_4)$ by applying complex, non-linear transformations to $Z$ (e.g., involving exponential, polynomial, and interaction terms). This setup ensures that models based on $Z$ are correctly specified, while models using only $X$ are misspecified. In particular, the $X$'s are defined as
\begin{equation}
X_{1} = \exp\left(\frac{Z_{1}}{2}\right)\,, \quad
X_{2} = \frac{Z_{2}}{1+\exp(Z_{1})} + 10\,,\quad
X_{3} = \left(\frac{Z_{1}Z_{3}}{25} + 0.6\right)^3\,, \quad
X_{4} = (Z_{2} + Z_{4} + 20)^2\,.
\end{equation}
• Treatment. The treatment assignment indicator, $A$, is generated from a Bernoulli distribution with a probability (propensity score) that is a logistic function of the true covariates $Z$, that is
\begin{equation}
logit(\pi(Z)) = logit\left(\mathbbmss{P}\left[A=1 \mid Z\right]\right) = -Z_{1} + 0.5Z_{2} - 0.25Z_{3} - 0.1Z_{4}\,.
\end{equation}
• Outcomes. Pre-treatment outcomes are generated as linear functions of $Z$ plus standard normal error terms, that is
\begin{equation}
Y_0 = 210 + 27.4Z_{1} + 13.7Z_{2} + 13.7Z_{3} + 13.7Z_{4} + \varepsilon\,, \quad where \quad \varepsilon \sim \mathcal{N}\left(0,1\right)\,.
\end{equation}
Post-treatment potential outcomes are then generated as
\begin{equation}
\begin{split}
&Y_1^{(0)} = Y_0 + \varepsilon\,, \quad \text{where} \quad \varepsilon \sim \mathcal{N}\left(0,1\right)\,, \\
&Y_1^{(1)} = Y_1^{(0)} + A\theta^\star\,,
\end{split}
\end{equation}
which are then summarized, by consistency assumption, into the observed post-treatment outcome $Y_1 = A Y_1^{(1)} + (1-A) Y_1^{(0)}$.
• \textbf{Missingness mechanism.} To align with our theoretical framework, we simulate two distinct missingness patterns for the pre-treatment outcome $Y_0$. Under Assumption (ref), the missingness indicator $R_0$ is drawn from a Bernoulli distribution where the probability is a logistic function of the true covariates $Z$ and the treatment status $A$, i.e.
\begin{equation}
\text{logit}(\textcolor{black}{\gamma}(Z,A)) = \text{logit}\left(\mathbbmss{P}\left[R_0=1 \mid Z,A\right]\right) = - 0.25 Z_{1} - 0.1 Z_{2} - 0.5 Z_{3} + 0.3 Z_{4} - 0.2 A\,.
\end{equation}
Under Assumption (ref), this probability additionally depends on the post-treatment outcome $Y_1$, creating a more complex logistic dependency:
\begin{equation}
\text{logit}(\textcolor{black}{\gamma}(Z,A,\textcolor{black}{Y_1})) = \text{logit}\left(\mathbbmss{P}\left[R_0=1 \mid Z,A,\textcolor{black}{Y_1}\right]\right) = - 0.25 Z_{1} - 0.1 Z_{2} - 0.5 Z_{3} + 0.3 Z_{4} - 0.2 A + \textcolor{black}{0.01} Y_1 \,.
\end{equation}
For each of the 500 simulated datasets, we estimate the ATT using the influence function-based estimators in Equations (ref) and (ref). The nuisance functions are estimated using standard logistic and ordinary least squares models. To test the multiple robustness property, we cycle through all possible combinations of correctly specified and misspecified models for the nuisance functions. A model is correctly specified if it uses the true covariates $Z$ as predictors and misspecified if it uses the observed, non-linear covariates $X$.
We evaluate the performance of our estimators across these different specifications using two standard metrics: bias, which is the difference between estimated ATT $\hat\theta$ and true ATT $\theta^\star$; and root mean squared error (RMSE), which is the square root of the average squared difference between the estimate and the true value. RMSE penalizes large errors and captures both bias and variance, providing a comprehensive measure of estimator quality. \textcolor{black}{In Appendix Section (ref), we also display simulation results for a third evaluation metric, empirical coverage.}
figure[figure omitted — 905 chars of source]
The results, displayed in Figure (ref) and Appendix Tables (ref) and (ref), provide strong evidence for the theoretical properties of our estimators. As predicted by theory, the bias and the RMSE are negligible for both estimators in all scenarios where the conditions for consistency are met. This holds true, for example, when the outcome model is correctly specified. Conversely, the estimators show a clear bias and higher RMSE in the theoretically inconsistent scenarios. This demonstrates the estimators' breaking point\textcolor{black}{, which we further investigate with additional simulations in Appendix Section (ref)}. Overall, the simulation results strongly support the validity and robustness of the proposed estimators.
Conclusions
The Difference-in-Differences (DiD) method is a cornerstone of applied research, yet its validity is often threatened by the practical challenge of missing outcome data -- a problem that can introduce significant selection bias and invalidate standard estimators. This paper addresses this critical gap by developing a rigorous and comprehensive framework for DiD estimation when pre or post-treatment outcomes are missing at random (MAR). Drawing on semiparametric theory, we make several key contributions. First, we establish nonparametric identification of the Average Treatment Effect on the Treated (ATT) under two distinct and plausible MAR mechanisms: one where missingness is independent of the outcome conditional on covariates, and another where it may depend on the post-treatment outcome. For each setting, we derive the semiparametric efficiency bound, establishing a formal benchmark for asymptotic precision. We then propose novel estimators that achieve these bounds, ensuring asymptotic semiparametric efficiency. A critical feature of our estimators is their multiple robustness, which guarantees consistency as long as a subset of the nuisance function models is correctly specified, providing a layer of protection against model misspecification in practice.
The implications of this work are both theoretical and practical. Our framework provides applied researchers in economics, public health, and social sciences with a principled and efficient toolkit to conduct credible DiD analysis using incomplete panel data. By formally accounting for missing data, our estimators enhance the reliability of causal claims drawn from real-world observational studies where complete data is the exception rather than the rule.
While this paper focuses on the canonical two-group, two-period setting for clarity, the principles developed here open several avenues for future research. A natural next step is the extension of this framework to more complex scenarios, such as the staggered treatment adoption settings that have been the focus of much recent literature. Further investigation into the performance of different machine learning methods for the nuisance components, particularly the nested regression function, would also be valuable. Finally, incorporating our efficient estimators into standard DiD software would greatly facilitate their adoption by the broader research community, strengthening the quality and credibility of causal inference across disciplines.