{did2s}: Two-Stage Difference-in-Differences

Recent work has highlighted the difficulties of estimating difference-in-differences models when treatment timing occurs at different times for different units. This article introduces the R package did2s which implements the estimator introduced in Gardner (2021). The article provides an approachable review of the underlying econometric theory and introduces the syntax for the function did2s. Further, the package introduces a function, event_study, that provides a common syntax for all the modern event-study estimators and plot_event_study to plot the results of each estimator.


Introduction
A rapidly growing econometric literature has identified difficulties in traditional difference-in-differences estimation when treatment turns on at different times for different groups and when the effects of treatment vary across groups and over time (Callaway and Sant'Anna, 2020;Sun and Abraham, 2020;Goodman-Bacon, 2018;Borusyak et al., 2021;de Chaisemartin and D'Haultfoeuille, 2019).Gardner (2021) proposes an estimator of the two-way fixed-effects model that is quick and intuitive.The estimator relies on the standard two-way fixed-effect model (see the following section) and forms an intuitive estimate: the average difference in outcomes between treated and untreated units after removing fixed unit-and time-invariant shocks.
This article first discusses the modern difference-in-differences theory in an approachable way and second discusses the software package, did2s, which implements the two-stage estimation approach proposed by Gardner (2021) to estimate robustly the two-way fixed-effects (TWFE) model.There are two notable technical features of this package.First, did2s utilizes the incredibly fast package, fixest (Bergé, 2018), which can estimate regressions with a high number of fixed-effects very quickly.Since there are a few alternative TWFE event-study estimators implemented in R, each with their own syntax and data formatting requirements, the package also has a set of functions that allow quick estimation and plotting of every alternative event study estimator using a standardized syntax.This allows for easy comparison between the results of different methods.

Difference-in-Differences Theory
Researchers commonly use the difference-in-differences (DiD) methodology to estimate the effects of treatment in the case where treatment is non-randomly assigned.Instead of random assignment giving rise to identification, the DiD method relies on the so-called "parallel trends" assumption, which asserts that outcomes would evolve in parallel between the treated and untreated groups in a world where the treated were untreated.This is formalized with the two-way fixed-effects (TWFE) model.In a static setting where treatment effects are constant across treatment groups and over time, researchers estimate the static TWFE model: where y igt is the outcome variable of interest, i denotes the individual, t denotes time, and g denotes group membership where "group" is defined as all units that start treatment at time g. 1 µ g is a vector of time-invariant group characteristics, η t is a vector of shocks in a given time period that is experienced by all individuals equally, and D gt is an indicator for whether initial-treatment group g is receiving treatment in period t, i.e.D gt ≡ 1(g ≤ t).The coefficient of interest is τ, which is the (constant) average effect of the treatment on the treated (ATT).If it is indeed true that the treatment effect is constant across groups and over time, then the estimate formed by estimating the static TWFE model will be consistent for τ.
However, treatment effects are not constant in most settings.The magnitude of a unit's treatment effect can differ based on group status g (e.g. if groups that benefit more from a policy implement it earlier) and treatment duration (e.g. if treatment effects grow as the policy has been in place for longer periods).Therefore to enrich our model, we allow heterogeneity in treatment effects across g and t by introducing the group-time average treatment effect, τ gt .Correspondingly, we modify the TWFE model as follows: The key difference is that treatment effects are allowed to differ based on group status g and time period t.Estimating any individual τ gt may not be desirable since there would be too few observations.Instead, researchers aggregate group-time average treatment effects into the overall average treatment effect, τ, which averages across τ gt : where N gt denotes the number of (post-treatment) group-time pairs, {g, t}, observed in our sample.The natural question is, "does the static TWFE model produce a consistent estimate for the overall average treatment effect?"Except for a few specific scenarios, the answer is no (Sun and Abraham, 2020;Goodman-Bacon, 2018;Borusyak et al., 2021;de Chaisemartin and D'Haultfoeuille, 2019).
One way of thinking about this disappointing result is through the Frisch-Waugh-Lovell (FWL) theorem (Frisch and Waugh, 1933).This theorem says that estimating the Static TWFE model is equivalent to estimating where Dgt denotes the residuals from regressing D gt on µ g and η t , μg and ηt are estimates for the group and time fixed-effects, respectively.The left-hand side of this equation, under a parallel trends restriction on the error term ε it , is our estimate for τ gt .Therefore, the FWL theorem tells us estimating the static TWFE model is equivalent to estimating2 τgt = τ Dgt + εgt The resulting estimate for τ can be written as: where w gt is the weight put on the corresponding τgt .Results of Gardner (2021), Borusyak et al. (2021), and de Chaisemartin and D'Haultfoeuille (2019) all characterize the weights w gt from this regression.
There are only two cases where the τ is a consistent estimate for the overall average treatment effect.First, when treatment occurs at the same time for all treated units, then w gt is equal to 1/N gt for all {g, t} and therefore τ is a consistent estimate for the overall average treatment effect.The other scenario when τ estimates the overall average treatment effect is when τ gt is constant across group and time, i.e. τ gt = τ.Since the weights, w gt , always sum to one, we have that τ = ∑ w gt τgt → ∑ w gt τ = τ.The above cases are not the norm in research.If there is heterogeneity in group-time treatment effects and differences in when units get treated, then τ is not a consistent estimate for the average treatment effect τ. Instead, τ will be a weighted average of group-time treatment effects with some weights, w gt , being potentially negative.This yields a treatment effect estimate that does not provide a good summary of the "average" treatment effect.It is even possible for the sign of τ to differ from that of the overall average treatment effect.This would occur, for example, if negative weights are placed primarily on the largest (in magnitude) group-time treatment effects.
To summarize the modern literature, the fundamental problem faced in estimating the TWFE model is the potential negative weighting.The proposed methodology in Gardner (2021) is based on the fact that if τgt is regressed on D gt , instead of Dgt , the resulting weights would be exactly equal to 1/N gt and the coefficient of D gt would estimate the overall average treatment effect.

