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.
61,858 characters · 12 sections · 29 citation commands
Doubly Robust Estimation of Treatment Effects in Staggered Difference-in-Differences with Time-Varying Covariates
Difference-in-differences (DiD) is a fundamental econometric technique that estimates the causal effect by comparing changes in outcomes over time between treated and control groups. The classical two-group and two-period setup has been extensively studied in the past decades dimick2014methods, roth2023s, wang2024advances. By assuming parallel trends for the counterfactual outcomes under control, the mean potential outcomes are imputed for treated units. The average treatment effect on the treated (ATT) is then identified by restricting the target population to the treated group. Regression and weighting methods can be used to estimate the ATT. To increase model flexibility, researchers proposed semiparametric models abadie2005semiparametric, athey2006identification, wooldridge2010econometric. Doubly robust and locally efficient estimators are also proposed based on the semiparametric theory, which provide protection for estimation consistency against model misspecification sant2020doubly, deng2025improved.
In studies with multiple periods, treatment can be initiated in different periods, referred to as a staggered design athey2006identification, goodman2021difference, athey2022design. The two-way fixed effects (TWFE) model provides the most straightforward estimator and is therefore most commonly used in empirical studies de2020two, imai2021use. The coefficient associated with the treatment indicator is interpreted as the average treatment effect under the correct model specification. When the treatment effect is heterogeneous across groups and periods, the coefficient estimated from the TWFE model is a weighted average of group-period treatment effects, possibly involving negative weights goodman2021difference. As a result, the interpretation of ATT has sparked extensive discussion. To relax model restrictions, various estimators have been proposed in staggered difference-in-differences, including event studies sun2021estimating, de2023two, roth2023efficient, borusyak2024revisiting. The consistency of estimators relies on correct model specification. de2020two also observed that modeling the change in outcomes between adjacent periods may improve robustness. However, model-based methods conflate science and tools. In fact, the treatment effect is a quantity defined by counterfactuals, while models are merely tools for learning science holland1986statistics.
To avoid model-dependent estimands, a principled approach is to define the treatment effect at the group-period level. The group-period treatment effects can then be aggregated in various ways, yielding groupwise, periodwise, dynamic, and overall average treatment effects callaway2021difference, athey2022design. These aggregated effects evaluate the treatment effect from multiple dimensions. These aggregated effects are linear combinations of group-period treatment effects. Therefore, accurately estimating group-period treatment effects is a central objective in such analyses.
However, existing ATT estimators either lack efficiency or cannot incorporate time-varying covariates. In staggered adoption settings, all groups that have not yet been treated can serve as control groups to identify parallel trends. Since the proportions of these not-yet-treated groups vary over time, the information available for assessing parallel trends is time-dependent, leading to heterogeneity in estimation precision across periods. Existing doubly robust methods use a “generalized propensity score” to weight the time-invariant control group relative to the treated group, thereby incurring information loss about the time-varying not-yet-treated groups callaway2021difference. Although efficient estimators have been found in staggered designs, they cannot incorporate time-varying covariates chen2025efficient. Time-varying covariates partially explain the increase in outcomes in each period; therefore, the parallel trends of potential outcomes under control should be conditional on these time-varying covariates. However, time-varying covariates may depend on the observed outcome history, so their distributions differ across groups and time periods. Since ATT is defined on treated group-period cells, the distributional shift in time-varying covariates should be accounted for in identification and estimation. Current attempts to account for time-varying covariates mainly focus on two-way fixed effect models, which may suffer from bad control groups caetano2022difference, caetano2024difference. Double robustness can be achieved via augmented inverse probability weighting (AIPW), but, as with callaway2021difference, these estimators lose information about time-varying, not-yet-treated groups.
Motivated by the inverse variance weighting (IVW) technique commonly used in meta-analysis lee2016comparison, we propose an augmented inverse variance weighting (AIVW) estimator for treatment effects. Our estimator accommodates time-varying covariates and efficiently pools information across groups, offering increased efficiency over existing methods and benefiting from double robustness. In the homoskedastic case, the AIVW estimator reduces to an augmented inverse probability weighting (AIPW) estimator, which is more statistically efficient than the doubly robust estimator proposed by callaway2021difference. This AIPW estimator for the overall ATT can be readily computed as a weighted average of the residuals from the outcome regression model for treated units. The finite-sample efficiency of our estimator over existing methods is demonstrated through Monte Carlo simulations. As an empirical application, we investigate the impact of the parallel admission mechanism relative to immediate admission in the China National College Entrance Examination (NCEE). In our dataset, different provinces transitioned from immediate admission to parallel admission at different times, resulting in a staggered adoption design. For the outcome of interest, we use justified envy as a measure of fairness. The proposed estimator detects a significant improvement in fairness among students following the adoption of the parallel admission mechanism.
The remainder of this article is organized as follows. In Section (ref), we define the group-period ATTs within the potential outcomes framework and aggregate these cell ATTs into summary ATTs. Next, we present assumptions for identifying the treatment effects. In Section (ref), we propose doubly-robust estimators in the forms of augmented inverse variance weighting and augmented inverse probability weighting, accompanied by inference methods. We also illustrate how to implement the proposed estimators using standard working models. In Section (ref), we conduct simulation studies to compare the proposed estimators with existing methods and find that they can be more robust and efficient. In Section (ref), we study the effect of parallel admission reform on admission fairness in the National College Entrance Examinations of China. Finally, we discuss possible extensions in Section (ref).
In the setting of staggered difference-in-differences, suppose that there are $T+1$ periods. In the pre-treatment period ($t=0$), none of the units is treated. A unit may initiate treatment in any period $t \in \{1, \ldots, T\}$ after the pre-treatment period or never receive treatment throughout the experiment. Let $D_t$ be the treatment received in period $t$. We denote $G=g$ if the unit initiates treatment in period $g$, i.e., $D_0 = \cdots = D_{g-1} = 0$ and $D_g = 1$. Once the unit is exposed to treatment, we assume it will remain treated. If a unit never receives treatment, we denote either $G = T+1$ or $G = \infty$. Therefore, we can write $D_t = I(G \geq t)$. Both the group assignment $G$ and the vector $\bar{D} = (D_0, \ldots, D_T)$ can describe the entire treatment path. For each individual, let $Y_t(g)$ be the potential outcome in period $t$ if it were to initiate treatment in period $g$. We denote the potential outcome of never being treated as $Y_t(\infty)$. The individual treatment effect in period $t$ if the individual were assigned to group $g$ is $Y_t(g) - Y_t(\infty)$.
As a fundamental problem of causal inference, we cannot observe $Y_t(g)$ and $Y_t(\infty)$ simultaneously for any $g \leq T$ on a single individual. We aim to utilize available information across the entire population to identify the average treatment effect in a specific subpopulation. A common approach to identifying the treatment effect is to envision a counterfactual outcome under control for the treated units. Although the units can exhibit different responses to the exposure, it is natural to assume that there is a common trend in the outcomes if left untreated. In this way, the target estimand is the average treatment effect on the treated (ATT). We define the ATT in group $g$ and period $t$ as \[ \tau_{g,t} = E\{Y_t(g)-Y_t(\infty) \mid G=g\}. \]
Note that group-period ATTs are only defined in a triangular array. To aggregate the group-period ATT, we define the ATT in group $g$ as \[ \tau_g^{\text{grp}} = \frac{1}{T-g+1} \sum_{t=g}^{T} \tau_{g,t}, \] the ATT in period $t$ as \[ \tau_t^{\text{prd}} = \frac{1}{\sum_{g=1}^{t} P(G=g)} \sum_{g=1}^{t} P(G=g) \tau_{g,t}, \] the dynamic ATT of duration $s$ as \[ \tau_s^{\text{dyn}} = \frac{1}{\sum_{g=1}^{T-s+1} P(G=g)} \sum_{g=1}^{T-s+1} P(G=g) \tau_{g,g+s-1}, \] and the overall ATT as \[ \tau = \frac{1}{\sum_{g=1}^{T} (T-g+1)P(G=g)} \sum_{g=1}^{T}\sum_{t=g}^{T} P(G=g) \tau_{g,t}. \] Figure (ref) illustrates how these ATTs are aggregated and weighted by the proportion of individuals in each group-period cell of the target population. Suppose there are four periods ($t=0,\ldots,3$) and four groups ($G=1,\ldots,4$). The target population in the groupwise ATT $\tau_2^{\text{grp}}$ is $\{G=2\}$, the target population in the periodwise ATT $\tau_2^{\text{prd}}$ is $\{G\leq2\}$, and the target population in the dynamic ATT $\tau_1^{\text{dyn}}$ is $\{G\leq3\}$. From a super-population perspective, the treated group on which the overall ATT is defined consists of all unit-period pairs under treatment, $\{D_t = 1: t=1, \ldots, T\}$.
Let $X_t$ be the time-varying covariates in period $t$, which may include baseline covariates. The observed data in the sample $\{O_i = (G_i,X_{0i},Y_{0i},\ldots,X_{Ti},Y_{Ti}): i=1,\ldots,n\}$ consist of $n$ independent and identically distributed units with observation $O = (G, X_0, Y_0, \ldots, X_T, Y_T)$.
We define the increase in potential outcomes in period $t$ as $\Delta Y_t(g) = Y_t(g) - Y_{t-1}(g)$. We assume conditional no anticipation, saying that the conditional mean of the potential outcome under control does not depend on the time of future treatment imitation.
No anticipation is similar to the stable unit treatment value assumption (STUVA) in causal inference, which states that there is only one version of control. In a staggered design, no anticipation means that the increase in potential outcomes under control does not depend on when treatment will be received in the future.
Due to unmeasured confounding between the treatment assignment and potential outcomes, the distribution of $Y_t(g)$ may differ across groups $G$, even conditional on observed covariates. We assume parallel trends for the potential outcomes under control.
The parallel trends assumption is implied by sequential randomization, in which the potential outcome $Y_t(g)$ is independent of the treatment assignments $G$ conditional on baseline covariates. If $X_t \equiv X$ is time-invariant, then the parallel trends imply that $E\{Y_t(\infty) - Y_s(\infty) \mid X, G\}$ does not depend on $G$, for any $s,t\in\{0,\ldots,T\}$. Compared with the parallel trends assumption in the literature, our assumption allows time-varying covariates. The time-varying covariates $X_t$ in period $t$ can either include only measurements of covariates in period $t$ or aggregate historical information from the baseline to period $t$. They may partially explain the time trends in outcomes beyond an earlier period; therefore, parallel trends are applied to adjacent periods. For example, suppose that $Z$ is a baseline covariate and $Z_t$ is a time-varying covariate caetano2024difference. Consider a structural causal model \[ Y_t(\infty) = f_t(Z, Z_t) + U + \varepsilon_t, \] where $U$ is a time-invariant unmeasured confounder and $\varepsilon_t$ is a completely exogenous random error. Then the change in potential outcomes under control \[ \Delta Y_t(\infty) = f_t(Z,Z_t) - f_{t-1}(Z,Z_{t-1}) + \varepsilon_t - \varepsilon_{t-1} \] does not depend on the unmeasured confounder $U$. This implies that $E\{\Delta Y_t(\infty) \mid Z,Z_t,Z_{t-1},G\}$ does not depend on $G$. The set of time-invariant and time-varying covariates $X_t = (Z,Z_t,Z_{t-1})$ adjusts the parallel time trends.
Positivity states that there is a positive probability that each unit is assigned to each group. The support for time-varying covariates is identical across groups, so the information on time trends can be borrowed from untreated groups to treated groups, conditional on $X_t$. Consistency says that we can observe the potential outcome $Y_t(g)$ if the unit is in group $g$.
The no anticipation and parallel trends assumptions provide an opportunity to impute the mean potential outcome under control for the units receiving treatment. If $X_t \equiv X$ is time-invariant, it is well known that the group-period ATT is identifiable under the assumptions above, \[ \tau_{g,t} = E(Y_t - Y_0 \mid G=g) - E\{E(Y_t - Y_0 \mid X, G=\infty) \mid G=g\}. \] The first term is the mean of the observed increase in outcomes in group $g$, reflecting the joint effect of exposure and time trends. The second term captures the pure effect of time trends, identified from the control group after covariate adjustment. Note that the inner expectation gives the time trends conditional on covariates, and the outer expectation is the average time trends for units in group $g$. Since the group-period cells before exposure share the same parallel trends, $\tau_{g,t}$ is also identical to \[ \tau_{g,t} = E(Y_t - Y_{g-1} \mid G=g) - E\{E(Y_t - Y_{g-1} \mid X, G=\infty) \mid G=g\}. \] In the presence of time-varying covariates, to identify the time trends in group $g$ with a profile of time-varying covariates $(X_0,\ldots,X_t)$, the distribution of time-varying covariates should be shifted from control groups to group $g$ for each period. We present the identifiability in the following theorem.
The proof is given in Supplementary Material A.
Let $\mathbb{P}$ denote the measure of the true data generating process $\mathbb{P}(h) = \int h(o)f(o)do$, where $f(o)$ is the density function of the observed data $o$. Let $\mathbb{P}_n$ denote the empirical measure in the sample $\mathbb{P}_n(h) = n^{-1}\sum_{i=1}^{n}h(O_i)$. Let $\pi_g(X_t) = P(G=g \mid X_t)$ be the propensity score as a function of time-varying covariates. Let $\mu_{g,t}(X_t) = E(Y_t \mid X_t, G=g)$ be the mean outcome in group $g$ and period $t$ with covariates $X_t$. Let $\delta_{g,t}(X_t) = E(\Delta Y_t \mid X_t, G=g)$ be the mean increase in $\Delta Y_t = Y_t - Y_{t-1}$ outcomes between adjacent periods.
It is standard to use parametric or semiparametric models to fit the unknown models. For instance, we may estimate the propensity score using a proportional odds or ordinal probit model and estimate mean outcomes using linear regression. Suppose that the propensity score $\pi_g(X_t)$, mean outcome $\mu_{g,t}(X_t)$, mean change in outcomes $\delta_{g,t}(X_t)$ by $\widehat\pi_g(X_t)$, $\widehat\mu_{g,t}(X_t)$, and $\widehat\delta_{g,t}(X_t)$, respectively. Motivated by the identification formula, a regression-based (imputation-based) estimator is
Under the parallel trends assumption, the counterfactual change in mean outcomes under control $\delta_{\infty,k}(X_k)$ is identical across groups. Obviously, using the never-treated group $G=\infty$ to estimate $\delta_{\infty,k}(X_k)$ is inefficient because it does not fully exploit the information to assess parallel trends. Therefore, we can use all not-yet-treated data to estimate $\delta_{\infty,k}(X_k)$. Note that $E\{\delta_{G,k}(X_k) \mid X_k, G=g\} = E\{\delta_{\infty,k}(X_k) \mid X_k, G=g\}$ for $k < g$, so the regression-based estimator can also be written as
where $s$ is an arbitrary period $0 \leq s < g$.
In addition to regression-based estimators, we can also construct weighted estimators by appropriately accounting for the covariates' shift from control groups to treated groups,
where $s$ is an arbitrary period $0 \leq s < g$. Still, using the never-treated group as the reference is inefficient. In earlier periods, more groups had not yet been treated, so information from these groups can be used to identify parallel trends. The not-yet-treated groups $G>k$ in each period $k$ can be considered as a whole. By using not-yet-treated groups instead of the never-treated group, another weighted estimator is given by
Despite using all available untreated groups to identify parallel trends, the regression and weighted estimators mentioned above remain inefficient. This is because the parallel trends assumption only imposes a first-order moment condition on the mean outcomes. Suppose one of the not-yet-treated groups has a smaller variance in outcomes. In that case, it is more accurate to identify parallel trends using this group compared to other groups with larger variances. In meta-analysis, inverse variance weighting (IVW) was proposed to combine summary statistics from multiple sources, as it yields an efficient estimator when the source-specific estimators are obtained via maximum likelihood estimation, due to the additivity of information matrices. Another issue with the regression and weighted estimators is how to make inferences. Since the estimators involve fitted models, the uncertainty in fitted models may induce additional variation in the resulting estimators, which is hard to quantify.
Let $\sigma^2_{g,t}(X_t) = \operatorname{var}(\Delta Y_t \mid X_t, G=g)$ be the variance of the increased outcomes in group $g$ and period $t$ conditional on time-varying covariates. Inspired by IVW, we define weights \[ W_{l,k}(X_k) = \left[\sum_{s>k}\frac{\pi_s(X_k)}{\sigma^2_{s,k}(X_k)}\right]^{-1} \frac{\pi_l(X_k)}{\sigma^2_{l,k}(X_k)}. \] Let $\widehat{W}_{l,k}(X_k)$ be an estimate of $W_{l,k}(X_k)$ by plugging in the estimated $\pi_{\cdot}(X_k)$ and $\sigma^2_{\cdot,k}(X_k)$. We propose an augmented inverse variance weighting (AIVW) estimator of $\tau_{g,t}$ as
By aggregating $\widehat\tau_{g,t}$ according to groups, periods, or durations of treatment, we estimate the aggregated ATTs as
In practice, we can use a pooled model to estimate the outcome and variance models. A possible simplification of the estimator is to pretend a constant variance across groups and periods, $\sigma^2_{l,k}(X_k) = \sigma^2$ for all $l$ and $k$. Under this homoskedastic working model, the estimator is reduced to an augmented inverse probability weighting (AIPW) estimator
It is common that the models $\{\widehat\delta_{\cdot,\cdot}(\cdot), \widehat\pi_{\cdot}(\cdot), \widehat\sigma^2_{\cdot,\cdot}(\cdot)\}$ are not too complex. Without loss of generality, we assume $Y_t$ and $X_t$ have bounded total variation. Then it is natural to expect that $\{\widehat\delta_{\cdot,\cdot}(\cdot), \widehat\pi_{\cdot}(\cdot),\widehat\sigma^2_{\cdot,\cdot}(\cdot)\}$ are bounded. To establish consistency of the estimators, a condition is that \[ I(G=g) \sum_{k=g}^{t} \big\{ \Delta Y_k - \widehat\delta_{\infty,k}(X_k) \big\} - \sum_{k=g}^{t} I(G>k)\frac{\widehat\pi_g(X_k)}{\widehat\pi_G(X_k)} \widehat{W}_{G,k}(X_k) \big\{\Delta Y_k - \widehat\delta_{\infty,k}(X_k)\big\}, \] which appears in the expression of $\widehat\tau_{g,t}$, belongs to a Glivenko--Cantelli class, so that the law of large numbers can be applied van2013weak. Since this function is Lipschitz continuous with respect to $\{\widehat\delta_{\cdot,\cdot}(\cdot), \widehat\pi_{\cdot}(\cdot), \widehat\sigma^2_{\cdot,\cdot}(\cdot)\}$, it suffices to make sure these models are Glivenko--Cantelli. Trivially, parametric models are Glivenko--Cantelli. A wide range of models, such as generalized linear models, transformation models, kernel regression, and splines can make the Glivenko--Cantelli condition, and even the stronger Donsker condition, hold lu2007estimation, athey2019generalized, kuchibhotla2020efficient, martinez2023efficient. A notable property of the proposed AIVW estimator (and also the AIPW estimator) is its double robustness.
The consistency of estimators does not require the correct specification of $\sigma^2_{\cdot,\cdot}(\cdot)$. Let
The influence function of $\widehat\tau_{g,t}$, $\widehat\tau_g^{\text{grp}}$, $\widehat\tau_t^{\text{prd}}$, $\widehat\tau_s^{\text{dyn}}$, and $\widehat\tau$ are
respectively. The asymptotic distribution of the proposed estimator is given in the following theorem.
As a regularity condition, the Donsker condition requires that the working models are not too complex. Typical parametric, semiparametric, and nonparametric models (such as kernel regression and splines) can satisfy this condition. If machine learning approaches are employed, the Donsker condition can be replaced with sample splitting chernozhukov2018double, chernozhukov2022locally. The proofs of Theorems 2 and 3 are given in Supplementary Material A. Let $\widehat\varphi_{g,t}$, $\widehat\varphi_g^{\text{grp}}$, $\widehat\varphi_t^{\text{prd}}$, $\widehat\varphi_s^{\text{dyn}}$, and $\widehat\varphi$ be the fitted influence functions. To make an inference, we estimate the variances of these estimators by the sample variance of the fitted influence functions divided by the sample size, $\widehat\operatorname{var}(\widehat\tau_{g,t}) = n^{-1}\mathbb{P}_n\widehat\varphi_{g,t}^2$, $\widehat\operatorname{var}(\widehat\tau_g^{\text{grp}}) = n^{-1}\mathbb{P}_n\widehat\varphi_{g}^{\text{grp}2}$, $\widehat\operatorname{var}(\widehat\tau_t^{\text{prd}}) = n^{-1}\mathbb{P}_n\widehat\varphi_{t}^{\text{prd}2}$, $\widehat\operatorname{var}(\widehat\tau_s^{\text{dyn}}) = n^{-1}\mathbb{P}_n\widehat\varphi_{s}^{\text{dyn}2}$, and $\widehat\operatorname{var}(\widehat\tau) = n^{-1}\mathbb{P}_n\widehat\varphi^2$.
The AIVW estimator can be obtained by weighting the residuals from the outcome regression. For the group-period ATT, we define
Then the influence function of $\tau_{g,r}$ can be written as
To make the outcome regression more flexible, the treatment indicator term in the two-way fixed effect model can be split into several dummy variables representing the length of exposure, thereby capturing heterogeneous effects across exposure lengths sun2021estimating. Furthermore, we add the interaction terms between covariates and treatment as well as between covariates and periods, so the working model for outcome regression is
Note that interactions between groups and periods should not be included in the model to respect parallel trends in observed outcomes for untreated units. The coefficient $\alpha$ is no longer the treatment effect. Under this working model, the conditional mean outcome in group $g$ and period $t$ is \[ \mu_{g,t}(X_t) = \sum_{k=0}^{T-1} \alpha_k I(t-g=k) + \lambda_{t} + \gamma_{g} + \beta_1^{{\mathrm{\scriptscriptstyle T} }} X_t + \beta_2^{{\mathrm{\scriptscriptstyle T} }} X_tI(g\leq t) + \beta_3^{{\mathrm{\scriptscriptstyle T} }} X_tt, \] and the counterfactual mean outcome under control is \[ \mu_{g,t}^0(X_t) = \lambda_t + \gamma_g + \beta_1^{{\mathrm{\scriptscriptstyle T} }} X_t + \beta_3^{{\mathrm{\scriptscriptstyle T} }} X_tt. \] We approximate the conditional mean increase in outcomes $\delta_{\infty,k}(X_k)$ that appears in the influence functions by $\mu^0_{g,k}(X_k) - \mu^0_{g,k-1}(X_{k-1})$, which does not rely on $g$, in line with the parallel trends assumption. If there is no missingness, the sample for outcome regression includes $n(T+1)$ observations indexed by $\{(i,t): i=1,\ldots,n, t=0,\ldots,T\}$, which we pretend to be independent when fitting the model. Let $\widehat\mu_{g,t}^0(X_t)$ be the fitted counterfactual mean outcome under control in group $g$ and period $t$.
Let $\widehat\varepsilon_t$ be the residual of the outcome regression model. To calculate the weights, we further fit a regression for residuals, \[ \log\{(\widehat\varepsilon_t-\widehat\varepsilon_{t-1})^2\} = \lambda_t^* + \gamma_G^* + \beta^{*{\mathrm{\scriptscriptstyle T} }} X_t + \epsilon_t \] based on the sample $\{(i,t): D_{ti}=0, t>0\}$. Then we estimate the variance of increased outcomes by $\widehat\sigma_{g,t}^2(X_t) = \exp(\widehat\lambda_t^* + \widehat\gamma_g^* + \widehat\beta^{*{\mathrm{\scriptscriptstyle T} }} X_t)$. The $\widehat\sigma_{g,t}^2(X_t)$ is involved in estimating the weights $W_{G,t}(X_t)$ in the AIVW estimator. However, it is computationally intensive because we need to calculate $\widehat\sigma_{g,t}^2(X_t)$ for every $g$ and $t$ subject to $0<t<g\leq T$. To improve computational efficiency, we can assume independent and homoskedastic error terms, $\varepsilon_{t} \sim N(0, \sigma_{\varepsilon}^2)$, with an unknown $\sigma_{\varepsilon}^2$. Therefore, the working variance of increased outcomes $\sigma^2_{g,t}(X_t) = 2\sigma_{\varepsilon}^2$ is a constant, and the AIVW estimator is reduced to the AIPW estimator. Additionally, we fit the time-varying propensity score $\pi_g(X_t)$ by ordinal logistic regression, \[ \frac{P(G \leq k \mid X_t)}{P(G > k \mid X_t)} = \frac{\sum_{s\leq k}\pi_s(X_t)}{\sum_{s>k}\pi_s(X_t)} = \exp(\zeta_{kt0} + \zeta_{kt}^{{\mathrm{\scriptscriptstyle T} }}X_t), \] denoted by $\widehat\pi_g(X_t)$.
The fitted models $\widehat\pi_g(X_t)$ and $\widehat\sigma^2_{g,t}(X_t)$ are parametric, so the fitted influence function $\varphi_{g,r}$ belongs to a Donsker class. Under correct model specification, these fitted models converge at a rate of $O_p(n^{1/2})$. By plugging in fitted models $\widehat\pi_g(X_t)$ and $\widehat\sigma^2_{g,t}(X_t)$ into $H_{G,t}^{g,r}$, the estimator $\widehat\tau_{g,r}$ is obtained by solving the empirical estimating equation $\mathbb{P}_n\widehat\varphi_{g,r} = 0$, which is the empirical average of $\widehat{H}_{G,t}^{g,r}\{Y_t-\widehat\mu_{G,t}^0(X_t)\}$ in the sample $\{(i,t): G_i=g, t=r\}$. The asymptotic variance of $\widehat\tau_{g,r}$ is estimated by the sample variance of $\widehat\varphi_{g,r}$ in all units.
For the overall ATT, we define \[ H_{G,t} = \sum_{r=1}^{T}\sum_{g=1}^{r} H_{G,t}^{g,r}. \] In particular, if AIPW is used,
The function $H_{G,t}$ involves the propensity score through the time-varying odds ratio. Then the influence function of $\tau$ can be written as
Here $P(D_t=1) = \sum_{g=1}^{T}(T-g+1)P(G=g)/(T+1)$ is the probability of being treated for all $(i,t)$ pairs. By plugging in fitted models, the estimator $\widehat\tau$ is obtained by solving the empirical estimating equation $\mathbb{P}_n\widehat\varphi = 0$, which is the empirical average of $\widehat{H}_{G,t}\{Y_t-\widehat\mu_{G,t}^0(X_t)\}$ in the sample $\{(i,t): D_{ti}=1\}$. The asymptotic variance of $\widehat\tau$ is estimated by the sample variance of $\widehat\varphi$ in all units.
In this section, we conduct simulation studies to demonstrate the finite-sample performance of the proposed AIPW and AIVW estimators. We consider three competing methods: (1) the TWFE model, (2) the doubly robust estimator that uses never-treated units as the control group (DRnt), and (3) the doubly robust estimator that uses not-yet-treated units as the control group. The latter two estimators were proposed in callaway2021difference, where the overall ATT is the simple aggregation of group-period ATTs.
Suppose that there are five periods, $t=0,\ldots,4$. We generate two independent time-invariant covariates, $Z_{1}$ and $Z_{2}$, following the standard normal distribution. In addition, we generate an exogenous time-varying covariate $Z_{3,t}$ following a normal distribution $N(0,(1+0.1t)^2)$, which is independent across periods and units. We denote $X_t=(Z_{1},Z_{2},Z_{3,t},Z_{3,t-1})$ with complementarily defining $Z_{3,-1}=0$.The period to initiate treatment is determined at baseline, with the propensity score generated from an ordinal logistic model,
where $(\zeta_1,\zeta_2,\zeta_3,\zeta_4) = (-1.5,-0.5,0,1)$. The probability of never-treated is $\pi_{\infty}(X_0) = \mathrm{P}(G>4 \mid X_0)$. The distribution of $(Z_{1},Z_{2})$ is different across groups, while the distribution of $(Z_{3,t},Z_{3,t-1})$ is identical across groups.
We consider two scenarios for generating the potential outcomes. The first scenario is the homogeneous case, where the treatment effect does not depend on group or time,
where $\xi \sim N(0,1)$ is an individual random effect and $u_t \sim N(0,1)$ is an error term independent across periods and units. The potential outcome $Y_t(g)$ is not independent of the treatment assignment $G$, but the parallel trends assumption holds, as the conditional change in potential outcomes under control $\{Y_t(\infty) - Y_{t-1}(\infty) \mid X_t, G\} = 0.5Z_{1}+0.2$ does not depend on $G$. The true ATT is $1$, the coefficient of $I(g\leq t)$. The second scenario is the heterogeneous case, where the treatment effect depends on group and time,
where $\xi \sim N(0,1)$ is an individual random effect and $u_t \sim N(0,1)$ is an error term independent across periods and units. The parallel trends assumption holds, as the conditional change in potential outcomes under control $E\{Y_t(\infty) - Y_{t-1}(\infty) \mid X_t, G\} = 0.5Z_{1}+0.2$ does not depend on $G$. The treatment effect is larger for a longer exposure. Noting that $E(Z_{3,t}\mid G)=0$, the true ATT is $\sum_{g=1}^{4}\pi_g\sum_{t=g}^{4}(1+0.5(t-g))/\sum_{g=1}^{4}(5-g)\pi_g$, in which the proportion of each group $\pi_g = \int_{-\infty}^{+\infty}\pi_g(x)\exp(-2x^2)/\sqrt{\pi/2}dx$ can be numerically calculated.
Let the sample size $n \in \{100, 500, 2000\}$. We replicate the data-generating process 1000 times and estimate the overall ATT by five methods in each replicate. Panel (A) of Table (ref) shows the average bias, standard deviation (SD), average standard error (SE), and the empirical coverage percentage (CP) of the nominal 95% confidence interval. The standard deviation of the estimator from the TWFE model is small because the model is oversimplified. However, TWFE is biased due to model misspecification. In line with the theory, AIPW and AIVW are both asymptotically unbiased. The coverage percentages of the confidence intervals associated with AIPW and AIVW are close to the nominal level. Since the AIVW estimator involves more models, the fitted models in AIVW exhibit greater finite-sample variation than those in AIPW. The standard deviations of AIPW and AIVW are smaller than those of DRnt and DRny. Although DRnt and DRny can incorporate time-varying covariates in outcome regression, these methods do not provide a theoretical guarantee for consistent estimation. When treatment effects are heterogeneous, DRnt and DRny introduce substantial bias.
To further compare the estimators, we consider four alternative data-generating mechanisms. In Scenarios 3 and 4, we generate the error term $u_t \sim N(0, \exp(0.2Z_{2}+0.3t-0.3G))$, while the mean outcomes remain the same as in Scenarios 1 and 2. The bias, standard deviation, average standard error, and empirical coverage percentage are shown in Panel (B) of Table (ref). Consistent estimation of ATT does not require modeling the outcome variance for the proposed methods, so AIPW and AIVW are both asymptotically unbiased. The outcomes in the never-treated group have a smaller variance than those in the not-yet-treated groups. DRny uses the never-treated group as the reference, which entails observations with greater variability and therefore reduces efficiency compared with DRnt. Since AIPW estimates time trends period-by-period, the larger variation in earlier periods will not lead to substantial variability in the estimation of time trends in later periods. Beyond AIPW, AIVW weights observations according to their variability, allowing observations with lower variability to contribute more to the estimation of time trends.
In Scenarios 5 and 6, we generate the error term $u_t = \sum_{k=0}^{t} u^*_k$, where the $u^*_k$'s are independent with $u^*_k \sim N(0, \exp(0.2Z_{2}+0.3t-0.3G))$. The mean outcomes remain the same as in the previous scenarios. The bias, standard deviation, average standard error, and empirical coverage percentage are shown in Panel (C) of Table (ref). The proposed AIPW and AIVW are asymptotically unbiased and exhibit smaller variance than DRnt and DRny. The data in these two scenarios are generated according to the mechanism described in Remark (ref), so AIVW is the most efficient. The confidence intervals of AIPW and AIVW have coverage percentages close to the nominal level. However, AIVW requires more computational time than AIPW because an additional variance model must be fitted.
To illustrate the inference for cell ATTs, we set the sample size $n=500$ in Scenario 4. We plot the average estimate, standard deviation (SD), and average standard error (SE) for the group-period ATTs in the first row of Figure (ref). Note that the true $\tau_{g,t} = 1+0.5(t-g)$ for $g \ge t$. In addition, we plot the groupwise ATTs, periodwise ATTs, and dynamic ATTs over time. Groupwise ATT in group $g$ averages the group-period ATTs for $t=g,\ldots,T$. Periodwise ATT in period $t$ averages the group-period ATTs for $g=1,\ldots,t$ according to group proportions $\pi_g$. Dynamic ATT after $r$ periods of exposure averages the group-period ATTs for $G-t=r$. The true values are plotted in black points. The vertical lines represent the mean estimate $\pm$ 1.96 SD or SE. The average estimate, standard deviation, and average standard error for ATTs based on AIPW and AIVW are very close. We find that the bias is close to zero, and the average standard error is close to the standard deviation, indicating that the inference is accurate.
Education plays a pivotal role in shaping individual labor market outcomes and promoting social mobility. Around the world, many countries employ centralized matching mechanisms to assign students to schools and universities. One of the most well-known of these is the Immediate Acceptance (IA) mechanism, also known as the Boston Mechanism, in which students apply to schools in order of their preferences abdulkadirouglu2003school. In each round, schools immediately and irrevocably accept their highest-priority applicants (based on factors like exam scores) until all available spots are filled. In recent decades, China has transitioned from using IA to a parallel mechanism, a variant of the Deferred Acceptance (DA) mechanism chen2020empirical, ha2020college. Students list several universities within each choice band. Within a choice band, the parallel mechanism allows temporary assignments. Previous research suggests that the parallel mechanism improves stability and reduces manipulability kang2020matching. However, existing studies rely on two-way fixed effects models or event studies, which limit their statistical efficiency and robustness.
The staggered adoption of the parallel mechanism across different provinces provides a unique natural experiment for evaluating its effects relative to the implementation of IA. Leveraging administrative data on the National College Entrance Examination (NCEE, also known as Gaokao) in China from 2007 to 2011, we aim to analyze the impact of these staggered provincial reforms on admission fairness. We retain 27 provinces in our sample after excluding those that neither used IA nor parallel admission. In 2007, all provinces used IA. The number of provinces reforming to a parallel admission mechanism was 3, 10, 6, and 2 in 2008, 2009, 2010, and 2011, respectively. Additionally, six provinces had not been reformed by 2011. We split the stem tracks (stem or non-stem) in each province, resulting in a sample size of $n=54$. We call each province-track in each year a game.
We use the justified envy (JE) as the outcome to quantify fairness. JE is a standard measure of fairness in the school choice literature balinski1999tale, abdulkadirouglu2003school, kamada2024fair. We say that student $i$ justifiably envies student $j$ for school $s$ if $i$ would rather be assigned to school $s$, where some student $j$, who has a lower priority (i.e., lower score) than $i$, is assigned. In this case, student $i$ is a blocking student; student $i$ and school $s$ are a blocking pair; student $i$, $j$, and school $s$ are a blocking triplet. We consider four numerical outcomes for justified ency: (1) BS, the total number of blocking students in a game; (2) BP, the total number of blocking pairs in a game; (3) TE, the total tridimensional envy of blocking triplets in a game, where the tridimensional envy is defined as the diagonal distance of student quantiles and university tiers in a blocking triplet; (4) BT, the total number of blocking triplets in a game. Formal definitions of these measures can be found in kang2020matching.
In the analysis, we controlled province-level log GDP per capita and population as baseline covariates in the propensity score model. For the outcome regression model, we include province-level log GDP per capita, population, and game-level track (stem or non-stem). Time-varying covariates include a dummy variable indicating whether to submit the preference after taking the exam, a dummy variable indicating whether to submit the preference after knowing the scores, and the number of students in a game. We check the parallel trends assumption in Supplementary Material B by plotting the residuals. We do not find evidence against parallel trends.
Table (ref) presents the estimated overall ATTs for these four measures of justified envy. Although the estimated effects are significant under TWFE, they are unreliable due to potential model misspecification. The doubly robust methods (DRnt and DRny) are potentially biased, as indicated by the simulation. The proposed AIVW yields significant estimated ATTs for BS, BP, TE, and BT at the 0.05 significance level, indicating that the parallel mechanism effectively reduces justified envy relative to IA. The absolute point estimates from AIPW are slightly smaller than those from AIVW. Based on AIPW, the treatment effects on BP, TE, and BT are significant at the 0.05 significance level. Despite the slight difference between AIPW and AIVW, the substantial conclusion remains the same.
The estimated overall ATTs by the proposed methods indicate that the extent of justified envy decreases after the implementation of the parallel mechanism. Figure (ref) shows the periodwise ATTs since 2007 based on AIPW and AIVW. The point estimates and confidence intervals by AIPW and AIVW are similar. In the first year of introducing the parallel mechanism in NCEE, the total number of blocking students and blocking triplets increased, but the magnitudes are insignificant. After the parallel mechanism is implemented, parents and students can learn from the experience of other provinces. They become more familiar with the process and perform better on college applications year by year. The treatment effects on justified envy become larger after two years of implementing the parallel mechanism, reflecting learning from other provinces' experience. In Supplementary Material B, we present groupwise and dynamic ATTs. We do not find obvious patterns for groupwise and dynamic ATTs. In Supplementary Material B, we also consider standardized justified envy, defined as the observed justified envy divided by the maximum possible justified envy in each game. All methods indicate that the treatment has significant effects on all groupwise, periodwise, dynamic, and overall ATTs.
Time-varying covariates are common in observational studies, while existing studies mainly focus on identifying and estimating treatment effects with baseline covariates. In this article, we propose an augmented inverse variance estimator (AIVW) for the ATT in staggered difference-in-differences with time-varying covariates. Compared to existing estimators, the proposed estimator can be more efficient because it leverages information from multiple not-yet-treated groups to determine parallel trends on a period-by-period basis. Under homoskedasticity, the AIVW estimator reduces to the augmented inverse probability weighting (AIPW) estimator. Based on group-period ATTs, we aggregate the cell treatment effects into periodwise, groupwise, and dynamic effects. For AIPW and AIVW, the influence functions of these aggregated ATTs are weighted sums of the influence functions of the cell influence functions. Statistical inference is straightforward using influence functions.
The AIPW estimator is computationally efficient because it requires only fitting an outcome regression model and a propensity score model. The overall ATT can be easily calculated as a weighted average of the fitted model residuals across treated unit-period pairs, as shown in Equation (ref). The finite sample of treated unit-period pairs corresponds to the treated portion of group-period cells in the super-population. Simulation studies indicate that the AIPW estimator performs comparably to the AIVW estimator. In small sample sizes, AIPW suffers less from finite-sample variation compared to AIVW because AIPW avoids additional uncertainty in fitted models. Other studies have observed that simpler models can outperform estimators based on efficient influence functions in terms of standard deviation when data are limited wang2023model. Given the computational efficiency and finite-sample performance of the AIPW estimator, we recommend using the linear model with equal variances for the error terms as the working model.
Depending on the data-generating mechanism, the AIVW or AIPW estimator is not necessarily the most efficient. One may be interested in finding the best estimator for estimating the ATT. However, this is not a simple task because the parallel trends assumption imposes restrictions on observed data. The model for observed data is semiparametric rather than nonparametric. Since the parallel trends assumption does not directly restrict the data-generating mechanism but instead imposes a moment condition, the efficient influence function for ATT is challenging to derive. A possible approach is to limit the semiparametric model space by only considering a subset of estimators, such as linear predictors. Another possible approach is to refine the parallel trends assumption on the data-generating mechanism, for example, by assuming conditional parallel trends given all history, so that the parallel trends directly restrict the likelihood.
The data that supports the findings of this work is available from the corresponding author upon request.