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.
62,745 characters · 13 sections · 27 citation commands
Difference-in-Differences with a Mediator
In the economic and health sciences, studying the mechanisms by which treatment affects outcomes has been both interesting and challenging. The treatment can exert its effects either directly on the outcome or indirectly through a mediator. Mediation analysis was initially proposed within linear models baron1986moderator. Later, mediation analysis was formalized using potential outcomes, in which the mediator is a potential outcome of the treatment and the primary outcome is a potential outcome of the treatment and the mediator rubin2004direct, goetgeluk2008estimation. By controlling the mediator at the natural level, the total effect is decomposed into a natural indirect and a natural direct effect. Identifiability assumptions have been proposed to identify the natural indirect and direct effects imai2010general, imai2010identification. In general, mediation analysis methods require that the treatment assignment, mediator, and primary outcome be unconfounded, which is a stringent condition in observational studies.
As a quasi-experimental design, difference-in-differences (DiD) was proposed to estimate treatment effects when the treatment is confounded with the outcome using panel data heckman1997matching, abadie2005semiparametric. By leveraging pre-treatment outcomes, the change in potential outcomes from pre- to post-treatment is assumed to follow a common time trend across treatment groups. The pre-treatment outcome serves as a negative control variable for the unmeasured confounding zivich2023introducing. Although the average treatment effect on the entire population is not identifiable due to unmeasured confounding, the mean potential outcome under control can be imputed for the treated units. Under appropriate assumptions, the average treatment effect in the treated group is identifiable. Doubly robust and locally efficient estimators have been proposed to estimate this treatment effect by deriving the efficient influence function chang2020double, sant2020doubly, callaway2021difference, deng2025improved.
Conducting mediation analysis in a difference-in-differences design can help address non-randomized treatment assignment. Existing work has followed several distinct identification strategies. deuchert2019direct makes a decomposition of the total effect into natural indirect and direct effects in experiments with non-compliance, assuming a monotonic effect of treatment on the binary mediator and imposing a parallel trends assumption within principal strata. huber2022direct propose an alternative identification strategy based on a changes-in-changes framework that avoids selection-on-observables assumptions and bypasses reliance on instrumental variables, relying instead on distributional restrictions for binary mediators and continuous outcomes. In contrast, hsia2025causal extended the causal mediation analysis to panel data using linear regression models, in which case the natural indirect effect can be inferred from the product of coefficients. However, the linear model imposes strong restrictions and is not flexible enough to accommodate interactions. Complementary to these works on natural effects, blackwell2025estimating considered efficient estimation of controlled indirect and direct effects when the treatment is randomized with a discrete mediator. They relied on the assumption of no modified direct effect for identification. Another limitation is that their estimation cannot be applied if the mediator is continuous. The controlled effect does not yield a heuristic decomposition of the total effect, especially when the mediator is multi-level or even continuous. huber2026difference identified the treatment effect conditional on an observed mediator in the treated group, and defined natural effects by integrating over the distribution of the mediator. However, their estimand can only be interpreted as a conditional treatment effect rather than a controlled effect that intervenes in the mediator if the mediator is not randomized.
Motivated by the Job Corps Study, we aim to study the effect of the job training program on earnings schochet2001national. Although the assignment was randomized, the actual treatment was not, because participants could decide whether to join training programs regardless of their treatment assignment. Some previous work used instrumental variables or principal stratification to account for non-compliance; other work constructed bounds for the treatment effect schochet2008does, lee2009training, zhang2009likelihood, flores2010learning, flores2012estimating, blanco2013bounds, chen2015bounds. These analyses indicated that job training has a significant long-term effect but an insignificant short-term effect on earnings. This is possibly because training takes up a large proportion of time, so the wage is lower than that of full-time employees. Therefore, it is essential to quantify the treatment effect mediated by work time while accounting for non-randomized treatment.
In this paper, we consider the identification and estimation of mediated effects in two-period difference-in-differences. Unlike the existing mediation literature in a difference-in-differences design, our strategy is built within a standard DiD framework and allows for unmeasured confounding between treatment and outcomes. Under certain assumptions, we show that the natural direct and indirect effects are identifiable in the treated group. These identifiability assumptions essentially require that there be no unmeasured confounding between the mediator and the treatment, and between the mediator and the change in potential outcomes. The presence of unmeasured confounding between treatment and outcomes is allowed. These assumptions are reasonable in the context of the Job Corps Study. Using semiparametric theory, we derive efficient influence functions for the natural and total effects. We then construct estimators that are multiply robust and locally nonparametrically efficient. The estimators are consistent if either two of the propensity score, mediator distribution, and the outcome regression models are correctly specified. We further consider estimation and inference for the controlled effect. To address the challenge of pathwise differentiability failure when the mediator is continuous, we use kernel smoothing to approximate the efficient influence functions. The estimator for controlled effects has a lower convergence rate but maintains multiple robustness. By analyzing data from the Job Corps Study, we find that participation in the training program significantly increases short-term earnings, after controlling for the proportion of weeks employed.
The remainder of this paper is organized as follows. Section (ref) presents the mediation DiD framework, defines the natural indirect and direct effects, and establishes their identifiability. Section (ref) develops semiparametric estimators based on efficient influence functions. A practical estimation strategy is provided by adopting generalized linear models as working models. Section (ref) extends the estimation and inference to controlled effects. Section (ref) investigates the performance of the proposed estimators through simulation studies. Section (ref) analyzes data from the Job Corps Study. Section (ref) concludes with a discussion.
In a two-group and two-period setting, let $G \in \{0,1\}$ be the treatment group assignment and $t \in \{0,1\}$ be the period indicator. Let $D_t$ be the treatment received in period $t$. In the pre-treatment period, all units are unexposed to treatment, $D_0=0$; in the post-treatment period, the units in the treated group received active treatment, $D_1=G$. Therefore, $D_t = Gt$; that is, $D_t=1$ only if $G=1$ and $t=1$. Let $Y_0(g)$ denote the potential outcome in the pre-treatment period under treatment assignment $g$. Suppose there is a mediator that may mediate the treatment effect on the outcome. Let $M(g)$ denote the potential mediator under treatment assignment $g$. The mediator can be either discrete or continuous. Let $Y_1(g,m)$ denote the potential outcome in the post-treatment period under treatment assignment $g$ given the mediator $m$. When the mediator is at the natural value $M(g)$ under treatment assignment $g$, the potential outcome in the post-treatment period is $Y_1(g,M(g))$, which we denote as $Y_1(g)$. Sometimes the mediator can be measured in the pre-treatment period, denoted as $M_0(g)$. This pre-treatment mediator should not be affected by treatment assignment, as the unit has not been exposed to the treatment. Therefore, it is reasonable to assume no anticipation for the mediator $M_0(0)=M_0(1)$ and simply treat $M_0(g)$ as a baseline covariate.
A key question in mediation analysis is the extent to which the treatment effect is mediated by the mediator. In our setting, treatment may not be randomized, so the distribution of unmeasured features may differ across treatment groups. Consistent with the difference-in-differences literature, we adopt the treated group $G=1$ as the target population to minimize the identification requirements. The total effect is
By switching the treatment associated with the mediator, we define the natural indirect effect (NIE) as \[ \tau_{IE} = E\{Y_1(0,M(1))-Y_1(0,M(0)) \mid G=1\}. \] By switching the treatment associated with the outcome, we define the natural direct effect (NDE) as \[ \tau_{DE} = E\{Y_1(1,M(1))-Y_1(0,M(1)) \mid G=1\}. \] The natural indirect effect measures the effect delivered through the mediator, and the natural direct effect measures the effect delivered without being mediated by the mediator. Accordingly, the total effect can be decomposed into the natural indirect and direct effects.
Let $X$ be baseline covariates measured prior to treatment, which may include pre-treatment mediator. Let $Y_0$ and $Y_1$ be the observed outcomes in the pre-treatment and post-treatment periods. The observed data consist of $n$ independent and identically distributed copies of $O=(G,X,Y_0,M,Y_1)$. To identify the natural indirect, direct, and total effects for the treated group, we adopt a difference-in-differences framework and impose the following assumptions.
Assumption (ref) rules out anticipation effects, requiring that future treatment assignment does not affect the potential outcome in the pre-treatment period. In other words, the potential outcome is a function of actual treatment, $Y_0(g) = Y_0(d)$ and $Y_1(g,m) = Y_1(d,m)$. Assumption (ref) extends the standard parallel trends assumption in difference-in-differences to accommodate a mediator, requiring that, conditional on baseline covariates and the effect on the mediator, the time trend in untreated potential outcomes is identical across groups. Although treatment assignment may be confounded with potential outcomes, taking differences over time removes time-invariant unobserved confounding. Compared with other literature on mediational difference-in-differences blackwell2025estimating, huber2026difference, we do not condition on the observed post-treatment variable $M$ in the parallel trends assumption, thereby avoiding potential collider bias.
Sequential ignorability is a common assumption in causal mediation analysis, excluding pairwise unmeasured confounding between the treatment, mediator, and outcome. The first part of sequential ignorability requires that there be no unmeasured confounding between treatment assignment and the mediator. The second part of sequential ignorability requires that there be no unmeasured confounding between the mediator $M(1)$ and the change in potential outcomes, so that the change in potential outcomes is independent of the observed mediator in the treated group. In other words, the observed mediator should not modify the underlying parallel trends. This assumption is weaker than the classical sequential ignorability assumption in mediation analysis because it does not require the mediator $M(1)$ to be randomized, i.e., conditionally independent of $Y_1(1,m)$, in the treated group.
To provide intuition for these identification assumptions, we consider the following structural causal model (SCM):
where $\epsilon_G$, $\epsilon_0$, $\epsilon_M$, and $\epsilon_1$ are external random errors. The unmeasured confounder $U$ influences $G$, $Y_0$, and $Y_1$, but not $M$, implying the first part of sequential ignorability. The change in potential outcomes under control $Y_1(0,m)-Y_0(0) = f_1(X,0,m)-f_0(X)+\epsilon_1-\epsilon_0$ does not depend on $G$ (as a function of $U$ and $\epsilon_G$) and $M(0)$ (as a function of $\epsilon_M$) conditional on $X$, implying the parallel trends. The change in potential outcomes under control and mediator $m$, $Y_1(0,m)-Y_0(0) = f_1(X,0,m)-f_0(X)+\epsilon_1-\epsilon_0$ does not depend on $M(g)$ (as a function of $\epsilon_M$) conditional on $X$ and $G$, implying the second part of sequential ignorability. Figure (ref) shows the directed acyclic graph (DAG) and single-world intervention graph (SWIG) of the mediator and primary outcomes. Due to the unmeasured confounding $U$ between $G$ and $Y_0$, the distribution of $Y_0$ is different between groups. However, since $U$ does not confound $M$ and $(Y_0,Y_1)$, the change in potential outcomes $Y_1(g,m)-Y_0(0)$ is independent of $M(g)$ and the cross-world $M(1-g)$. The mediator opens a “front door” between $G$ and $Y_1$, unaffected by unmeasured confounding pearl1995causal.
Another structural causal model that satisfies these assumptions is:
The treatment assignment $G$ is randomized, but the mediator is influenced by the unmeasured confounder $U$. The change in potential outcomes under control, $Y_1(0,m)-Y_0(0) = f_1(X,0,m)-f_0(X)+\epsilon_1-\epsilon_0$, does not depend on $G$ (as a function of $\epsilon_G$) and $M(0)$ (as a function of $U$ and $\epsilon_M$) conditional on $X$, implying the parallel trends. The potential mediator $M(g) = m(X,g,U)+\epsilon_M$ is independent of $G$ conditional on $X$, implying the first part of sequential ignorability. The change in potential outcomes under control and mediator $m$, $Y_1(0,m)-Y_0(0) = f_1(X,0,m)-f_0(X)+\epsilon_1-\epsilon_0$ does not depend on $M(g)$ (as a function of $U$ and $\epsilon_M$) conditional on $X$ and $G$, implying the second part of sequential ignorability.
Positivity ensures that there are individuals in both groups. The distribution of covariates and mediators should overlap between groups. Consistency links the potential outcomes and observed data.
To facilitate identification, we define $\delta(g,m,x) = E(Y_1-Y_0 \mid G=g,M=m,X=x)$, which characterizes the conditional mean of the change in observed outcomes between periods. The identification problem then reduces to expressing $\tau_{IE}$, $\tau_{DE}$, and $\tau$ in terms of the observable function $\delta(g,m,x)$. The following theorem formalizes this representation and provides an explicit identification of the natural indirect, direct, and total effects.
There are two expectations for $\delta(g,M,X)$, a function of random variables $M$ and $X$, in the expression of $\tau(g,g^*)$. The argument $g$ corresponds to the level of treatment, and the argument $g^*$ corresponds to the level of mediator. The inner expectation integrates $M$ out, transforming the mediator $M$ from the distribution in the group $g$ to the group $g^*$, giving a function of $X$. The outer expectation integrates $X$ out, transforming the covariates $X$ from the distribution in group $g^*$ to the treated group $G=1$. Among these quantities, $\tau(1,1)$ is directly identifiable from the treated group, since $\tau(1,1) = E\{\delta(1,M,X) \mid G=1\}$, whereas $\tau(0,0)$ and $\tau(0,1)$ require transporting both mediator and covariate distributions across groups, posing additional challenges for estimation. To represent this shift, we define $\nu(g^*,X) = E\{\delta(0,M,X) \mid G=g^*,X\}$.
Motivated by the identifiability result, we observe that $\tau(g,g^*)$ can be estimated by imputation. Specifically, we first regress $\Delta Y := Y_1-Y_0$ on $(G,M,X)$ and generate predicted pseudo-outcomes $\widehat\delta(0,M,X)$ based on the regression for each unit given observed $(M,X)$. Suppose that we estimate the conditional density of the mediator given $G=g^*$ and $X=x$ by $\widehat{f}(m \mid g^*,x)$. Next, we estimate $E\{\delta(0,M,X)\mid G=g^*,X\}$ by $\widehat\nu(g^*,X) = \int \widehat\delta(0,m,X)\widehat{f}(m\mid G=g^*,X)dm$. Finally, the empirical average of $\widehat\nu(g^*,X)$ in the treated group is an estimate for $\tau(0,g^*)$. However, the imputation-based estimator is prone to model misspecification. If either one of the models $\delta(0,m,x)$ and $f(m\mid 0,x)$ is misspecified, the imputation-based estimator is not consistent. In addition, the imputation-based estimator is inefficient because it does not fully utilize all the information available in the observed data.
To improve robustness and efficiency, we derive the efficient influence function (EIF) of $\tau(g,g^*)$, where $(g,g^*) \in \{(1,1), (0,0), (0,1)\}$. To achieve this, we need to find the tangent space, the collection of all infinitesimal directions in which the data-generating distribution can be perturbed while remaining inside the statistical model. The tangent space of the model space is shown in Lemma 1 of the Supplementary Material. For the estimand $\tau(g,g^*)$, there exists a unique regular and asymptotically linear (RAL) estimator whose influence function is in the tangent space bickel1993efficient. The influence function in the tangent space is called the efficient influence function (EIF). Deriving the efficient influence function enables the construction of estimators that are both robust to nuisance model misspecification and asymptotically efficient. Let $\pi(x) = P(G=1 \mid X=x)$ be the propensity score, $f(m \mid g,x)$ be the conditional density (or probability) of the mediator, and $\delta(g,m,x)$ be the conditional mean of the change in observed outcomes. The following theorem provides the closed-form expressions of the efficient influence function for $\tau(g,g^*)$.
The efficient influence function provides a lower bound for the variance of regular and asymptotically linear estimators. To show $\varphi(g,g^*)$ is the EIF of $\tau(g,g^*)$, it suffices to show that $\varphi(g,g^*)$ is in the tangent space $\dot{\mathcal{T}}$ and $\varphi(g,g^*)$ is a valid influence function for $\tau(g,g^*)$. The proof is given in the Supplementary Material.
The efficient influence functions involve three models: propensity score, conditional density of mediator, and outcome change model. In practice, we evaluate the influence functions at estimated nuisance functions $\widehat\pi(x)$, $\widehat{f}(m\mid g,x)$, and $\widehat\delta(g,m,x)$. However, modeling the conditional density is computationally infeasible because it has too much freedom. A poor estimation of $f(m\mid g,x)$ renders the fitting of $\varphi(0,0)$ and $\varphi(0,1)$ inaccurate. First, the imputation-based method estimates $\nu(0,X) = E\{\delta(0,M,X) \mid G=0,X\}$ by $\widehat\nu(0,X) = \int \widehat\delta(0,m,X)\widehat{f}(m\mid 0,X)dm$. Second, $f(m\mid g,x)$ directly appears in $\varphi(0,1)$. To avoid the curse of dimensionality in estimating $f(m \mid g,x)$, we note that by Bayes' formula,
so we can summarize the ratio of propensity scores and conditional mediator densities into a pseudo-propensity score $\varpi(m,x) = P(G=1\mid M=m, X=x)$. The response variable $G$ is binary, so modeling $\varpi(m,x)$ is much easier than modeling the conditional density. For example, it can be fitted by logistic regression, denoted by $\widehat\varpi(m,x)$. The pseudo-propensity score is well defined regardless of whether $M$ is continuous or discrete. Even if $M$ follows a mixture distribution, Equation (ref) is valid provided that the Radon-Nikodym derivative of the conditional distribution function of $M$ in $G=1$ with respect to $G=0$ exists. In addition, $\nu(0,x)$ can be estimated by imputing $m$ by the estimated conditional mean $E(M \mid G=0,X=x)$ in the regression model $\delta(0,m,x)$, respecting the dependence structure of $Y$ on both $M$ and $X$ without involving $f(m\mid 0,x)$. In this way, we denote the estimate of $\nu(0,x)$ by $\widehat\nu(0,x)$.
For any measurable function of data $q(o)$, let $\mathbb{P}(q) = \int q(o)f(o)do$ be the measure with respect to the true data-generating mechanism and $\mathbb{P}_n(q) = n^{-1}\sum_{i=1}^{n} q(O_i)$ be the empirical measure. We denote the $L_1$-norm of $q$ as $\|q\|_{L_1} = \mathbb{P}(|q|)$ and the $L_2$-norm as $\|q\|_{L_2} = \{\mathbb{P}(q^2)\}^{1/2}$. In practice, $\tau(g,g^*)$ is estimated by solving the empirical estimating equation $\mathbb{P}_n\varphi(g,g^*) = 0$, which leads to
The natural indirect, direct, and total effects are then estimated as \[ \widehat\tau_{IE} = \widehat\tau(0,1)-\widehat\tau(0,0), \quad \widehat\tau_{DE} = \widehat\tau(1,1)-\widehat\tau(0,1), \quad \widehat\tau = \widehat\tau(1,1)-\widehat\tau(0,0). \] Next, we establish a general consistency result for the proposed estimators.
The proposed estimator relies on four models: $\pi(x)$, $\varpi(m,x)$, $\delta(0,m,x)$, and $\nu(0,x)$. The Glivenko--Cantelli condition requires that the working models are not too complex vaart2023empirical. Common models satisfy this condition, including (generalized) linear models and kernel regression models. Note that $\nu(0,x) = \int_m \delta(0,m,x)f(m\mid 0,x) dm$ is a convex hull of $\delta(0,m,x)$ with weight $f(m\mid 0,x)$ and $\varpi(m,x)$ is a function of $\pi(x)$ and $f(m \mid g,x)$, so at least three working models should be specified: $\pi(x)$, $f(m\mid g,x)$, and $\delta(0,m,x)$. With these three working models, this theorem indicates multiple robustness: if either two of $\pi(x)$, $f(m\mid g,x)$, and $\delta(0,m,x)$ are correctly specified, then $\tau_{IE}$, $\tau_{DE}$, and $\tau$ are all consistent in probability. In practice, $\varpi(m,x)$ and $\nu(0,x)$ can be modeled separately. In this case, if either $\{\pi(x),\delta(0,m,x)\}$, $\{\pi(x),\varpi(m,x)\}$, $\{\nu(0,x),\delta(0,m,x)\}$, or $\{\varpi(m,x),\nu(0,x)\}$ are correctly specified, then $\tau_{IE}$, $\tau_{DE}$, and $\tau$ are all consistent in probability.
The Donsker condition requires that the working models are not too complex. Standard parametric and semiparametric models are Donsker. The fitted models are allowed to converge at rates lower than the parametric rate (for example, when using kernel smoothing as the working model). Nevertheless, the estimators for $\tau(g,g^*)$ are regular and asymptotically linear. The influence function in the linear expansion of $\widehat\tau(g,g^*)$ is exactly $\phi(g,g^*)$, the EIF of $\tau(g,g^*)$. Since influence functions are additive, the influence functions of $\widehat\tau_{IE}$, $\widehat\tau_{DE}$, and $\widehat\tau$ are $\varphi_{IE}$, $\varphi_{DE}$, and $\varphi$, respectively. Complex machine learning methods such as random forests and neural networks are not Donsker in general. To remove the Donsker condition, cross-fitting can be applied, where a training set is used to fit the nonparametric model, for example, using random forests or neural networks, and an inference set is used to estimate the treatment effects based on fitted influence functions chernozhukov2018double. A final estimator is obtained by averaging the estimates after switching the training set and inference set. Cross-fitting allows more flexible models; however, when the sample size is small, it may lead to greater finite-sample variation. Therefore, we illustrate our method using simple models satisfying the Donsker condition throughout this paper.
Under the conditions of multiple robustness, if part of the working models is misspecified, say $\widehat\pi(x) \to \pi^*(x)$, $\varpi(m,x) \to \varpi^*(m,x)$, $\widehat\delta(0,m,x) \to \delta^*(0,m,x)$, and $\widehat\nu(0,x) \to \nu^*(x)$ in the $L_2$-norm, then the influence function would consist of two parts: the EIF with the true models replaced with their limiting models, and a remainder term sue to the uncertainty of the estimated models. The asymptotic variance is not necessarily the smallest among regular and asymptotically linear estimators.
In Supplementary Material A, we outline practical estimation and inference procedures using generalized linear models. The standard errors of these estimators can be calculated from the variances of the fitted influence functions. Let $\widehat\varphi(g,g^*)$ be the fitted efficient influence function plugging in fitted models and $\widehat\tau(g,g^*)$. Then the standard errors are given by
The $P$-values are calculated based on the normal approximation.
Studying the marginal treatment effect calls for controlling the mediator at a specific value. We define the controlled direct effect (CDE) \[ \tau_{DE}(m) = E\{Y_1(1,m)-Y_1(0,m) \mid G=1\}. \] The controlled indirect effect reflects the effect of switching treatment when intervening in the mediator at $m$. If we define $\bar\tau(g,m) = E\{Y_1(g,m)-Y_0(0) \mid G=1\}$, then $\tau_{DE}(m) = \bar\tau(1,m)-\bar\tau(0,m)$. Under the identification assumptions, $\bar\tau(g,m)$ is identifiable and has the efficient influence function in the following theorem. Let $\boldsymbol{\delta}(\cdot)$ be the Dirac delta function.
When the mediator is discrete, the corresponding estimator can be constructed by replacing the unknown nuisance functions with their estimates, which we denote as $\widehat\varphi(0,m)$ and $\widehat\varphi(1,m)$. However, if the mediator is continuous, there will be no efficient influence function as in a Dirac mass at $M=m$ to estimate the conditional mean for the pathwise derivative, which does not belong to $L^2(\mathbb{P}_0)$. In such case, we use kernel smoothing to replace the Dirac mass at $m$. Let $K(x)$ be a kernel function and $K_h(x) = h^{-1}K(x/h)$ be a kernel with the bandwidth $h$. We come up with the smoothing efficient influence functions as
The corresponding estimated smoothing efficient influence functions are denoted as $\widehat\varphi_h(g,m)$. Define the corresponding estimator as
solving the estimation equation $\mathbb{P}_n \widehat\varphi(g,m) = 0 $ or $\mathbb{P}_n \widehat\varphi_h(g,m) = 0 $, depending on whether $M$ is discrete or continuous). In Supplementary Material A, we outline practical estimation procedures based on kernel smoothing.
We use $f(m \mid g, x)$ to represent either the conditional probability $P(M=m\mid G=g,X = x)$ or the conditional density. In addition, we assume $f(m \mid g, x)$ and $\delta(g, m, x)$ are at least two-order smooth with respect to $m$. The following theorem states that the above estimation based on Theorem (ref) achieves the near-optimal asymptotic properties regardless of whether $m$ is discrete or continuous, provided that $M$ is scalar.
Note that kernel smoothing breaks the Donsker condition for $\bar\varphi_h(g,m)$, so the convergence rate of $\widehat\tau_h(g,m)$ is smaller than $O_p(n^{-1/2})$. Nevertheless, up to the scale $h$, the influence function for $\widehat\tau_h(g,m)$ over $(g,m)$ can be Donsker. Although Theorem (ref) is stated for the univariate mediator case, the extension to multivariate $M$ is straightforward by appropriately modifying the bandwidth conditions in the continuous case. We refer to the result as near-optimal because, as shown in Theorem (ref), the proposed estimator enjoys double robustness regardless of whether $M$ is discrete or continuous, and it attains the lowest possible asymptotic variance among all regular and asymptotically linear estimators when $M$ is discrete.
In this section, we conduct simulation studies to compare our proposed method with the regression-based method in hsia2025causal, which uses regression coefficients to represent treatment effects as an extension to baron1986moderator.
Let the sample size $n \in \{200, 1000, 5000\}$. We generate two independent baseline covariates $X_i = (X_{1i},X_{2i})^{\top}$, with each element following the standard normal distribution. In addition, we generate an individual-level random effect $u_i \sim N(0,1)$. We consider two settings for generating the mediator. In Setting 1, the potential mediator is continuous, generated as \[ M_i(g) = 0.6X_{1i} - 0.3X_{2i} + g + \varepsilon_i, \] where $\varepsilon_i \sim N(0,1)$. In Setting 2, the potential mediator is binary, generated from \[ P(M_i(g)=1 \mid X_0) = \Phi(0.6X_{1i} - 0.3X_{2i} + g), \] where $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution. The potential outcomes are generated as
where $\epsilon_{0i}$ and $\epsilon_{1i}$ are independent random errors following $N(0,0.5^2)$. Finally, the treatment assignment mechanism is \[ P(G_i=1 \mid X_i) = \frac{\exp(0.3+0.4X_1+0.5X_2)}{1+\exp(0.3+0.4X_1+0.5X_2)}. \] The observed mediator and outcomes $M_i=M_i(G_i)$, $Y_{0i}=Y_{0i}(0)$, and $Y_{1i}=Y_{1i}(G_i,M_i)$.
The propensity score is modeled by logistic regression, and the change in observed outcomes is modeled by linear regression with interaction between $G$, $M$, and $X$, both of which are correctly specified. The conditional density of the mediator is modeled using the transformation technique in Equation (ref). We replicate the data generation 1000 times. The Panel (O) of Table (ref) shows the average bias, average standard error (SE), standard deviation (SD), and coverage percentage (CP) of 95% confidence intervals for the estimated natural indirect, direct, and total effects. We observe that the regression-based method exhibits considerable bias in estimating $\tau_{IE}$ and $\tau$, because the linear model does not account for the interaction between $M$ and $X$. In contrast, the bias of the proposed method is minor. The average standard error and standard deviation are in close agreement. The coverage percentage of the 95% confidence intervals is close to the nominal level.
To demonstrate the potential bias in estimation and inference when models are misspecified, we consider two alternative data-generating mechanisms. First, we assume the potential outcomes are generated as
where $\epsilon_{0i}$ and $\epsilon_{1i}$ are independent random errors following $N(0,0.5^2)$, so that the outcome model is misspecified. Second, we assume the treatment assignment is generated from \[ P(G_{i}=1 \mid X_i) = \Phi(0.3+0.4X_{1i}+0.5X_{2i}+0.3X_{1i}X_{2i}), \] so that the propensity score is misspecified.
The Panels (A) and (B) of Table (ref) show the average bias, average standard error (SE), standard deviation (SD), and coverage percentage (CP) of 95% confidence intervals for the estimated natural indirect, direct, and total effects when either the outcome model or propensity score is misspecified. The regression-based method shows considerable bias because the linear model does not correctly specify the data-generating mechanism. Even if a working model is severely misspecified, the bias of the proposed method remains small, and the coverage percentage of the 95% confidence intervals is near the nominal level. In this case, the proposed method is not efficient.
Finally, we also compare our proposed estimator with regression-based methods for estimating the controlled effects at $m = 0$, considering both continuous and binary mediators under the same setting as above. For the continuous mediator setting, we employ the kernel-weighted local polynomial with degree two to estimate the conditional outcome, utilizing the Gaussian kernel with the bandwidth selected using Silverman's rule of thumb silverman2018density, to localize the estimation around the target mediator. The results are reported in Table (ref).
The U.S. Job Corps experimental study is a job training program for disadvantaged youths schochet2001national. Participants in the treated group are provided access to academic and vocational training. Approximately 60.4% of the 9240 participants were assigned to the treated group. Although the treatment assignment was randomized, not all participants complied with the assignment. In the first year after assignment, only 71.1% actually enrolled in education and/or vocational training. As a result, the realized treatment is no longer randomized. In particular, individuals with more pessimistic expectations about their labor market prospects were more likely to take up the training. The primary outcome of interest is weekly earnings. Previous studies using instrumental variables methods find that job training has a significant effect on earnings three years after assignment, whereas the effect in the first two years is not statistically significant.
Participants in the training program spent a substantial amount of time on training, leaving less time available for work. Consequently, earnings tend to be lower when employees cannot work full-time. Existing studies did not consider this fact and reckon that job training has only long-term effects schochet2008does. In this study, we aim to adjust for the treatment effect on earnings mediated by work time. We let $G$ be the actual enrollment in education and/or vocational training. The pre-treatment outcome $Y_0$ is the average weekly gross earnings at assignment, after $\log(1+x)$ transformation. The post-treatment $Y_1$ is the average weekly earnings in the second year after assignment, after $\log(1+x)$ transformation. The mean non-zero income is 151 at pre-treatment and 175 at post-treatment, much larger than 1; thus, the log transformation is insensitive to the constant. Such a transformation ensures that zero income is always zero before and after the transformation. The mediator $M$ is the proportion of weeks employed during the second year after assignment. We control for four baseline covariates in $X$: sex, age, year of education, and race.
In the control group ($G=0$), 23.7% of the participants did not work in the second year after the assignment, and 19.8% worked full-time. In the treated group ($G=1$), 25.1% of the participants did not work in the second year after the assignment, and 14.8% worked full-time. The average proportion of weeks employed is 49.0% in the control group and 45.0% in the treated group. The average difference in $Y_1$ and $Y_0$ (average ratio of weekly earnings between post-treatment and pre-treatment) is 3.07 in the control group and 3.00 in the treated group. After adjusting for covariates, the average weeks employed is 1.4% lower in the treated group, and the increase in logarithmic earnings is 4.5% larger in the treated group, but the effects are not significant.
For each individual, baseline earnings do not depend on whether he/she will enroll in the training program in the future, so the no anticipation assumption is reasonable. We assume that there is a common trend for the change in logarithmic earnings. The proportion of weeks employed is unlikely to be confounded with treatment or outcome, as the work time mainly depends on how long the participant spends on training. Therefore, we assume sequential ignorability. Positivity and consistency naturally hold.
Based on the proposed method, the total effect (average treatment effect in the treated group) is estimated at 0.0609 (s.e. 0.0623, 95%CI [-0.0612, 0.1830]). The natural indirect effect is estimated at -0.0447 (s.e. 0.0436, 95%CI [-0.1302, 0.0408]), indicating that joining the training program slightly reduces earnings by reducing work time, while this effect is insignificant. The natural direct effect is estimated at 0.1055 (s.e. 0.0445, 95%CI [0.0183, 0.1927]), indicating that joining the training program significantly leads to higher earnings if the participant had not reduced work time ($P=0.0178$).
Similarly, we take the actual training in the second year as the treatment, the weekly earnings in the second year as the pre-treatment outcome, the weekly earnings in the third year as the post-treatment outcome, and the proportion of weeks employed in the third year as the mediator. After adjusting for covariates, the average weeks employed is 1.5% higher in the treated group, and the increase in logarithmic earnings is 41.8% higher in the treated group. Based on the proposed method, the total effect is estimated at 0.4129 (s.e. 0.0452, 95%CI [0.3243, 0.5015]), the natural indirect effect is estimated at 0.0295 (s.e. 0.0143, 95%CI [0.0014, 0.0575]), and the natural direct effect is estimated at 0.3834 (s.e. 0.0429, 95%CI [0.2993, 0.4675]). All the effects are significant. The direct effect indicates that training has a positive effect on earnings. The indirect effect indicates that participants can work with ease in the third year after getting trained.
Since a large proportion of $M$ is observed at 0 and 1, we further conduct a sensitivity analysis by categorizing $M$ into four levels: 0 for unemployed ($M=0$), 1 for part-time employed ($0<M\leq0.5$), 2 for almost employed ($0.5<M<1$), and 3 for full-time employed ($M=1$). In the first year after training, the total effect is estimated at 0.0621 (s.e. 0.0623, 95%CI [-0.0600, 0.1841]), the natural indirect effect is estimated at -0.0704 (s.e. 0.0438, 95%CI [-0.1562, 0.0155]), and the natural direct effect is estimated at 0.1324 (s.e. 0.0435, 95%CI [0.0472, 0.2176]). In the second year after training, the total effect is estimated at 0.4135 (s.e. 0.0452, 95%CI [0.3249, 0.5020]), the natural indirect effect is estimated at 0.0270 (s.e. 0.0139, 95%CI [-0.0003, 0.0543]), and the natural direct effect is estimated at 0.3865 (s.e. 0.0430, 95%CI [0.3022, 0.4708]). The substantial conclusion remains the same. Table (ref) lists the estimates obtained using the proposed method.
To further investigate how the treatment effect varies with the intensity of employment, we estimate the controlled direct effect (CDE) curve $\tau_{DE}(m)$ as a function of $m$ ranging from 0 to 1. This represents the expected causal effect of training on potential earnings growth for the treated population had they fixed their employment level at $m$. Because earnings are zero when unemployed, the CDE is zero at $m=0$, reflecting the absence of an economically meaningful growth effect when no employment occurs. Figure (ref) presents the resulting CDE curves for both continuous and ordinal measures of employment intensity over two post-treatment years. In Year 1, the estimated CDE is generally small and fluctuates around zero across most values of $m$. For the continuous mediator, the confidence interval covers the null line. A similar pattern appears in the ordinal specification. Nevertheless, the natural direct effect is significant in an average sense. In Year 2, CDE exhibits clearer evidence of a positive direct effect of training. For both the continuous and ordinal mediators, the estimated CDE is significantly positive across a broad range of employment levels. This pattern indicates that, in the longer run, the Job Corps program enhances earnings capacity conditional on employment, consistent with delayed human capital accumulation beyond its impact on employment participation alone.
In this paper, we studied the causal mediation analysis in difference-in-differences. Under a set of assumptions, we showed that the total effect, natural indirect effect, and natural direct effect are identifiable in the treated group. Based on semiparametric theory, we derived the efficient influence function for the causal estimand and then proposed estimators that are multiply robust and nonparametrically efficient. Furthermore, we illustrated a practical strategy for estimation and inference using parametric models. The proposed framework addresses the challenges of mediation analysis in observational studies by considering difference-in-differences. The proposed framework can be extended to staggered designs. The parallel trends assumption should condition on the counterfactual mediator process. The sequential ignorability can be generalized by considering adjacent periods. However, nonparametrically efficient estimation is very challenging because every period has multiple working models. To reduce the computational burden, model-based methods would be an alternative, with a trade-off between statistical efficiency and model flexibility.
A limitation of causal mediation analysis is that sequential ignorability is difficult to interpret and impossible to verify, which means the mediator must be randomized. There have been approaches to relax sequential ignorability, such as imposing structural models for the potential outcomes, where the identifiability of treatment effects is reduced to parameter identifiability zheng2015causal. Alternative frameworks to natural effects have also been proposed, such as randomized interventional effects, in which the mediator is assumed to be manipulable according to a given distribution vanderweele2017mediation. When auxiliary data are available, instrumental variables help to identify mediated effects frolich2017direct, rudolph2024using.
Another limitation of the proposed estimation method is that the conditional density of the mediator is hard to estimate. Although we transformed the ratio of conditional densities into the ratio of pseudo-propensity scores, the model does not reflect the data-generating mechanism because it involves post-treatment variables. Therefore, the propensity score model and the pseudo-propensity score model may not be compatible. Complex nonparametric models, on the other hand, may not satisfy the Donsker condition, so cross-fitting is a possible approach at a higher computational cost and a larger finite-sample variation.
In the Job Corps Study, some participants' earnings are zero, a phenomenon known as the truncation-by-death problem zhang2009likelihood. Zero earnings may result from self-selection in the job market, with a higher expected wage than the offered wage. Our analysis did not address this problem.
The authors declare no conflict of interest.
The data that support the findings of this paper are publicly available in the R package “causalweight”.