Event-study Estimates
Researchers have attempted to model treatment effect heterogeneity by allowing treatment effects to change over time.To do this, they introduce a (dynamic) event-study TWFE model: where D k gt are lags/leads of treatment (k periods from initial treatment date).The coefficients of interests are the τ k , which represent the average effect of being treated for k periods.For negative values of k, τ k are known as "pretrends," and represent the average deviation in outcomes for treated units k periods away from treatment, relative to their value in the reference period.These pre-trend estimates are commonly used as a test of the parallel counterfactual trends assumption.
Our goal is to estimate the average treatment effect of being exposed for k periods, an average of τ gt for only the set of {g, t} where k periods have elapsed since g, i.e. t − g = k: where the sum is over {g, t} with t − g = k and N k gt is the count of {g, t} pairs that satisfy that condition.The results of Sun and Abraham (2020) show that even though we allow for our average treatment effects to vary over time τ k , the negative weighting problems would arise if units are treated at different times and there is group-heterogeneity in treatment effects.Similar to the static TWFE model, the estimates of τ k from running the event-study model form non-intuitively weighted averages of τ gt with w k gt = N k gt .Even worse, the group-time treatment effects for t − g = k will be included in the estimate of τk .Hence, the need for a robust difference-in-differences estimator remains even in the event-study model.Gardner (2021) proposes an estimator to resolve the problem with the two-way fixed-effects approaches.Rather than attempting to estimate the group and time effects simultaneously with the ATT (causing D it to be residualized), Gardner's approach proceeds from the observation that, under parallel trends, the group and time effects are identified from the subsample of untreated/not-yet-treated observations (D gt = 0).This suggests a simple two-stage difference-in-differences estimator:

