Tutorial for Causal Panel Analysis
Ziyi Liu and Yiqing Xu
In recent years, researchers have proposed various robust estimators for heterogeneous treatment effects (HTE) in causal panel analysis as an alternative to traditional two-way fixed effects (TWFE) models. Examples of these estimators include those proposed by Sun and Abraham (2021), Callaway and Sant’Anna (2021), Liu, Wang, and Xu (2022), Borusyak, Jaravel, and Spiess (2021), and Imai, Kim, and Wang (2021). This tutorial will guide you through the implementation of these HTE-robust estimators, as well as TWFE, using R. It will also provide instructions on how to create event study plots that display estimated dynamic treatment effects. Additionally, this tutorial will present a recommended pipeline for analyzing panel data, covering data exploration, estimation, result visualization, and diagnostic tests. To begin, you will need to install necessary packages from CRAN and Github. The tutorial will then provide two empirical examples, based on Hainmueller and Hangartner (2019) and Grumbach and Sahn (2020), to illustrate these methods.
Installing Packages
# install packages from CRAN
packages <- c("dplyr", "readstata13", "fixest", "did", "fect",
"panelView", "PanelMatch", "ggplot2", "bacondecomp")
install.packages(setdiff(packages, rownames(installed.packages())))
# install "paneltools"
if ("paneltools" %in% rownames(installed.packages()) == FALSE) {
devtools:: install_github("xuyiqing/paneltools")
}library(dplyr)
library(readstata13)
library(fixest)
library(did)
library(fect)
library(panelView)
library(PanelMatch)
library(ggplot2)
library(bacondecomp)
library(paneltools)Example 1: w/o Treatment Reserveal
We begin with an empirical example from Hainmueller and Hopkins (2019), in which the authors investigate the effects of indirect democracy (versus direct democracy) on naturalization rates in Switzerland using municipality-year level panel data from 1991 to 2009. The study shows that a switch from direct to indirect democracy increased naturalization rates by an average of 1.22 percentage points (Model 1 of Table 1 in the paper).
data <- read.dta13(paste0(path, "jens_ajps2019_clean.dta"))
head(data)
#> bfs year nat_rate_ord indirect
#> 1 1 1991 0.000000 0
#> 2 1 1992 0.000000 0
#> 3 1 1993 0.000000 0
#> 4 1 1994 3.448276 0
#> 5 1 1995 0.000000 0
#> 6 1 1996 0.000000 0Data Visualization
First, we will examine the evolution of treatment status using the
panelView package. The variables bfs and
year are the unit and time indicators, respectively. The
treatment variable is indirect, and the outcome is
nat_rate_ord. In the baseline analysis, we will not include
any covariates. We will set the option by.timing = TRUE to
sort units by the timing of receiving the treatment.
panelview(nat_rate_ord ~ indirect, data = data, index = c("bfs","year"),
xlab = "Year", ylab = "Unit", display.all = T,
gridOff = TRUE, by.timing = TRUE)
The treatment status plot shows that the treatment rolls out in different years for different municipalities and there are no treatment reversals (staggered adoption). Following the definition in Sun and Abraham (2021), we define the units that take the treatment in the same year as a cohort. We can plot the average outcomes for each cohort using panelView.
panelview(data = data,Y='nat_rate_ord',
D='indirect',index=c("bfs","year"),
by.timing = T, display.all = T,
type = "outcome",by.cohort = T)
#> Number of unique treatment history: 17
TWFE
The staggered adoption setup allows us to implement several estimators to obtain the treatment effect estimates. We will first estimate a two-way fixed-effects (TWFE) model, as the authors do in the paper: \[Y_{it} = \alpha_i + \xi_t + \delta^{TWFE}D_{it} + \epsilon\] We can implement this estimator using the command feols in the **fixest package. The estimated coefficient using TWFE is 1.339, with a standard error of 0.187.
model.twfe.0 <- feols(nat_rate_ord~indirect|bfs+year,
data=data, cluster = "bfs") #use the clustered standard error
print(model.twfe.0)
#> OLS estimation, Dep. Var.: nat_rate_ord
#> Observations: 22,971
#> Fixed-effects: bfs: 1,209, year: 19
#> Standard-errors: Clustered (bfs)
#> Estimate Std. Error t value Pr(>|t|)
#> indirect 1.33932 0.186525 7.18039 1.2117e-12 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> RMSE: 4.09541 Adj. R2: 0.152719
#> Within R2: 0.005173Goodman-Bacon Decomposition
Goodman-Bacon (2021) demonstrates that
the two-way fixed-effects (TWFE) estimator in a staggered adoption
setting can be represented as a weighted average of all possible 2x2
difference-in-differences (DID) estimates between different cohorts.
However, when treatment effects change over time heterogeneously across
cohorts, the “forbidden” comparisons that use post-treatment data from
early adopters as controls for late adopters may introduce bias in the
TWFE estimator. We employ the procedure outlined in Goodman-Bacon (2021) and decompose the TWFE
estimate using the command bacon in the
bacondecomp package.
data.complete <- data[which(!is.na(data$nat_rate_ord)),] # bacon requires no missingness in the data
df_bacon <- bacon(nat_rate_ord~indirect,
data = data.complete,
id_var = "bfs",
time_var = "year")
#> type weight avg_est
#> 1 Earlier vs Later Treated 0.17605 1.97771
#> 2 Later vs Always Treated 0.31446 0.75233
#> 3 Later vs Earlier Treated 0.05170 0.32310
#> 4 Treated vs Untreated 0.45779 1.61180
ggplot(df_bacon) +
aes(x = weight, y = estimate, shape = factor(type), color = factor(type)) +
labs(x = "Weight", y = "Estimate", shape = "Type", color = 'Type') +
geom_point()
print(aggregate(df_bacon$estimate * df_bacon$weight, list(df_bacon$type), FUN=sum))
#> Group.1 x
#> 1 Earlier vs Later Treated 0.34817362
#> 2 Later vs Always Treated 0.23657303
#> 3 Later vs Earlier Treated 0.01670477
#> 4 Treated vs Untreated 0.73787343The decomposition of the TWFE estimator in a staggered adoption setting by Goodman-Bacon (2021) shows that the estimates from the DIDs comparing ever-treated cohorts switching into treatment with the never-treated (purple crosses labeled “Treated vs Untreated”) contribute the most to the TWFE estimate. The DIDs comparing ever-treated cohorts switching into treatment and other ever-treated cohorts that are still in their pretreatment periods (the red dots labeled “Earlier vs Later Treated”) and the “forbidden” comparisons between ever-treated cohorts and the always-treated (green triangles labeled “Later vs Always Treated”) rank second and third in terms of their contribution. The “forbidden” DIDs comparing ever-treated cohorts switching into treatment and other ever-treated cohorts that are already treated (the blue squares labeled “Later vs Earlier Treated”) contribute the least to the TWFE estimate.
Additionally, to further address the issue of heterogeneous treatment effects across cohorts, we can use methods specifically designed to handle heterogeneous treatment effects, such as the methods outlined in Sun and Abraham (2021), Callaway and Sant’Anna (2021), and Liu, Wang, and Xu (2022). These methods can provide robust causal estimates that take into account the heterogeneous nature of treatment effects across different cohorts. By using these methods, we can better account for the potential bias introduced by the “forbidden” comparisons and obtain more accurate estimates of the treatment effects.
# drop always treated units
df <- as.data.frame(data %>% group_by(bfs) %>% mutate(treatment_mean = mean(indirect,na.rm = TRUE)))
df.use <- df[which(df$treatment_mean<1),]
# Re-estimate TWFE on this Sub-sample
model.twfe.1 <- feols(nat_rate_ord~indirect|bfs+year,
data=df.use, cluster = "bfs")
print(model.twfe.1)
#> OLS estimation, Dep. Var.: nat_rate_ord
#> Observations: 17,594
#> Fixed-effects: bfs: 926, year: 19
#> Standard-errors: Clustered (bfs)
#> Estimate Std. Error t value Pr(>|t|)
#> indirect 1.60858 0.195296 8.23661 6.009e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> RMSE: 4.43613 Adj. R2: 0.143692
#> Within R2: 0.007415After removing the always-treated units, the estimated average treatment effect (ATT) increases to 1.609, with a standard error of 0.195. However, to fully address the bias caused by heterogeneous treatment effects (HTE), it is necessary to eliminate the “forbidden comparisons” through the use of HTE-robust estimators. Before implementing these estimators, we will demonstrate how to create an event study plot using a TWFE model with a dynamic specification.
Event Study Plot for TWFE
In order to create an event study plot, we must first generate the
cohort index and the relative time period for each observation in
relation to treatment. These details are essential for various
HTE-robust estimators, such as Sun and Abraham
(2021) and Callaway and Sant’Anna
(2021). We can use the function get.cohort in the
package paneltools to generate these indices.
df.use <- get.cohort(df.use,D="indirect",index=c("bfs","year"))
head(df.use[,-5],19)
#> bfs year nat_rate_ord indirect FirstTreat Cohort Time_to_Treatment
#> 1 1 1991 0.000000 0 2006 Cohort:2006 -15
#> 2 1 1992 0.000000 0 2006 Cohort:2006 -14
#> 3 1 1993 0.000000 0 2006 Cohort:2006 -13
#> 4 1 1994 3.448276 0 2006 Cohort:2006 -12
#> 5 1 1995 0.000000 0 2006 Cohort:2006 -11
#> 6 1 1996 0.000000 0 2006 Cohort:2006 -10
#> 7 1 1997 0.000000 0 2006 Cohort:2006 -9
#> 8 1 1998 0.000000 0 2006 Cohort:2006 -8
#> 9 1 1999 0.000000 0 2006 Cohort:2006 -7
#> 10 1 2000 0.000000 0 2006 Cohort:2006 -6
#> 11 1 2001 2.380952 0 2006 Cohort:2006 -5
#> 12 1 2002 0.000000 0 2006 Cohort:2006 -4
#> 13 1 2003 0.000000 0 2006 Cohort:2006 -3
#> 14 1 2004 0.000000 0 2006 Cohort:2006 -2
#> 15 1 2005 0.000000 0 2006 Cohort:2006 -1
#> 16 1 2006 1.612903 1 2006 Cohort:2006 0
#> 17 1 2007 10.937500 1 2006 Cohort:2006 1
#> 18 1 2008 14.285715 1 2006 Cohort:2006 2
#> 19 1 2009 9.677419 1 2006 Cohort:2006 3The variable FirstTreat indicates the first period when
a unit gets treated (for never-treated units, FirstTreat is
set to NA), and the variable Time_to_Treatment
indicates the relative period to the treatment’s onset (for
always-treated and never-treated units, Time_to_Treatment
is set to NA). The period -1 corresponds to
the last pre-treatment period.
We first use the dynamic TWFE regression, which is the workhorse model to estimate dynamic treatment effects in event studies, to create the event study plot. This model includes a series of interaction terms between a dummy that indicates whether a unit is a treated unit and each lead (lag) indicator relative to the treatment in a TWFE regression. This specification allows for effects to vary over time and usually uses the period immediately preceding a switch into the treatment as the reference period. If the regression is saturated and there is no heterogeneous treatment effects across cohorts, this dynamic TWFE can consistently estimate the dynamic treatment effect.
We still use the function feols to implement this
estimator. The variable Time_to_Treatment can serve as the
lead (lag) indicator (we need to replace the NA in
Time_to_Treatment with an arbitrary number). Note that the
period -1 in Time_to_Treatment corresponds to
the last pre-treatment period and is the reference period.
# TWFE Dynamic
df.twfe <- df.use
df.twfe$treat <- as.numeric(df.twfe$treatment_mean>0) # drop always treated units
df.twfe[which(is.na(df.twfe$Time_to_Treatment)),'Time_to_Treatment'] <- 0 # can be an arbitrary value
twfe.est <- feols(nat_rate_ord ~ i(Time_to_Treatment, treat, ref = -1)| bfs + year,
data = df.twfe, cluster = "bfs")
twfe.output <- as.matrix(twfe.est$coeftable)
print(twfe.output)
#> Estimate Std. Error t value Pr(>|t|)
#> Time_to_Treatment::-18:treat -0.32065076 0.6327017 -0.5067961 6.124187e-01
#> Time_to_Treatment::-17:treat -0.11854037 0.3911477 -0.3030578 7.619139e-01
#> Time_to_Treatment::-16:treat -0.39563213 0.4029210 -0.9819098 3.264010e-01
#> Time_to_Treatment::-15:treat -0.41097577 0.3438989 -1.1950482 2.323745e-01
#> Time_to_Treatment::-14:treat 0.24451225 0.3691273 0.6624062 5.078759e-01
#> Time_to_Treatment::-13:treat -0.04431012 0.3540169 -0.1251639 9.004210e-01
#> Time_to_Treatment::-12:treat -0.54589060 0.3151718 -1.7320415 8.359954e-02
#> Time_to_Treatment::-11:treat -0.29313843 0.3468354 -0.8451804 3.982287e-01
#> Time_to_Treatment::-10:treat -0.42398107 0.3193748 -1.3275344 1.846595e-01
#> Time_to_Treatment::-9:treat -0.46711352 0.3010171 -1.5517839 1.210560e-01
#> Time_to_Treatment::-8:treat -0.11563334 0.3336591 -0.3465614 7.289997e-01
#> Time_to_Treatment::-7:treat -0.72428992 0.3057269 -2.3690752 1.803716e-02
#> Time_to_Treatment::-6:treat -0.57295020 0.3171533 -1.8065400 7.115891e-02
#> Time_to_Treatment::-5:treat -0.55499321 0.3092458 -1.7946670 7.303318e-02
#> Time_to_Treatment::-4:treat 0.06620907 0.3451034 0.1918529 8.478996e-01
#> Time_to_Treatment::-3:treat -0.50953012 0.3162958 -1.6109292 1.075363e-01
#> Time_to_Treatment::-2:treat -0.19509150 0.3302199 -0.5907926 5.548037e-01
#> Time_to_Treatment::0:treat 0.52271487 0.3596460 1.4534151 1.464477e-01
#> Time_to_Treatment::1:treat 1.66883434 0.3679186 4.5358788 6.490737e-06
#> Time_to_Treatment::2:treat 1.72423014 0.4713628 3.6579681 2.685814e-04
#> Time_to_Treatment::3:treat 1.17833624 0.4005586 2.9417322 3.345229e-03
#> Time_to_Treatment::4:treat 0.75632268 0.4786966 1.5799624 1.144573e-01
#> Time_to_Treatment::5:treat 2.43822935 0.7054512 3.4562694 5.726291e-04
#> Time_to_Treatment::6:treat 2.84859038 0.8053149 3.5372380 4.244244e-04
#> Time_to_Treatment::7:treat 2.40935895 0.7837161 3.0742752 2.172064e-03
#> Time_to_Treatment::8:treat 2.27042846 0.7978687 2.8456167 4.530459e-03
#> Time_to_Treatment::9:treat 2.39267761 0.8632219 2.7717990 5.686658e-03
#> Time_to_Treatment::10:treat 2.00283605 0.9536545 2.1001695 3.598421e-02
#> Time_to_Treatment::11:treat 0.84259034 0.5602476 1.5039607 1.329328e-01
#> Time_to_Treatment::12:treat 3.83191775 2.3004225 1.6657452 9.610276e-02
#> Time_to_Treatment::13:treat 4.43643668 2.3223177 1.9103487 5.639730e-02
#> Time_to_Treatment::14:treat 4.75697725 3.4855231 1.3647814 1.726536e-01
#> Time_to_Treatment::15:treat 1.54196720 1.4152835 1.0895112 2.762123e-01
#> Time_to_Treatment::16:treat 7.16894343 4.8440161 1.4799586 1.392248e-01
#> Time_to_Treatment::17:treat 6.49085407 2.7265262 2.3806315 1.748463e-02
#> attr(,"type")
#> [1] "Clustered (bfs)"The paneltools package offers a useful function,
esplot(), that can be utilized to effectively visualize
results. In this instance, we utilize the function to examine the
dynamic effects over a 13 period pre-treatment span and a 10 period
post-treatment span. The data reveals a significant impact of the
treatment, as well as a subtle pre-trend discrepancy in a number of
periods preceding the treatment. It is important to note that
esplot() considers period 0 to be the final pre-treatment
period, so it is necessary to shift the period index by one unit.
twfe.output <- as.data.frame(twfe.output)
twfe.output$Time <- c(c(-18:-2),c(0:17))+1 #
p.twfe <- esplot(twfe.output,Period = 'Time',Estimate = 'Estimate',
SE = 'Std. Error', xlim = c(-12,10))
p.twfe
Stacked DID
The study, Cengiz et al. (2019), examines the effect of variations in minimum wage on low-wage employment across 138 state-level minimum wage changes in the United States between 1979 and 2016, utilizing a stacked difference-in-differences (DID) approach. A similar estimator is also implemented in the Stata package STACKEDEV, as per Bleiberg (2021). This stacked DID method eliminates the bias resulted from “forbidden” comparisons in the Goodman-Bacon decomposition. However, it still relies on the implicit weights chosen by OLS, making it less robust in terms of heterogeneous treatment effects (HTE) and less interpretable compared to other HTE-robust estimators.
The study employs a stacked difference-in-differences approach to estimate the impact of minimum wage changes on low-wage jobs across a series of 138 state-level minimum wage changes between 1979 and 2016 in the United States. For each ever-treated cohort, a cohort-specific dataset is created, comprising of the cohort and all never-treated units. These cohort-specific datasets are then stacked in order to calculate an average effect across all cohorts by utilizing fixed effects regression, while controlling for stack-unit interaction fixed effects and stack-year interaction fixed effects. With this stacked cohort-specific dataset, dynamic treatment effects can be obtained in a similar way as with the TWFE estimator.
df.st <- NULL
target.cohorts <- setdiff(unique(df.use$Cohort),"Control")
k <- 1
for(cohort in target.cohorts){
df.sub <- df.use[which(df.use$Cohort%in%c(cohort,"Control")),]
df.sub$stack <- k
df.st <- rbind(df.st,df.sub)
k <- k + 1
}
df.st$st_unit <- as.numeric(factor(paste0(df.st$stack,'-',df.st$bfs)))
df.st$st_year <- as.numeric(factor(paste0(df.st$stack,'-',df.st$year)))
model.st <- feols(nat_rate_ord~indirect|st_unit+st_year,
data=df.st, cluster = "st_unit")
print(model.st)
#> OLS estimation, Dep. Var.: nat_rate_ord
#> Observations: 127,186
#> Fixed-effects: st_unit: 6,694, st_year: 285
#> Standard-errors: Clustered (st_unit)
#> Estimate Std. Error t value Pr(>|t|)
#> indirect 1.66861 0.18237 9.14958 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> RMSE: 3.9095 Adj. R2: 0.130512
#> Within R2: 0.001783Event Study Plot
df.st$treat <- as.numeric(df.st$treatment_mean>0)
df.st[which(is.na(df.st$Time_to_Treatment)),'Time_to_Treatment'] <- 0 # can be arbitrary value
st.est <- feols(nat_rate_ord ~ i(Time_to_Treatment, treat, ref = -1)| st_unit + st_year,data = df.st,cluster = "st_unit")
print(st.est$coeftable)
#> Estimate Std. Error t value Pr(>|t|)
#> Time_to_Treatment::-18:treat -0.32876531 0.6246666 -0.52630524 5.986936e-01
#> Time_to_Treatment::-17:treat -0.08402861 0.3726755 -0.22547391 8.216175e-01
#> Time_to_Treatment::-16:treat -0.23704732 0.3890202 -0.60934446 5.423168e-01
#> Time_to_Treatment::-15:treat -0.33039222 0.3358986 -0.98360689 3.253444e-01
#> Time_to_Treatment::-14:treat 0.21942491 0.3595513 0.61027425 5.417009e-01
#> Time_to_Treatment::-13:treat -0.13145445 0.3443222 -0.38177744 7.026386e-01
#> Time_to_Treatment::-12:treat -0.58822817 0.3113593 -1.88922655 5.890467e-02
#> Time_to_Treatment::-11:treat -0.26344604 0.3437730 -0.76633722 4.435027e-01
#> Time_to_Treatment::-10:treat -0.54375457 0.3218685 -1.68936845 9.119542e-02
#> Time_to_Treatment::-9:treat -0.56120024 0.3065223 -1.83086265 6.716548e-02
#> Time_to_Treatment::-8:treat -0.09265692 0.3389478 -0.27336634 7.845801e-01
#> Time_to_Treatment::-7:treat -0.58880594 0.3070364 -1.91770753 5.519062e-02
#> Time_to_Treatment::-6:treat -0.47516514 0.3257192 -1.45881816 1.446621e-01
#> Time_to_Treatment::-5:treat -0.39993315 0.3112587 -1.28488986 1.988752e-01
#> Time_to_Treatment::-4:treat 0.23387391 0.3587515 0.65191065 5.144812e-01
#> Time_to_Treatment::-3:treat -0.24046110 0.3308839 -0.72672341 4.674208e-01
#> Time_to_Treatment::-2:treat -0.01674545 0.3466380 -0.04830817 9.614721e-01
#> Time_to_Treatment::0:treat 0.63988908 0.3639778 1.75804440 7.878566e-02
#> Time_to_Treatment::1:treat 1.63436478 0.3698039 4.41954479 1.004695e-05
#> Time_to_Treatment::2:treat 1.78039699 0.4911171 3.62519828 2.908847e-04
#> Time_to_Treatment::3:treat 1.10329058 0.4123146 2.67584672 7.472221e-03
#> Time_to_Treatment::4:treat 0.85952703 0.4895732 1.75566599 7.919120e-02
#> Time_to_Treatment::5:treat 2.44794424 0.7113218 3.44140212 5.822459e-04
#> Time_to_Treatment::6:treat 2.80029784 0.8125259 3.44641044 5.715789e-04
#> Time_to_Treatment::7:treat 2.29340055 0.7978800 2.87436774 4.061131e-03
#> Time_to_Treatment::8:treat 2.38558084 0.7977891 2.99023997 2.797728e-03
#> Time_to_Treatment::9:treat 2.20907842 0.8509387 2.59604885 9.450873e-03
#> Time_to_Treatment::10:treat 1.90401019 0.9458353 2.01304633 4.414972e-02
#> Time_to_Treatment::11:treat 0.51499900 0.5518266 0.93326236 3.507182e-01
#> Time_to_Treatment::12:treat 3.95439028 2.2940011 1.72379617 8.479083e-02
#> Time_to_Treatment::13:treat 4.39078917 2.3089816 1.90161286 5.726468e-02
#> Time_to_Treatment::14:treat 4.61626419 3.4608264 1.33386181 1.822945e-01
#> Time_to_Treatment::15:treat 1.18977420 1.3883525 0.85696841 3.914930e-01
#> Time_to_Treatment::16:treat 6.75219353 4.8161984 1.40197579 1.609689e-01
#> Time_to_Treatment::17:treat 5.95993087 2.7217298 2.18975843 2.857616e-02
#> attr(,"type")
#> [1] "Clustered (st_unit)"st.output <- as.data.frame(st.est$coeftable)
st.output$Time <- c(c(-18:-2),c(0:17))+1 #
p.st <- esplot(st.output,Period = 'Time',Estimate = 'Estimate',
SE = 'Std. Error', xlim = c(-12,10))
p.st
Interaction Weighted
In the study Sun and Abraham (2021), an
interaction-weighted (iw) estimator is proposed as a means of estimating
heterogeneous treatment effects. This estimator is a weighted average of
the average treatment effect (ATT) estimates for each cohort, obtained
from a TWFE regression that includes cohort dummies fully interacted
with indicators of relative time to the treatment’s onset. The iw
estimator is robust to heterogeneous treatment effects (HTE) and can be
implemented using the command sunab in the package
fixest. The variable FirstTreat serves as
the cohort indicator, and any missing values in this variable should be
replaced with an arbitrary number. The estimated ATT obtained using the
iw estimator is 1.331, with a standard error of 0.288.
df.sa <- df.use
df.sa[which(is.na(df.sa$FirstTreat)),"FirstTreat"] <- 1000 # replace NA with an arbitrary number
model.sa.1 <- feols(nat_rate_ord~sunab(FirstTreat,year)|bfs+year,
data = df.sa, cluster = "bfs")
summary(model.sa.1,agg = "ATT")
#> OLS estimation, Dep. Var.: nat_rate_ord
#> Observations: 17,594
#> Fixed-effects: bfs: 926, year: 19
#> Standard-errors: Clustered (bfs)
#> Estimate Std. Error t value Pr(>|t|)
#> ATT 1.3309 0.287971 4.62163 4.3474e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> RMSE: 4.34795 Adj. R2: 0.163888
#> Within R2: 0.046484Event Study Plot
The results of this estimation are saved in the object
coeftable for further analysis.
sa.output <- as.matrix(model.sa.1$coeftable)
print(sa.output)
#> Estimate Std. Error t value Pr(>|t|)
#> year::-18 -0.26588766 0.7883091 -0.33728860 7.359759e-01
#> year::-17 0.17803769 0.5029562 0.35398250 7.234327e-01
#> year::-16 0.11405908 0.5185987 0.21993705 8.259687e-01
#> year::-15 -0.10395137 0.4074327 -0.25513753 7.986736e-01
#> year::-14 0.41028405 0.4131183 0.99313944 3.209017e-01
#> year::-13 0.09073535 0.3972375 0.22841588 8.193734e-01
#> year::-12 -0.56191688 0.3689700 -1.52293385 1.281170e-01
#> year::-11 -0.32332911 0.3836954 -0.84267130 3.996301e-01
#> year::-10 -0.50966995 0.3531793 -1.44309135 1.493333e-01
#> year::-9 -0.46891460 0.3436058 -1.36468774 1.726831e-01
#> year::-8 0.06132399 0.3645946 0.16819776 8.664645e-01
#> year::-7 -0.50027176 0.3084227 -1.62203273 1.051371e-01
#> year::-6 -0.39220920 0.3289499 -1.19230663 2.334468e-01
#> year::-5 -0.22991949 0.3217699 -0.71454628 4.750697e-01
#> year::-4 0.41857502 0.3512752 1.19158713 2.337288e-01
#> year::-3 -0.19712944 0.3199730 -0.61608158 5.379922e-01
#> year::-2 -0.02786173 0.3420857 -0.08144663 9.351044e-01
#> year::0 0.71412228 0.3755381 1.90159764 5.753426e-02
#> year::1 1.58464318 0.3615507 4.38290760 1.305219e-05
#> year::2 1.63374738 0.5093517 3.20750354 1.384925e-03
#> year::3 0.90899520 0.4276159 2.12572845 3.379026e-02
#> year::4 0.40908313 0.5464691 0.74859332 4.542927e-01
#> year::5 2.11063582 0.7741950 2.72623267 6.527297e-03
#> year::6 2.48424347 0.8628637 2.87906816 4.080461e-03
#> year::7 2.38015298 0.7861370 3.02765689 2.532923e-03
#> year::8 2.56023855 0.8504950 3.01029215 2.680822e-03
#> year::9 2.00995330 0.9579247 2.09823736 3.615491e-02
#> year::10 1.63910238 0.9721952 1.68598068 9.213676e-02
#> year::11 0.02976808 0.5346014 0.05568277 9.556066e-01
#> year::12 0.69417518 2.7317900 0.25411001 7.994670e-01
#> year::13 1.13548295 1.9972762 0.56851575 5.698228e-01
#> year::14 1.36246531 3.1273735 0.43565801 6.631865e-01
#> year::15 -2.06989292 1.4381238 -1.43930094 1.504036e-01
#> year::16 2.27432076 3.0784404 0.73878992 4.602220e-01
#> year::17 0.10722911 4.5738694 0.02344385 9.813013e-01
#> attr(,"type")
#> [1] "Clustered (bfs)"We can visualize the estimated dynamic treatment effects using the same method as we do for the TWFE estimator. The results indicate a significant treatment effect, but no pre-treatment trend difference is observed. The minimal difference between the results obtained from the TWFE and iw estimators suggests that the weak pre-treatment trend difference may be a result of heterogeneous treatment effects rather than a violation of the parallel trend assumption.
sa.output <- as.data.frame(sa.output)
sa.output$Time <- c(c(-18:-2),c(0:17))+1
p.sa <- esplot(sa.output,Period = 'Time',Estimate = 'Estimate',
SE = 'Std. Error', xlim = c(-12,10))
p.sa
CSDID
Callaway and Sant’Anna (2021) proposed a doubly-robust estimator, referred to as csdid, that incorporates pre-treatment covariates by using either never-treated or not-yet-treated units as the comparison group. In this tutorial, we will only be utilizing the outcome model instead of the double-robust model to calculate the ATT and dynamic treatment effects.
First, we use the never-treated units as the comparison group. The
variable FirstTreat serves as the cohort indicator, where
we replace the missing values in FirstTreat with 0. The
estimated ATT is numerically the same as the iw estimator. Additionally,
we set the option bstrap to FALSE in order to
report the analytical standard error, which is 0.302.
df.cs <- df.use
df.cs[which(is.na(df.cs$FirstTreat)),"FirstTreat"] <- 0 # replace NA with 0
cs.est.1 <- att_gt(yname = "nat_rate_ord",
gname = "FirstTreat",
idname = "bfs",
tname = "year",
xformla = ~1,
control_group = "nevertreated",
allow_unbalanced_panel = TRUE,
data = df.cs,
est_method = "reg")
cs.est.att.1 <- aggte(cs.est.1, type = "simple", na.rm=T, bstrap = F)
print(cs.est.att.1)
#>
#> Call:
#> aggte(MP = cs.est.1, type = "simple", na.rm = T, bstrap = F)
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015>
#>
#>
#> ATT Std. Error [ 95% Conf. Int.]
#> 1.3309 0.3019 0.7392 1.9226 *
#>
#>
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> Control Group: Never Treated, Anticipation Periods: 0
#> Estimation Method: Outcome RegressionEvent Study Plot
The csdid estimator utilizes the command aggte to
aggregate the cohort-period average treatment effects. We set the option
bstrap to FALSE to report analytical standard
errors and the option cband to FALSE to report
period-wise confidence intervals. This allows for easy comparison with
other estimators. Note that the aggte command also has the
capability to compute uniform confidence bands by setting
cband and bstrap to TRUE, but in
this case we chose to report period-wise intervals.
cs.att.1 <- aggte(cs.est.1, type = "dynamic",
bstrap=FALSE, cband=FALSE, na.rm=T)
print(cs.att.1)
#>
#> Call:
#> aggte(MP = cs.est.1, type = "dynamic", na.rm = T, bstrap = FALSE,
#> cband = FALSE)
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015>
#>
#>
#> Overall summary of ATT's based on event-study/dynamic aggregation:
#> ATT Std. Error [ 95% Conf. Int.]
#> 1.2205 0.6314 -0.0171 2.458
#>
#>
#> Dynamic Effects:
#> Event time Estimate Std. Error [95% Pointwise Conf. Band]
#> -17 0.1350 0.6467 -1.1324 1.4025
#> -16 -0.0839 0.3583 -0.7861 0.6183
#> -15 0.1467 0.3964 -0.6303 0.9236
#> -14 0.7509 0.3219 0.1200 1.3818 *
#> -13 -0.2318 0.3392 -0.8967 0.4331
#> -12 -0.5851 0.2814 -1.1367 -0.0335 *
#> -11 0.2656 0.2705 -0.2647 0.7958
#> -10 -0.2386 0.2805 -0.7883 0.3111
#> -9 0.0576 0.2702 -0.4721 0.5872
#> -8 0.5244 0.3165 -0.0959 1.1448
#> -7 -0.6282 0.3132 -1.2420 -0.0144 *
#> -6 0.1081 0.2565 -0.3946 0.6108
#> -5 0.1623 0.2890 -0.4041 0.7287
#> -4 0.6485 0.3380 -0.0139 1.3109
#> -3 -0.6161 0.3674 -1.3362 0.1040
#> -2 0.1505 0.3355 -0.5071 0.8081
#> -1 0.0279 0.3425 -0.6435 0.6992
#> 0 0.7141 0.3754 -0.0217 1.4500
#> 1 1.5846 0.3681 0.8632 2.3060 *
#> 2 1.6337 0.5118 0.6306 2.6369 *
#> 3 0.9090 0.4396 0.0473 1.7707 *
#> 4 0.4091 0.5725 -0.7130 1.5311
#> 5 2.1106 0.7976 0.5473 3.6739 *
#> 6 2.4842 0.9078 0.7051 4.2634 *
#> 7 2.3802 0.8538 0.7068 4.0535 *
#> 8 2.5602 0.8723 0.8506 4.2699 *
#> 9 2.0100 0.9769 0.0953 3.9246 *
#> 10 1.6391 1.0257 -0.3712 3.6494
#> 11 0.0298 0.6100 -1.1658 1.2253
#> 12 0.6942 2.8043 -4.8021 6.1904
#> 13 1.1355 2.3284 -3.4281 5.6991
#> 14 1.3625 3.4131 -5.3270 8.0519
#> 15 -2.0699 1.7160 -5.4331 1.2933
#> 16 2.2743 3.6182 -4.8172 9.3658
#> 17 0.1072 4.5338 -8.7789 8.9934
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> Control Group: Never Treated, Anticipation Periods: 0
#> Estimation Method: Outcome RegressionThe estimated dynamic treatment effects during the post-treatment period are numerically the same as those obtained from the iw estimator. However, the estimates for the pre-treatment period may differ as the csdid estimator does not rely on a single reference period. We present the results in a similar manner and it can be observed that there is a significant treatment effect present, with no clear pre-trend difference.
cs.output <- cbind.data.frame(Estimate = cs.att.1$att.egt,
SE = cs.att.1$se.egt,
time = cs.att.1$egt + 1)
p.cs.1 <- esplot(cs.output,Period = 'time',Estimate = 'Estimate',
SE = 'SE', xlim = c(-12,10))
p.cs.1
We can also use the not-yet-treated units as the
comparison group by setting
control_group = "notyettreated". The choice of the
comparison group depends on the specific version of the parallel trend
assumption that researchers need to justify. The estimated ATT is 1.292,
with a standard error of 0.308.
cs.est.2 <- att_gt(yname = "nat_rate_ord",
gname = "FirstTreat",
idname = "bfs",
tname = "year",
xformla = ~1,
control_group = "notyettreated",
allow_unbalanced_panel = TRUE,
data = df.cs,
est_method = "reg")
cs.est.att.2 <- aggte(cs.est.2, type = "simple",na.rm=T, bstrap = F)
print(cs.est.att.2)
#>
#> Call:
#> aggte(MP = cs.est.2, type = "simple", na.rm = T, bstrap = F)
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015>
#>
#>
#> ATT Std. Error [ 95% Conf. Int.]
#> 1.2923 0.308 0.6886 1.896 *
#>
#>
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> Control Group: Not Yet Treated, Anticipation Periods: 0
#> Estimation Method: Outcome RegressionWe can calculate the dynamic treatment effects using the same method.
cs.att.2 <- aggte(cs.est.2, type = "dynamic",
bstrap=FALSE,cband=FALSE,na.rm=T)
print(cs.att.2)
#>
#> Call:
#> aggte(MP = cs.est.2, type = "dynamic", na.rm = T, bstrap = FALSE,
#> cband = FALSE)
#>
#> Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015>
#>
#>
#> Overall summary of ATT's based on event-study/dynamic aggregation:
#> ATT Std. Error [ 95% Conf. Int.]
#> 1.2359 0.6337 -0.0062 2.478
#>
#>
#> Dynamic Effects:
#> Event time Estimate Std. Error [95% Pointwise Conf. Band]
#> -17 -0.0629 0.6363 -1.3100 1.1841
#> -16 -0.2656 0.3626 -0.9763 0.4452
#> -15 0.2242 0.4072 -0.5738 1.0222
#> -14 0.8526 0.3273 0.2111 1.4940 *
#> -13 -0.2011 0.3436 -0.8746 0.4724
#> -12 -0.6488 0.2829 -1.2033 -0.0943 *
#> -11 0.1963 0.2735 -0.3396 0.7323
#> -10 -0.0609 0.2733 -0.5966 0.4748
#> -9 0.0233 0.2663 -0.4985 0.5452
#> -8 0.3539 0.3225 -0.2781 0.9859
#> -7 -0.7367 0.3253 -1.3743 -0.0992 *
#> -6 0.1658 0.2615 -0.3468 0.6784
#> -5 0.0376 0.3036 -0.5574 0.6325
#> -4 0.7782 0.3527 0.0869 1.4695 *
#> -3 -0.7900 0.3849 -1.5444 -0.0356 *
#> -2 0.3214 0.3423 -0.3495 0.9923
#> -1 0.0906 0.3386 -0.5731 0.7543
#> 0 0.5757 0.3835 -0.1761 1.3274
#> 1 1.5347 0.3755 0.7986 2.2707 *
#> 2 1.5818 0.5134 0.5756 2.5880 *
#> 3 0.8832 0.4451 0.0108 1.7556 *
#> 4 0.3647 0.5792 -0.7705 1.5000
#> 5 2.1045 0.8005 0.5355 3.6734 *
#> 6 2.6091 0.9054 0.8346 4.3836 *
#> 7 2.4869 0.8481 0.8247 4.1492 *
#> 8 2.5860 0.8725 0.8759 4.2961 *
#> 9 2.1076 0.9814 0.1841 4.0311 *
#> 10 1.5353 1.0223 -0.4683 3.5389
#> 11 0.0751 0.6043 -1.1094 1.2596
#> 12 0.8129 2.8107 -4.6960 6.3218
#> 13 1.2552 2.3333 -3.3181 5.8284
#> 14 1.4144 3.4089 -5.2670 8.0958
#> 15 -2.0470 1.7125 -5.4034 1.3094
#> 16 2.2586 3.6230 -4.8425 9.3596
#> 17 0.1072 4.5338 -8.7789 8.9934
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#>
#> Control Group: Not Yet Treated, Anticipation Periods: 0
#> Estimation Method: Outcome Regressioncs.output <- cbind.data.frame(Estimate = cs.att.2$att.egt,
SE = cs.att.2$se.egt,
time = cs.att.2$egt + 1)
p.cs.2 <- esplot(cs.output,Period = 'time',Estimate = 'Estimate',
SE = 'SE', xlim = c(-12,10))
p.cs.2
PanelMatch
In Imai, Kim, and Wang (2021), a matching framework is proposed as a generalization of the difference-in-differences estimator. This framework allows researchers to match each treated observation for a given unit in a particular period with untreated observations from other units in the same period that have similar treatment, outcome, or covariate histories.
The PanelMatch command can be utilized to create these
matched sets. The option lag = 3 indicates that the
pre-treatment history periods to be matched on should have a length of
3, while the option lead = c(0:3) specifies that the lead
window should be four post-treatment periods. To assign equal weights to
all control units in each matched set, the option
refinement.method="none" is set. It is important to note
that by matching on treatment history and specifying the lead window,
the PanelMatch command only uses a subset of treated units
that have 3 pre-treatment periods and 4 post-treatment periods to
compute the average treatment effects.
The PanelMatch estimator is also equivalent to the \(DID_{M}\) estimator proposed by Chaisemartin and D’Haultfœuille (2020) if we only match on the last pre-treatment period and only estimate the treatment effect at the first post-treatment period.
df.pm <- df.use
# we need to convert the unit and time indicator to integer
df.pm[,"bfs"] <- as.integer(as.factor(df.pm[,"bfs"]))
df.pm[,"year"] <- as.integer(as.factor(df.pm[,"year"]))
df.pm <- df.pm[,c("bfs","year","nat_rate_ord","indirect")]
PM.results <- PanelMatch(lag=3,
time.id="year",
unit.id = "bfs",
treatment = 'indirect',
refinement.method = "none",
data = df.pm,
qoi = "att",
lead = c(0:3),
outcome.var = 'nat_rate_ord',
match.missing = TRUE)
## For pre-treatment dynamic effects
PM.results.placebo <- PanelMatch(lag=3,
time.id="year",
unit.id = "bfs",
treatment = 'indirect',
refinement.method = "none",
data = df.pm,
qoi = "att",
lead = c(0:3),
outcome.var = 'nat_rate_ord',
match.missing = TRUE,
placebo.test = TRUE)We can estimate the ATT and dynamic treatment effects using the
command PanelEstimate. To obtain the ATT, we set the option
pooled = TRUE. The standard error is calculated using a
block bootstrapping method.
# ATT
PE.results.pool <- PanelEstimate(PM.results, data = df.pm, pooled = TRUE)
summary(PE.results.pool)
#> Matches created with 3 lags
#>
#> Standard errors computed with 1000 Weighted bootstrap samples
#>
#> Estimate of Average Treatment Effect on the Treated (ATT) by Period:
#> $summary
#> estimate std.error 2.5% 97.5%
#> [1,] 1.177674 0.3632403 0.4308786 1.862727
#>
#> $lag
#> [1] 3
#>
#> $iterations
#> [1] 1000
#>
#> $qoi
#> [1] "att"Event Study Plot
We can also use the command PanelEstimate to estimate the dynamic
treatment effects at post-treatment periods. Additionally, for
pre-treatment periods, PanelMatch offers the command
placebo_test to compute the change in outcomes for each
pre-treatment period in the lag window compared to the last
pre-treatment period.
# Dynamic Treatment Effects
PE.results <- PanelEstimate(PM.results, data = df.pm)
PE.results.placebo <- placebo_test(PM.results.placebo, data = df.pm, plot = F)est_lead <- as.vector(PE.results$estimates)
est_lag <- as.vector(PE.results.placebo$estimates)
sd_lead <- apply(PE.results$bootstrapped.estimates,2,sd)
sd_lag <- apply(PE.results.placebo$bootstrapped.estimates,2,sd)
coef <- c(est_lag, 0, est_lead)
sd <- c(sd_lag, 0, sd_lead)
pm.output <- cbind.data.frame(ATT=coef, se=sd, t=c(-2:4))
p.pm <- esplot(data = pm.output,Period = 't',
Estimate = 'ATT',SE = 'se')
p.pm
Imputation Method
Liu, Wang, and Xu (2022) proposes a method for estimating the fixed effect counterfactual model by utilizing the untreated observations. The model, represented by the equation \(Y_{it} = \alpha_i + \xi_t + \delta^{TWFE}D_{it} + \epsilon\), estimates the parameters \(\alpha_{i}\) and \(\xi_{t}\) using only untreated observations and imputes the counterfactual using the estimated values of \(\hat{\alpha}{i}\) and \(\hat{\xi}{t}\). This method, independently proposed by Borusyak, Jaravel, and Spiess (2021) and Gardner (2021), is also referred to as the “imputation method” or “two-stage DID”, respectively. To implement this estimator, the fe (fixed effects) method in the fect package can be utilized. The estimated average treatment effect (ATT) is 1.506, with a standard error of 0.197.
out.fect <- fect(nat_rate_ord~indirect, data = df, index = c("bfs","year"),
method = 'fe', se = TRUE)
print(out.fect$est.avg)
#> ATT.avg S.E. CI.lower CI.upper p.value
#> [1,] 1.506416 0.2052678 1.104099 1.908734 2.156053e-13The fect function stores the estimated dynamic treatment effects that
are obtained using the counterfactual method in the object
est.att. The uncertainty estimates for these estimates are
obtained through a nonparametric clustered bootstrap procedure.
fect.output <- as.matrix(out.fect$est.att)
print(fect.output)
#> ATT S.E. CI.lower CI.upper p.value count
#> -17 -0.01985591 0.4878539 -0.976032054 0.936320242 9.675347e-01 89
#> -16 0.07269638 0.2096894 -0.338287228 0.483679982 7.288272e-01 157
#> -15 -0.21365489 0.2657927 -0.734599085 0.307289307 4.214893e-01 162
#> -14 -0.13153052 0.1707489 -0.466192133 0.203131099 4.411127e-01 350
#> -13 0.46902947 0.2230022 0.031953159 0.906105787 3.544390e-02 371
#> -12 0.24357664 0.1944239 -0.137487287 0.624640573 2.102741e-01 398
#> -11 -0.24399485 0.1217989 -0.482716314 -0.005273378 4.514941e-02 417
#> -10 -0.02765866 0.1773139 -0.375187452 0.319870132 8.760433e-01 434
#> -9 -0.09347047 0.1368538 -0.361698965 0.174758024 4.946099e-01 451
#> -8 -0.10638473 0.1493041 -0.399015371 0.186245919 4.761322e-01 464
#> -7 0.17476094 0.1930326 -0.203575925 0.553097806 3.652830e-01 468
#> -6 -0.40297604 0.1577457 -0.712151985 -0.093800092 1.063122e-02 503
#> -5 -0.24326436 0.1736601 -0.583631917 0.097103206 1.612717e-01 503
#> -4 -0.23687833 0.1633786 -0.557094604 0.083337935 1.470939e-01 503
#> -3 0.44244127 0.2235247 0.004340950 0.880541588 4.777274e-02 503
#> -2 -0.11927690 0.1840486 -0.480005472 0.241451671 5.169378e-01 508
#> -1 0.17612194 0.1798590 -0.176395132 0.528639013 3.274700e-01 512
#> 0 0.22491529 0.1884205 -0.144382200 0.594212770 2.326002e-01 514
#> 1 0.83813813 0.2705439 0.307881734 1.368394529 1.948485e-03 514
#> 2 1.81274486 0.3180371 1.189403617 2.436086098 1.199546e-08 425
#> 3 1.91941052 0.3998606 1.135698119 2.703122920 1.585080e-06 357
#> 4 1.19567464 0.3516167 0.506518669 1.884830611 6.726118e-04 352
#> 5 0.93520961 0.4804345 -0.006424742 1.876843963 5.158376e-02 164
#> 6 2.52517848 0.6692102 1.213550597 3.836806354 1.610564e-04 143
#> 7 2.84129736 0.7904270 1.292088977 4.390505751 3.248452e-04 116
#> 8 2.32288547 0.7119625 0.927464665 3.718306279 1.103751e-03 97
#> 9 2.19105962 0.8215706 0.580810849 3.801308382 7.655084e-03 80
#> 10 1.71382906 0.8486171 0.050570108 3.377088003 4.342955e-02 63
#> 11 1.07651548 1.0094509 -0.901971919 3.055002878 2.862263e-01 50
#> 12 -0.30828629 0.6336991 -1.550313785 0.933741197 6.266220e-01 46
#> 13 0.35079211 2.0289255 -3.625828722 4.327412944 8.627336e-01 11
#> 14 0.81696815 2.7207014 -4.515508670 6.149444969 7.639647e-01 11
#> 15 0.97427612 4.1629709 -7.184996838 9.133549086 8.149587e-01 11
#> 16 -2.49173503 1.7783901 -5.977315558 0.993845488 1.611787e-01 11
#> 17 1.42334105 5.0108692 -8.397782168 11.244464274 7.763715e-01 6
#> 18 -0.01384226 4.7741056 -9.370917213 9.343232690 9.976866e-01 2Event Study Plot
We can visualize the dynamic treatment effects in the same way.
fect.output <- as.data.frame(fect.output)
fect.output$Time <- c(-17:18)
p.fect <- esplot(fect.output,Period = 'Time',Estimate = 'ATT',
SE = 'S.E.',CI.lower = "CI.lower",
CI.upper = 'CI.upper',xlim = c(-12,10))
p.fect
The plot function included in fect, by
default, also visualizes the estimated dynamic treatment effects. It
incorporates a bar plot at the bottom of the plot to illustrate the
number of treated units for each period. In this instance, we set the
option proportion = 0.1 to only display those periods whose
number of treated observations exceeds ten percent of the maximum number
of treated observations within a period.
Compare with PanelMatch
To compare the results with those obtained via the
PanelMatch method, fect also provides the
capability to compute the ATT for the subset of the sample that includes
units with 3 pre-treatment periods and 4 post-treatment periods. This
can be achieved by specifying the option balance.period = c(-2,4), which
denotes the range of periods to be included in the analysis (-2, -1, and
0 for pre-treatment, and 1, 2, 3, and 4 for post-treatment).
out.fect.balance <- fect(nat_rate_ord~indirect, data = df, index = c("bfs","year"),
method = 'fe', se = TRUE, balance.period = c(-2,4))
#ATT
print(out.fect.balance$est.balance.avg)
#> ATT.avg S.E. CI.lower CI.upper p.value
#> [1,] 1.711554 0.2088427 1.30223 2.120878 2.220446e-16
#Dynamic Treatment Effectsfect.balance.output <- as.data.frame(out.fect.balance$est.balance.att)
fect.balance.output$Time <- c(-2:4)
p.fect.balance <- esplot(fect.balance.output,Period = 'Time',Estimate = 'ATT',
SE = 'S.E.',CI.lower = "CI.lower",
CI.upper = 'CI.upper')
p.fect.balance
Example 2: w/ Treatment Reversal
We utilize the data from Grumbach and Sahn (2020) to demonstrate the analysis procedure when the treatment can switch on and off. The study, which uses district-election year level TSCS data from US House general elections between 1980 and 2012, argues that the presence of Asian (Black/Latino) candidates in general elections leads to an increase in the proportion of campaign contributions by Asian (Black/Latino) donors. Specifically, we focus on the effects of Asian candidates as presented in the top left panel of Figure 5 in the paper.
data <- read.dta13(paste0(path, "grumbach-sahn_apsr2020_full.dta"))
data$cycle <- as.integer(as.numeric(data$cycle/2))
head(data)
#> cycle district_final general_sharetotal_A_all general_sharetotal_B_all
#> 1 991 AK1 0.000000000 0.00000000
#> 2 992 AK1 0.000000000 0.05061753
#> 3 993 AK1 0.036166365 0.07233273
#> 4 994 AK1 0.015651902 0.01043460
#> 5 995 AK1 0.000000000 0.01882505
#> 6 996 AK1 0.006508764 0.03319470
#> general_sharetotal_H_all cand_A_all cand_B_all cand_H_all dfe
#> 1 0.000000000 0 0 0 1
#> 2 0.000000000 0 0 0 1
#> 3 0.057866184 0 0 0 1
#> 4 0.005217301 0 0 0 1
#> 5 0.030926865 0 0 0 1
#> 6 0.019689011 0 0 0 1Data Visualization
In our analysis, the variable district_final and
cycle serve as the unit and time index, respectively. The
treatment variable is cand_A_all, the outcome of interest
is general_sharetotal_A_all, and the covariates included in
the analysis are cand_H_all and cand_B_all. To
gain insight into the treatment conditions and missing values in the
data, we utilize the panelView package. The option
by.timing = TRUE is set in order to organize the units
based on the timing of receiving the treatment. It should be noted that
this dataset includes missingness in the outcome variable and the
treatment may switch on and off.
y <- "general_sharetotal_A_all"
d <- "cand_A_all"
unit <- "district_final"
time <- "cycle"
controls <- c("cand_H_all", "cand_B_all")
index <- c("district_final", "cycle")
panelview(Y=y, D=d, X=controls, index = index, data = data, xlab = "Time Period", ylab = "Unit",
gridOff = TRUE, by.timing = TRUE, cex.legend=5, cex.axis= 5, cex.main = 10, cex.lab = 5)
TWFE
First, we obtain the ATT using the TWFE estimator.
model.twfe <- feols(general_sharetotal_A_all ~ cand_A_all + cand_H_all + cand_B_all|district_final + cycle,
data=data, cluster = "district_final") #use the clustered standard error
print(model.twfe)
#> OLS estimation, Dep. Var.: general_sharetotal_A_all
#> Observations: 6,847
#> Fixed-effects: district_final: 489, cycle: 17
#> Standard-errors: Clustered (district_final)
#> Estimate Std. Error t value Pr(>|t|)
#> cand_A_all 0.128098 0.025007 5.12256 4.3480e-07 ***
#> cand_H_all 0.012181 0.005274 2.30980 2.1316e-02 *
#> cand_B_all 0.010456 0.004594 2.27615 2.3269e-02 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> RMSE: 0.049535 Adj. R2: 0.275898
#> Within R2: 0.083964Event Study Plot
To estimate the dynamic treatment effects, we first utilize the
function get.cohort to obtain the relative period to the
treatment. We then generate the treated unit indicator
treat by following these steps: First, we set
treat = 0 for observations of units that have never been
treated. Second, for units that have been treated, if the treatment
switches on and off, we only set treat = 1 for observations
prior to the last treatment exit, i.e., if a unit has the treatment
status sequence as 0, 0, 0, 1, 1, 0, 0, we only set
treat = 1 for the first five observations. This procedure
assumes no carryover treatment effects, meaning that the treatment only
takes effect at the current period. Finally, we interact the dummy
variable treat with the relative period index in a two-way
fixed effects regression to estimate the dynamic treatment effects.
data_cohort <- get.cohort(data, index = index, D=d)
# Generate a dummy variable treat
data_cohort$treat <- 0
data_cohort[which(data_cohort$Cohort!='Control'),'treat'] <- 1
data_cohort[which(is.na(data_cohort$Time_to_Treatment)), "treat"] <- 0
# remove observations that starts with treated status
remove <- intersect(which(is.na(data_cohort$Time_to_Treatment)), which(data_cohort[,d]==1))
if(length(remove)>0){data_cohort <- data_cohort[-remove,]}
# replace missingness in Time_to_Treatment with an arbitrary number
data_cohort[which(is.na(data_cohort$Time_to_Treatment)), "Time_to_Treatment"] <- 999
twfe.est <- feols(general_sharetotal_A_all ~ i(Time_to_Treatment, treat, ref = -1) +
cand_H_all +cand_B_all | district_final + cycle,
data = data_cohort, cluster = "district_final")We then visualize the estimated dynamic treatment.
twfe.output <- as.matrix(twfe.est$coeftable[c(1:25),])
print(twfe.output)
#> Estimate Std. Error t value Pr(>|t|)
#> Time_to_Treatment::-16:treat -0.014940426 0.023438304 -0.6374363 5.241396e-01
#> Time_to_Treatment::-15:treat -0.038164998 0.011362006 -3.3590017 8.434783e-04
#> Time_to_Treatment::-14:treat -0.035172442 0.010524624 -3.3419189 8.958865e-04
#> Time_to_Treatment::-13:treat -0.032246260 0.010475484 -3.0782598 2.199328e-03
#> Time_to_Treatment::-12:treat 0.003439330 0.028648925 0.1200509 9.044922e-01
#> Time_to_Treatment::-11:treat -0.021413340 0.017987605 -1.1904497 2.344485e-01
#> Time_to_Treatment::-10:treat -0.019647640 0.012642996 -1.5540335 1.208247e-01
#> Time_to_Treatment::-9:treat -0.001535462 0.013341633 -0.1150880 9.084227e-01
#> Time_to_Treatment::-8:treat -0.002309351 0.012435639 -0.1857043 8.527538e-01
#> Time_to_Treatment::-7:treat -0.003408378 0.012370766 -0.2755188 7.830343e-01
#> Time_to_Treatment::-6:treat 0.016583110 0.018584831 0.8922927 3.726758e-01
#> Time_to_Treatment::-5:treat -0.014949213 0.012092481 -1.2362404 2.169639e-01
#> Time_to_Treatment::-4:treat -0.019672539 0.007123897 -2.7614858 5.970795e-03
#> Time_to_Treatment::-3:treat -0.019037443 0.008972166 -2.1218337 3.435495e-02
#> Time_to_Treatment::-2:treat -0.007612287 0.006002114 -1.2682677 2.053071e-01
#> Time_to_Treatment::0:treat 0.116367293 0.019291177 6.0321511 3.199576e-09
#> Time_to_Treatment::1:treat 0.133724736 0.045502782 2.9388255 3.450662e-03
#> Time_to_Treatment::2:treat 0.316864160 0.125811780 2.5185572 1.210235e-02
#> Time_to_Treatment::3:treat 0.198552123 0.043150455 4.6013912 5.357012e-06
#> Time_to_Treatment::4:treat 0.261498017 0.093699031 2.7908295 5.463286e-03
#> Time_to_Treatment::5:treat 0.155250352 0.027192716 5.7092625 1.972843e-08
#> Time_to_Treatment::6:treat 0.184621686 0.053799796 3.4316429 6.509094e-04
#> Time_to_Treatment::7:treat 0.192064488 0.024909842 7.7103857 7.095591e-14
#> Time_to_Treatment::8:treat 0.066716303 0.024755944 2.6949610 7.282539e-03
#> Time_to_Treatment::9:treat 0.042376250 0.024941526 1.6990239 8.995201e-02twfe.output <- as.data.frame(twfe.output)
twfe.output$Time <- c(c(-16:-2),c(0:9))+1
p.twfe <- esplot(twfe.output,Period = 'Time',Estimate = 'Estimate',
SE = 'Std. Error', xlim = c(-15,1))
p.twfe
PanelMatch
PanelMatch is also equipped to handle scenarios where the treatment can switch on and off. In this case, we will utilize a subset of the treated units that have two pre-treatment periods and one post-treatment period to calculate the average treatment effects.
df.pm <- data_cohort
# we need to convert the unit and time indicator to integer
df.pm[,"district_final"] <- as.integer(as.factor(df.pm[,"district_final"]))
df.pm[,"cycle"] <- as.integer(as.factor(df.pm[,"cycle"]))
df.pm <- df.pm[,c("district_final","cycle","cand_A_all","general_sharetotal_A_all")]
PM.results <- PanelMatch(lag=2,
time.id="cycle",
unit.id = "district_final",
treatment = "cand_A_all",
refinement.method = "none",
data = df.pm,
qoi = "att",
lead = 0,
outcome.var = "general_sharetotal_A_all",
match.missing = TRUE)
## For pre-treatment dynamic effects
PM.results.placebo <- PanelMatch(lag=2,
time.id="cycle",
unit.id = "district_final",
treatment = "cand_A_all",
refinement.method = "none",
data = df.pm,
qoi = "att",
lead = 0,
outcome.var = "general_sharetotal_A_all",
match.missing = TRUE,
placebo.test = TRUE)# ATT
PE.results.pool <- PanelEstimate(PM.results, data = df.pm, pooled = TRUE)
summary(PE.results.pool)$summary
#> Matches created with 2 lags
#>
#> Standard errors computed with 1000 Weighted bootstrap samples
#>
#> Estimate of Average Treatment Effect on the Treated (ATT) by Period:
#> estimate std.error 2.5% 97.5%
#> [1,] 0.1197584 0.0219985 0.07881666 0.1638193Event Study Plot
We obtain the dynamic treatment effects in the same way as we showed in the case of Hainmueller and Hangartner (2019).
# Dynamic Treatment Effects
PE.results <- PanelEstimate(PM.results, data = df.pm)
PE.results.placebo <- placebo_test(PM.results.placebo, data = df.pm, plot = F)est_lead <- as.vector(PE.results$estimates)
est_lag <- as.vector(PE.results.placebo$estimates)
sd_lead <- apply(PE.results$bootstrapped.estimates,2,sd)
sd_lag <- apply(PE.results.placebo$bootstrapped.estimates,2,sd)
coef <- c(est_lag, 0, est_lead)
sd <- c(sd_lag, 0, sd_lead)
pm.output <- cbind.data.frame(ATT=coef, se=sd, t=c(-1:1))
p.pm <- esplot(data = pm.output,Period = 't',
Estimate = 'ATT',SE = 'se')
p.pm
Imputation Method
The stacked DiD, iw estimator and csdid estimator do not account for treatment reversals. To estimate the average treatment effect (ATT) in scenarios where treatment can switch on and off, we employ the imputation method proposed by Liu, Wang, and Xu (2022), which accounts for treatment reversals and heterogeneous treatment effects.
model.fect <- fect(Y = y, D = d, X= controls, data = data, method = "fe",
index = index, se = TRUE, parallel = TRUE, seed = 1234, force = "two-way")
print(model.fect$est.avg)
#> ATT.avg S.E. CI.lower CI.upper p.value
#> [1,] 0.1277491 0.02469477 0.07934829 0.17615 2.302096e-07Event Study Plot
We visualize the estimated dynamic treatment for the counterfactual estimator in the same way.
fect.output <- as.matrix(model.fect$est.att)
print(fect.output)
#> ATT S.E. CI.lower CI.upper p.value count
#> -15 0.0007768005 0.016832653 -0.0322145930 0.033768194 9.631919e-01 11
#> -14 -0.0211704682 0.006041969 -0.0330125107 -0.009328426 4.584942e-04 15
#> -13 -0.0197937240 0.006866693 -0.0332521951 -0.006335253 3.944454e-03 21
#> -12 -0.0161604498 0.006170367 -0.0282541474 -0.004066752 8.817714e-03 24
#> -11 0.0173223287 0.025022821 -0.0317214996 0.066366157 4.887733e-01 26
#> -10 -0.0101671188 0.013181064 -0.0360015294 0.015667292 4.405038e-01 32
#> -9 -0.0086394046 0.009699092 -0.0276492747 0.010370465 3.730667e-01 36
#> -8 0.0091411874 0.011105172 -0.0126245488 0.030906924 4.104245e-01 39
#> -7 0.0077375741 0.009057584 -0.0100149645 0.025490113 3.929584e-01 43
#> -6 0.0076803525 0.009904005 -0.0117311416 0.027091847 4.380564e-01 48
#> -5 0.0260877842 0.016026689 -0.0053239483 0.057499517 1.035734e-01 53
#> -4 -0.0052389859 0.010522568 -0.0258628410 0.015384869 6.185680e-01 56
#> -3 -0.0091852656 0.005278075 -0.0195301034 0.001159572 8.181198e-02 57
#> -2 -0.0099003105 0.007631375 -0.0248575308 0.005056910 1.945222e-01 62
#> -1 0.0006410009 0.005232789 -0.0096150766 0.010897078 9.025054e-01 65
#> 0 0.0166019901 0.008119382 0.0006882944 0.032515686 4.088092e-02 73
#> 1 0.1233181317 0.019883053 0.0843480634 0.162288200 5.568890e-10 74
#> 2 0.1389600839 0.044001797 0.0527181458 0.225202022 1.588257e-03 13
#> 3 0.3300925450 0.145937099 0.0440610865 0.616124004 2.370467e-02 5
#> 4 0.2025269124 0.070553138 0.0642453035 0.340808521 4.097473e-03 3
#> 5 0.2659320021 0.126480846 0.0180340988 0.513829905 3.550533e-02 3
#> 6 0.1591613175 0.025561200 0.1090622864 0.209260349 4.764340e-10 3
#> 7 0.1738020810 0.093887870 -0.0102147628 0.357818925 6.414560e-02 2
#> 8 0.1082956579 0.003457040 0.1015199848 0.115071331 0.000000e+00 1
#> 9 -0.0188593819 0.002655091 -0.0240632652 -0.013655499 1.219913e-12 1
#> 10 -0.0438930771 0.002543391 -0.0488780322 -0.038908122 0.000000e+00 1fect.output <- as.data.frame(fect.output)
fect.output$Time <- c(-15:10)
p.fect <- esplot(fect.output,Period = 'Time',Estimate = 'ATT',
SE = 'S.E.',CI.lower = "CI.lower",
CI.upper = 'CI.upper', xlim = c(-15,1))
p.fect
The plot function shipped with fect can also visualize the estimated dynamic treatment effects.
plot(model.fect, stats = "F.p")
An \(F\)-test for the null hypothesis of zero residual averages in the pre-treatment periods yields a p-value greater than 0.1, providing insufficient evidence to reject the null hypothesis of zero pre-treatment average residuals.
We can visualize the period-wise ATTs relative to the
exit of treatments by setting
type = "exit".
plot(model.fect, stats = "F.p", type = 'exit')
Diagnostic Tests
Fect also offers a placebo test function, which
utilizes a subset of the data by removing observations in a specified
range for model fitting, and subsequently assessing whether the
estimated average treatment effect (ATT) in this range is statistically
significant from zero. In this instance, we have set
placebo.period = c(-2, 0), indicating that the placebo
periods will consist of the three periods preceding the onset of
treatment.
out.fect.p <- fect(Y = y, X = controls, D = d, data = data, index = index,
method = 'fe', se = TRUE, placeboTest = TRUE, placebo.period = c(-2,0))
p.placebo <- plot(out.fect.p, proportion = 0.1, stats = "placebo.p")
p.placebo
Test for (No) Carryover Effects
The concept of the placebo test can be extended to examine the
existence of carryover effects. Rather than obscuring a few periods
immediately preceding the onset of treatment, this method involves
obscuring a few periods immediately following the termination of
treatment. If carryover effects do not exist, one would expect the
average prediction error in these periods to be close to zero. To
perform the carryover test, the option carryoverTest = TRUE
is set. Additionally, a range of post-treatment periods can be specified
in the option carryover.period to remove observations from
this range for the purposes of model fitting, and then test whether the
estimated Average Treatment Effect (ATT) in this range is significantly
different from zero.
out.fect.c <- fect(Y = y, X = controls, D = d, data = data, index = index,
method = 'fe', se = TRUE, carryoverTest = TRUE, carryover.period = c(1,2))
p.carryover <- plot(out.fect.c, stats = "carryover.p")
p.carryover
Compare with PanelMatch
We can compare it with the estimated ATT and dynamic treatment effects on the “balanced” sample using fect.
out.fect.balance <- fect(Y = y, X = controls, D = d, data = data, index = index,
method = 'fe', se = TRUE, balance.period = c(-1,1))
#ATT
print(out.fect.balance$est.balance.avg)
#> ATT.avg S.E. CI.lower CI.upper p.value
#> [1,] 0.1330911 0.02246486 0.08906074 0.1771214 3.134189e-09
#Dynamic Treatment Effectsfect.balance.output <- as.data.frame(out.fect.balance$est.balance.att)
fect.balance.output$Time <- c(-1:1)
p.fect.balance <- esplot(fect.balance.output,Period = 'Time',Estimate = 'ATT',
SE = 'S.E.',CI.lower = "CI.lower",
CI.upper = 'CI.upper')
p.fect.balance