Examples with simulated data
Let’s start with a really simple example with simulated data. Here, there are going to be 4 time periods. There are 4000 units in the treated group that are randomly (with equal probability) assigned to first participate in the treatment (a group) in each time period. And there are 4000 ``never treated’’ units. The data generating process for untreated potential outcomes
\[ Y_{it}(0) = \theta_t + \eta_i + X_i'\beta_t + v_{it} \]
This is an example of a very simple model for untreated potential outcomes that is compatible with a conditional parallel trends assumption. In particular,
We consider the case where \(\eta_i\) can be distributed differently across groups. This means that comparisons of outcomes in levels between treated and untreated units will not deliver an average treatment effect parameter.
Next, notice that \[ \Delta Y_{it}(0) = (\theta_t - \theta_{t-1}) + X_i'(\beta_t - \beta_{t-1}) + \Delta v_{it} \]
so that the path of outcomes depends on covariates. And, in general, unconditional parallel trends is not valid in this setup unless either (i) the mean of the covariates is the same across groups (e.g., when treatment groups are independent of covariates), or (ii) \(\beta_t = \beta_{t-1} = \cdots = \beta_{1}\) (this is the case when the path of untreated potential outcomes doesn’t actually depend on covariates).
In order to think about treatment effects, we also use the following stylized model for treated potential outcomes \[ Y_{it}(g) = Y_{it}(0) + \mathbf{1}\{t \geq g\} (e+1) + ( u_{it} - v_{it}) \]
where \(e := t-g\) is the event time (i.e., the difference between the current time period and the time when a unit becomes treated), and the last term just allows for the error terms to be different for treated and untreated potential outcomes. It immediately follows that, in this particular example, \[ATT(g,t) = e+1\] for all post-treatment periods \(t \ge g\).
In other words, the average effect of participating in the treatment for units in group \(g\) is equal to their length of exposure plus one. Note that, for simplicity, we are considering the case where treatment effects are homogeneous across groups. Also, in this setup, units do not anticipate their treatment status, so the no-anticipation assumption is also satisfied.
For the simulations, we also set
- \(\theta_t = \beta_t = t\) for \(t=1,\ldots,4\)
- \(\eta_i \sim N(G_i, 1)\) where \(G_i\) is which group an individual belongs to,
- \(X_i \sim N( \mu_{D_i}, 1)\) where \(\mu_{D_i} = 1\) for never treated units and 0 otherwise,
- \(v_{it} \sim N(0,1)\), and \(u_{it} \sim N(0,1)\).
- Also, all random variables are drawn independently of each other.
- \(\theta_t = \beta_t = t\) for \(t=1,\ldots,4\)
Estimating Group-Time Average Treatment Effects
Building the dataset
# set seed so everything is reproducible
set.seed(1814)
# generate dataset with 4 time periods
time.periods <- 4
# add dynamic effects
te.e <- 1:time.periods
# generate data set with these parameters
# here, we dropped all units who are treated in time period 1 as they do not help us recover ATT(g,t)'s.
dta <- build_sim_dataset()
# How many observations remained after dropping the ``always-treated'' units
nrow(dta)
#> [1] 27952
#This is what the data looks like
head(dta)
#> G X id period Y treat
#> 1 2 0.530541 1 1 4.795125 1
#> 8001 2 0.530541 1 2 7.159789 1
#> 16001 2 0.530541 1 3 8.482624 1
#> 24001 2 0.530541 1 4 11.288278 1
#> 2 4 1.123550 2 1 4.375698 1
#> 8002 4 1.123550 2 2 9.285481 1Estimating Group-Time Average Treatment Effects
The main function to estimate group-time average treatment effects is the att_gt function. See documentation here. The most basic call to att_gt is given in the following
# estimate group-time average treatment effects using att_gt method
example_attgt <- att_gt(yname = "Y",
tname = "period",
idname = "id",
gname = "G",
xformla = ~X,
data = dta
)
# summarize the results
summary(example_attgt)
#>
#> Call:
#> att_gt(yname = "Y", tname = "period", idname = "id", gname = "G",
#> xformla = ~X, data = dta)
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Forthcoming at the Journal of Econometrics <https://arxiv.org/abs/1803.09015>, 2020.
#>
#> Group-Time Average Treatment Effects:
#> Group Time ATT(g,t) Std. Error [95% Simult. Conf. Band]
#> 2 2 0.9615 0.0593 0.8035 1.1196 *
#> 2 3 2.0424 0.0587 1.8859 2.1989 *
#> 2 4 3.0104 0.0601 2.8503 3.1705 *
#> 3 2 -0.0392 0.0551 -0.1862 0.1077
#> 3 3 1.0877 0.0617 0.9233 1.2522 *
#> 3 4 1.9938 0.0589 1.8369 2.1506 *
#> 4 2 -0.0537 0.0585 -0.2096 0.1022
#> 4 3 0.0382 0.0563 -0.1119 0.1883
#> 4 4 0.9410 0.0578 0.7869 1.0951 *
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> P-value for pre-test of parallel trends assumption: 0.79243
#> Control Group: Never Treated, Anticipation Periods: 0
#> Estimation Method: Doubly RobustThe summary of example_attgt provides estimates of the group-time average treatment effects in the column labeled att, and the corresponding bootstrapped-based standard errors are given in the column se. The corresponding groups and times are given in the columns group and time. Under the no-anticipation and parallel trends assumptions, group-time average treatment effects are identified in periods when \(t \geq g\) (i.e., post-treatment periods for each group). The table also reports pseudo group-time average treatment effects when \(t < g\) (i.e., pre-treatment periods for group \(g\)). These can be used as a pre-test for the parallel trends assumption (as long as we assume that the no-anticipation assumption indeed holds). In addition, the results of a Wald pre-test of the parallel trends assumption is reported in the summary of the results. A much more detailed discussion of using the did package for pre-testing is available here.
Next, we’ll demonstrate how to plot group-time average treatment effects. To plot these, use the ggdid function which builds off the ggplot2 package.
The resulting figure is one that contains separate plots for each group. Notice in the figure above, the first plot is labeled “Group 2”, the second “Group 3”, etc. Then, the figure contains estimates of group-time average treatment effects for each group in each time period along with a simultaneous confidence interval. The red dots in the plots are pre-treatment pseudo group-time average treatment effects and are most useful for pre-testing the parallel trends assumption. The blue dots are post-treatment group-time average treatment effects and should be interpreted as the average effect of participating in the treatment for units in a particular group at a particular point in time.
Other features of the did package
The above discussion covered only the most basic case for using the did package. There are a number of simple extensions that are useful in applications.
Adjustments for Multiple Hypothesis Testing
By default, the did package reports simultaneous confidence bands in plots of group-time average treatment effects with multiple time periods – these are confidence bands that are robust to multiple hypothesis testing [essentially, the idea here is to use the same standard errors but make an adjustment to the critical value to account for multiple testing – in the example in this section, the critical value for a 95% uniform confidence band is 2.66 instead of 1.96]. You can turn this off and compute analytical standard errors and corresponding figures with pointwise confidence intervals by setting bstrap=FALSE, cband=FALSE in the call to att_gt…but we don’t recommend it!
Aggregating group-time average treatment effects
In many applications, there can be a large number of groups and time periods. In this case, it may be infeasible to interpret plots of group-time average treatment effects. The did package provides a number of ways to aggregate group-time average treatment effects using the aggte function.
Simple Aggregation
One idea that is likely to immediately come to mind is to just return a weighted average of all group-time average treatment effects with weights proportional to the group size. This is available by calling the aggte function with type = simple.
agg.simple <- aggte(example_attgt, type = "simple")
summary(agg.simple)
#>
#> Call:
#> aggte(MP = example_attgt, type = "simple")
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Forthcoming at the Journal of Econometrics <https://arxiv.org/abs/1803.09015>, 2020.
#>
#>
#> Overall ATT:
#> ATT Std. Error [95% Conf. Int.]
#> 1.6557 0.0317 1.5937 1.7178 *
#>
#>
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> Control Group: Never Treated, Anticipation Periods: 0
#> Estimation Method: Doubly RobustThis sort of aggregation immediately avoids the negative weights issue that two-way fixed effects regressions can suffer from, but we often think that there are better alternatives. In particular, this simple aggregation tends to overweight the effect of early-treated groups simply because we observe more of them during post-treatment periods. We think there are likely to be better alternatives in most applications.
Dynamic Effects and Event Studies
One of the most common alternative approaches is to aggregate group-time effects into an event study plot. Group-time average treatment effects can immediately be averaged into average treatment effects at different lengths of exposure to the treatment using the following code:
agg.es <- aggte(example_attgt, type = "dynamic")
summary(agg.es)
#>
#> Call:
#> aggte(MP = example_attgt, type = "dynamic")
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Forthcoming at the Journal of Econometrics <https://arxiv.org/abs/1803.09015>, 2020.
#>
#>
#> Overall ATT:
#> ATT Std. Error [95% Conf. Int.]
#> 2.0086 0.0346 1.9408 2.0763 *
#>
#>
#> Dynamic Effects:
#> event time ATT Std. Error [95% Simult. Conf. Band]
#> -2 -0.0537 0.0593 -0.2031 0.0956
#> -1 -0.0011 0.0370 -0.0944 0.0923
#> 0 0.9988 0.0300 0.9231 1.0744 *
#> 1 2.0166 0.0403 1.9150 2.1181 *
#> 2 3.0104 0.0599 2.8594 3.1613 *
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> Control Group: Never Treated, Anticipation Periods: 0
#> Estimation Method: Doubly Robust
ggdid(agg.es)
In this figure, the x-axis is the length of exposure to the treatment. Length of exposure equal to 0 provides the average effect of participating in the treatment across groups in the time period when they first participate in the treatment (instantaneous treatment effect). Length of exposure equal to -1 corresponds to the time period before groups first participate in the treatment, and length of exposure equal to 1 corresponds to the first time period after initial exposure to the treatment.
As we would expect based on the data that we generated, it looks like parallel trends holds in pre-treatment periods and the effect of participating in the treatment is increasing with length of exposure of the treatment.
The Overall ATT here averages the average treatment effects across all lengths of exposure to the treatment.
Group-Specific Effects
Another idea is to look at average effects specific to each group. It is also straightforward to aggregate group-time average treatment effects into group-specific average treatment effects using the following code:
agg.gs <- aggte(example_attgt, type = "group")
summary(agg.gs)
#>
#> Call:
#> aggte(MP = example_attgt, type = "group")
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Forthcoming at the Journal of Econometrics <https://arxiv.org/abs/1803.09015>, 2020.
#>
#>
#> Overall ATT:
#> ATT Std. Error [95% Conf. Int.]
#> 1.4801 0.0293 1.4227 1.5375 *
#>
#>
#> Group Effects:
#> group ATT Std. Error [95% Simult. Conf. Band]
#> 2 2.0048 0.0493 1.8915 2.1181 *
#> 3 1.5407 0.0483 1.4297 1.6517 *
#> 4 0.9410 0.0580 0.8077 1.0743 *
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> Control Group: Never Treated, Anticipation Periods: 0
#> Estimation Method: Doubly Robust
ggdid(agg.gs)
In this figure, the y-axis is categorized by group. The x-axis provides estimates of the average effect of participating in the treatment for units in each group averaged across all time periods after that group becomes treated. In our example, the average effect of participating in the treatment is largest for group 2 – this is because the effect of treatment here increases with length of exposure to the treatment, and they are the group that is exposed to the treatment for the longest.
The Overall ATT averages the group-specific treatment effects across groups. In our view, this parameter is a leading choice as an overall summary effect of participating in the treatment. It is the average effect of participating in the treatment that was experienced across all units that participate in the treatment in any period. In this sense, it has a similar interpretation to the ATT in the textbook case where there are exactly two periods and two groups.
Calendar Time Effects
Finally, the did package allows aggregations across different time periods. To do this
agg.ct <- aggte(example_attgt, type = "calendar")
summary(agg.ct)
#>
#> Call:
#> aggte(MP = example_attgt, type = "calendar")
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Forthcoming at the Journal of Econometrics <https://arxiv.org/abs/1803.09015>, 2020.
#>
#>
#> Overall ATT:
#> ATT Std. Error [95% Conf. Int.]
#> 1.5028 0.0336 1.437 1.5685 *
#>
#>
#> Time Effects:
#> time ATT Std. Error [95% Simult. Conf. Band]
#> 2 0.9615 0.0591 0.8163 1.1067 *
#> 3 1.5651 0.0475 1.4482 1.6819 *
#> 4 1.9817 0.0407 1.8817 2.0817 *
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> Control Group: Never Treated, Anticipation Periods: 0
#> Estimation Method: Doubly Robust
ggdid(agg.ct)
In this figure, the x-axis is the time period and the estimates along the y-axis are the average effect of participating in the treatment in a particular time period for all groups that participated in the treatment in that time period.
Small Group Sizes
Small group sizes can sometimes cause estimation problems in the did package. To give an example, if there are any groups that have fewer observations than the number of covariates in the model, the code will error. The did package reports a warning if there are any groups that have fewer than 5 observations.
In addition, statistical inference, particularly on group-time average treatment may become more tenuous with small groups. For example, the effective sample size for estimating the change in outcomes over time for individuals in a particular group is equal to the number of observations in that group and asymptotic results are unlikely to provide good approximations to the sampling distribution of group-time average treatment effects when the number of units in a group is small. In these cases, one should be very cautious about interpreting the \(ATT(g,t)\) results for these small groups.
A reasonable alternative approach in this case is to just focus on aggregated treatment effect parameters (i.e., to run aggte(..., type = "group") or aggte(...,type = "dynamic")). For each of these cases, the effective sample size is the total number of units that are ever treated. As long as the total number of ever treated units is “large” (which should be the case for many DiD application), then the statistical inference results provided by the did package should be more stable.
Selecting Alternative Control Groups
By default, the did package uses the group of units that never participate in the treatment as the control group. In this case, if there is no group that never participates, then the did package will drop the last period and set units that do not become treated until the last period as the control group (this will also throw a warning). The other option for the control group is to use the “not yet treated”. The “not yet treated” include the never treated as well as those units that, for a particular point in time, have not been treated yet (though they eventually become treated). This group is at least as large as the never treated group though it changes across time periods. To use the “not yet treated” as the control, set the option control_group="notyettreated".
Repeated cross sections
The did package can also work with repeated cross section rather than panel data. If the data is repeated cross sections, simply set the option panel = FALSE. In this case, idname is also ignored. From a usage standpoint, everything else is identical to the case with panel data.
Unbalanced Panel Data
By default, the did package takes in panel data and, if it is not balanced, coerces it into being a balanced panel by dropping units with observations that are missing in any time period. However, if the user specifies the option allow_unbalanced_panel = TRUE, then the did package will not coerce the data into being balanced. Practically, the main cost here is that the computation time will increase in this case. We also recommend that users think carefully about why they have an unbalanced panel before proceeding this direction.
Alternative Estimation Methods
The did package implements all the \(2 \times 2\) DiD estimators that are in the DRDID package. By default, the did package uses “doubly robust” estimators that are based on first step linear regressions for the outcome variable and logit for the generalized propensity score. The other options are “ipw” for inverse probability weighting and “reg” for regression.
example_attgt_reg <- att_gt(yname = "Y",
tname = "period",
idname = "id",
gname = "G",
xformla = ~X,
data = data,
est_method = "reg"
)
summary(example_attgt_reg)The argument est_method is also available to pass in a custom function for estimating DiD with 2 periods and 2 groups. See its documentation for more details. Fair Warning: this is very advanced use of the did package and should be done with caution.