Two-stage Difference-in-Differences Estimator
1. Estimate the model using the subsample of untreated/not-yet-treated observations (i.e., all observations for which D gt = 0), retaining the estimated group and time effects to form the adjusted outcomes ỹigt ≡ y igt − μg − ηt .
2. Regress adjusted outcomes ỹigt on treatment status D gt or D k gt in the full sample to estimate treatment effects τ or τ k .
To see why this procedure works, note that parallel trends implies that outcomes can be expressed as is the average treatment effect for group g in period t3 and τ = E(τ gt |D gt = 1) is the overall average treatment effect4 .Note from parallel trends, E(ε igt |D gt , g, t) = 0. Rearranging, this gives Suppose you knew the time and group fixed-effects and were able to directly observe the left-hand side (later we will estimate the left-hand side).Regressing the adjusted y variable, on D gt will produce a consistent estimator for τ.To see this, note that E[(τ gt − τ)D gt | D gt ] = 0. Hence, the treatment dummy is uncorrelated with the omitted variable and the average treatment effect is identified in the second-stage.Since we are not able to directly observe µ g and η t , we estimate them using the untreated/not-yet-treated observations in the first-stage.However, standard errors need adjustment to account for the added uncertainty from the first-stage estimation.
This approach can be extended to dynamic models by replacing the second stage of the procedure with a regression of residualized outcomes onto the leads and lags of treatment status, D k gt , k ∈ {−L, . . . ,K}.Under parallel trends, the second-stage coefficients on the lags identify the overall average effect of being treated for k periods (where the average is taken over all units treated for at least that many periods).The second-stage coefficients on the leads identify the average deviation from predicted counterfactual trends among units that are k periods away from treatment, which under parallel trends should be zero for any pre-treatment value of k.Hence, the coefficients on the leads represent a test of the validity of the parallel trends assumption.

Inference
The standard variance-covariance matrix from the second-stage regression will be incorrect since it fails to account for the fact that the dependent variable is generated from the first-stage regression.However, this estimator takes the form of a joint generalized method of moments (GMM) estimator whose asymptotic variance is well understood (Newey and McFadden, 1986).Specifically, the estimator takes the form of a two-stage GMM estimator with the following two moment conditions: where X 1 is the matrix of group and time fixed-effects, X 10 corresponds to the matrix X 1 , but with rows corresponding to observations for which D gt = 1 replaced with zeros (as only observations with D gt = 0 are used in the first stage) and X 2 is the matrix of treatment variable(s).The first equation corresponds with the first stage and the second equation corresponds with the second stage.From Theorem 6.1 of Newey and McFadden (1986), the asymptotic variance of the two-stage estimator is where from our moment conditions, we have: This can be estimated using where and matrices indexed by g correspond to the gth cluster.

The did2s Package
The did2s package introduces two sets of functions.The first is the did2s command which implements the two-stage difference-in-differences estimator as described above.The second is the event_study and plot_event_study commands that allow individuals to implement alternative 'robust' estimators using a singular common syntax.

The did2s Command
The command did2s implements the two-stage difference-in-differences estimator following Gardner (2021).The general syntax is did2s(data, yname, first_stage, second_stage, treatment, cluster_var, weights = NULL, bootstrap = FALSE, n_bootstraps = 250, verbose = TRUE) and full details on the arguments is available in the help page, available by running ?did2s.There are a few arguments that are worth discussing in more detail.
The first_stage and second_stage arguments require formula arguments.These formulas are passed to the fixest::feols function from fixest and can therefore utilize two non-standard formula options that are worth mentioning (Bergé, 2018).First, fixed-effects can be inserted after the covariates, e.g.~x1 | fe_1 + fe_2, which will make estimation much faster than using factor(fe_1).with fixest::etable and event-study plots with fixest::iplot/fixest::coefplot.The fixest::i function is better integrated with these functions as we will see below.
The option treatment is the variable name of a 0/1 variable that denotes when treatment is active for a given unit, D gt in the above notation.Observations with D gt = 0 will be used to estimate the first stage, which removes the problem of treatment effects contaminating estimation of the unit and time fixed-effects.However, as an important note, if you suspect anticipation effects before treatment begins, the treatment variable should be shifted forward by x periods for observations to prevent the aforementioned contamination.For example, if you suspect that units could experience treatment effects 1 period ahead of treatment (a so-called anticipatory effect), then the treatment should begin one period ahead.These anticipation effects can be estimated, after adjusting the treatment variable, by using a reference year of say, t = −2 and looking at the estimate for relative year −1.
Example usage of did2s For basic usage, I will use the simulated dataset, df_het, that comes with the did2s package with the command data(df_het, package = "did2s") The data-generating process is displayed in Figure 1.The lines represent the mean outcome for each treatment group and the never-treated group.In the absence of treatment, each group is simulated to be on parallel trends.There is heterogeneity in treatment effects both within a treatment group over time and across treatment groups.
First, we will calculate a static difference-in-differences estimate using the did2s function.However, since there are dynamic treatment effects in this example, it is much better to estimate the dynamic effects themselves using an event-study specification.We will then plot the results using fixest::iplot, which plots coefficients corresponding to an i() variable.Note that rel_year is coded as Inf for never-treated units, so this has to be noted in the reference part of the formula.The event study estimates are found in Figure 2 and match closely to the true average treatment effects.For comparison to traditional OLS estimation of the event-study specification, Figure 3 plots point estimates from both methods.As pointed out by Sun and Abraham (2020), treatment effect heterogeneity between groups biases the estimated pre-trends.In the figure below, the OLS estimates appear to show violations of pre-trends even though the data was simulated under parallel pre-trends.# Legend legend(x=-20, y=3, col = c("steelblue", "#82b446", "grey60"), pch = c(20, 18, 20), legend = c("Two-stage estimate", "TWFE", "True Effect")) The event_study and plot_event_study command The command event_study presents a common syntax that estimates the event-study TWFE model for treatment-effect heterogeneity robust estimators recommended by the literature and returns all the estimates in a data.framefor easy plotting by the command plot_event_study.The general syntax is event_study( data, yname, idname, tname, gname, estimator, xformla = NULL, horizon = NULL, weights = NULL ) The option data specifies the data set that contains the variables for the analysis.The four other required options are all names of variables: yname corresponds with the outcome variable of interest; idname is the variable corresponding to the (unique) unit identifier, i; tname is the variable corresponding to the time period, t; and gname is a variable indicating the period when treatment first starts (group status).

Event study: Staggered treatment (comparison)
Relative time to treatment Estimate and 95% Conf.Int.There are five main estimators available and the choice is specified for the estimator argument and are described in Table 2. 5 The following paragraphs will aim to highlight the differences and commonalities between estimators.These estimators fall into two broad categories of estimators.First, did2s and didimputation (Butts, 2021) are imputation-based estimators as described above.Both rely on "residualizing" the outcome variable Ỹ = Y it − μg − ηt and then averaging those Ỹ to estimate the event-study average treatment effect τ k .These two estimators return identical point estimates, but differ in their asymptotic regime and hence their standard errors.
The second type of estimator, which we label 2x2 aggregation, takes a different approach for estimating event-study average treatment effects.The packages did (Callaway and Sant'Anna, 2021), fixest and staggered (Roth and Sant'Anna, 2021) first estimate τ gt for all group-time pairs.To estimate a particular τ gt , they use a two-period (periods t and g − 1) and two-group (group g and a "control group") difference-in-differences estimator, known as a 2x2 difference-in-differences. The particular "control group" they use will differ based on estimator and is discussed in the next paragraph.Then, the estimator manually aggregate τ gt across all groups that were treated for (at least) k periods to estimate the event-study average treatment effect τ k .These estimators do not all rely on the same underlying assumptions, so the rest of the table summarizes the primary differences between estimators.The comparison group column describes which units are utilized as comparison groups in the estimator and hence will determine which units need to satisfy a parallel trends assumption.For example, in some circumstances, treated units will look very different from never-treated units.In this case, parallel trends may only hold between ever-treated units and hence only these units should be used in estimation.In other cases, for example if treatment is assigned randomly, then it's reasonable to assume that both not-yet-and never-treated units would all satisfy parallel trends.
For estimators labeled "Not-yet-and/or never-treated", the default is to use both not-yet-and never-treated units in the estimator.However, if all never-treated units are dropped from the data set before using the estimator, then these estimators will use only not-yet-treated groups as the comparison group.did provides an option to use either the not-yet-treated or the never-treated group as a comparison group depending on which group a researcher thinks will make a better comparison group.staggered will automatically drop units that are never treated from the sample and hence only use not-yet-treated groups as a comparison group.
The next column, Main Assumptions, summarize concisely the main theoretical assumptions underlying each estimator.First, the assumptions about parallel trends match the previous discussion on the correct comparison group.The only estimator that doesn't rely on a parallel trends assumption is staggered which relies on the assumption that when a unit receives treatment is random.
The next assumption, that is common across all estimators, is that there should be "limited anticipation" of treatment.In general, anticipatory effects are when units respond to treatment before it is actually implemented.For example, this can be common if the news of a policy triggers behavior responses before the treatment is put in place."Limited anticipation" is when these anticipatory effects can only exist in a "few" pre-periods. 6In any of these cases, "treatment" should be manually moved back by the maximum number of periods where anticipation can occur.For example, if treatment starts in 2012 and anticipatory effects are reasonably only possible 2 years before, this units' "group" should be labeled as 2010 in the data.
The imputation-based estimators require an additional assumption that the parametric model of Y(0) = µ i + η t + ε it is correctly specified.This is because in the first stage, you have to accurately impute Y(0) when residualizing Y which relies on the correct specification of Y(0).The 2x2 aggregation models do not estimate a parametric form of Y(0) and hence only relies on a parallel trends assumption.The last column highlights that did allows for uniform inference of estimates.This addresses the problem that multiple hypotheses tests are being done by researchers (e.g.checking individually if all post periods are significant) by creating standard errors that adjust for multiple testing.
Example usage of event_study The result of event_study is a tibble in a tidy format (Robinson et al., 2021) that contains point estimates and standard errors for each relative time indicator for each estimator.The results of event_study are stored as a dataframe with event-study term, the estimate, standard error, and a column containing which estimator is used for that estimate.This output dataframe will in turn be passed to plot_event_study for easy comparison.plot_event_study will return a ggplot object (Wickham, 2016).We return to the df_het dataset to see example usage of these functions.

Conclusion
This article introduced the package did2s which provides a fast, memory-efficient, and treatment-effect heterogeneity robust way to estimate two-way fixed-effect models.The package also includes the event_study and plot_event_study functions to allow for a single syntax for the various estimators The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859

Figure 1 :
Figure 1: Example data with heterogeneous and dynamic treatment effects.Each line represents the average outcome in a given year for each group.In the absence of treatment, all three groups would exhibit parallel trends.

Figure 2 :
Figure 2: Event-study Estimate of TWFE Model.Standard Errors clustered at unit level.Estimated using the Two-Stage Difference-in-Differences proposed by Gardner (2021).

Figure 3 :
Figure 3: Event-study Estimate of TWFE Model.Standard Errors clustered at unit level.Estimated using the Two-Stage Difference-in-Differences proposed by Gardner (2021) and a Traditional TWFE model.
is random • Limited anticipation * This table summarizes the differences between various proposed event-study estimators in the econometric literature.* Anticipation can be accounted for by adjusting 'initial treatment day' back x periods, where x is the number of periods before treatment that anticipation can occur.The R Journal Vol.XX/YY, AAAA 20ZZ ISSN 2073-4859

Figure 4 :
Figure 4: Event-study Estimators.This figure plots the results from the event_study command on example data.
Since the returning object is a fixest object, all the accompanying output commands from fixest are available to use.For example, we can create regression tables:

Table 1 :
Estimate of Static TWFE Model Standard errors clustered at unit level.Estimated using Two-Stage Difference-in-Differences. proposed by Gardner (2021).