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", "didimputation",
"panelView", "PanelMatch", "ggplot2", "bacondecomp", "HonestDiD")
install.packages(setdiff(packages, rownames(installed.packages())),
repos = "https://cloud.r-project.org")
# install "HonestDID"
if ("HonestDID" %in% rownames(installed.packages()) == FALSE) {
devtools:: install_github("asheshrambachan/HonestDiD")
}
# update "fect"
devtools:: install_github("xuyiqing/fect", upgrade = "always")
# 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)
library(didimputation)
library(doParallel)
library(HonestDiD)Lood data:
data(paneltools)
ls()
#> [1] "gs2020" "hh2019" "packages"
summary(hh2019)
#> bfs year nat_rate_ord indirect
#> Min. : 1 Min. :1991 Min. : 0.0000 Min. :0.0000
#> 1st Qu.: 957 1st Qu.:1995 1st Qu.: 0.0000 1st Qu.:0.0000
#> Median :3219 Median :2000 Median : 0.3849 Median :0.0000
#> Mean :3182 Mean :2000 Mean : 2.4710 Mean :0.3411
#> 3rd Qu.:5229 3rd Qu.:2005 3rd Qu.: 3.3058 3rd Qu.:1.0000
#> Max. :6800 Max. :2009 Max. :46.6667 Max. :1.0000
summary(gs2020)
#> cycle district_final general_sharetotal_A_all
#> Min. :1980 Length:7007 Min. :0.00000
#> 1st Qu.:1988 Class :character 1st Qu.:0.00000
#> Median :1996 Mode :character Median :0.00473
#> Mean :1996 Mean :0.02305
#> 3rd Qu.:2004 3rd Qu.:0.02005
#> Max. :2012 Max. :1.00000
#> NA's :160
#> general_sharetotal_B_all general_sharetotal_H_all cand_A_all
#> Min. :0.00000 Min. :0.00000 Min. :0.00000
#> 1st Qu.:0.00410 1st Qu.:0.00000 1st Qu.:0.00000
#> Median :0.02819 Median :0.00367 Median :0.00000
#> Mean :0.05400 Mean :0.02536 Mean :0.01713
#> 3rd Qu.:0.06569 3rd Qu.:0.01694 3rd Qu.:0.00000
#> Max. :1.00000 Max. :1.00000 Max. :1.00000
#> NA's :160 NA's :160
#> cand_B_all cand_H_all dfe
#> Min. :0.000 Min. :0.0000 Min. : 1.0
#> 1st Qu.:0.000 1st Qu.:0.0000 1st Qu.:122.0
#> Median :0.000 Median :0.0000 Median :241.0
#> Mean :0.102 Mean :0.0598 Mean :241.3
#> 3rd Qu.:0.000 3rd Qu.:0.0000 3rd Qu.:363.0
#> Max. :1.000 Max. :1.0000 Max. :495.0
#> 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 <- hh2019
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 histories: 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. Notably,
in order to maintain consistency with canonical TWFE regressions, we
specify the start0 option as TRUE, thereby
designating period \(0\) as the first
post-treatment period, while period \(-1\) represents the last pre-treatment
period.
df.use <- get.cohort(df.use,D="indirect",index=c("bfs","year"),start0 = TRUE)
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 noteworthy that
esplot() considers period \(0\) to be the last pre-treatment period by
default, so it is necessary to shift the period index by one.
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
In addition, by setting the start0 option to
TRUE, we have the flexibility to designate period \(0\) as the first post-treatment period,
thereby rendering period \(-1\) as the
last pre-treatment period. Consequently, there is no requirement to
shift the period index by one in this scenario, as it aligns with the
conventional TWFE framework.
twfe.output <- as.data.frame(twfe.est$coeftable)
twfe.output$Time <- c(c(-18:-2),c(0:17))
p.twfe <- esplot(twfe.output,Period = 'Time',Estimate = 'Estimate',
SE = 'Std. Error', xlim = c(-12,10),start0 = TRUE)
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
The package also offers an option, ‘base_period = universal’, which sets the base period to always be the last pre-treatment period. This produces results that are quite similar to the iw estimator.
cs.est.1.u <- 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",
base_period = "universal")
cs.att.1.u <- aggte(cs.est.1.u, type = "dynamic",
bstrap=FALSE, cband=FALSE, na.rm=T)
cs.output.u <- cbind.data.frame(Estimate = cs.att.1.u$att.egt,
SE = cs.att.1.u$se.egt,
time = cs.att.1.u$egt + 1)
p.cs.1.u <- esplot(cs.output.u,Period = 'time',Estimate = 'Estimate',
SE = 'SE', xlim = c(-12,10))
p.cs.1.u
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
And we can set the reference period to be the last pre-treatment period.
cs.est.2.u <- 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",
base_period = "universal")
cs.att.2.u <- aggte(cs.est.1.u, type = "dynamic",
bstrap=FALSE, cband=FALSE, na.rm=T)
cs.output.u <- cbind.data.frame(Estimate = cs.att.2.u$att.egt,
SE = cs.att.2.u$se.egt,
time = cs.att.2.u$egt + 1)
p.cs.2.u <- esplot(cs.output.u,Period = 'time',Estimate = 'Estimate',
SE = 'SE', xlim = c(-12,10))
p.cs.2.u
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.3553834 0.4901864 1.852845
#>
#> $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.1996611 1.115088 1.897745 4.52971e-14The 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.4920597 -0.984275169 0.94456336 9.678120e-01 89
#> -16 0.07269638 0.2459788 -0.409413117 0.55480587 7.675820e-01 157
#> -15 -0.21365489 0.2725465 -0.747836219 0.32052644 4.330865e-01 162
#> -14 -0.13153052 0.1502720 -0.426058234 0.16299720 3.814200e-01 350
#> -13 0.46902947 0.2063744 0.064543137 0.87351581 2.304355e-02 371
#> -12 0.24357664 0.1911275 -0.131026371 0.61817966 2.025147e-01 398
#> -11 -0.24399485 0.1330149 -0.504699353 0.01670966 6.660331e-02 417
#> -10 -0.02765866 0.1789370 -0.378368771 0.32305145 8.771587e-01 434
#> -9 -0.09347047 0.1593950 -0.405879027 0.21893809 5.576016e-01 451
#> -8 -0.10638473 0.1791419 -0.457496427 0.24472697 5.526076e-01 464
#> -7 0.17476094 0.2100882 -0.237004365 0.58652625 4.054961e-01 468
#> -6 -0.40297604 0.1564322 -0.709577507 -0.09637457 9.993826e-03 503
#> -5 -0.24326436 0.1761838 -0.588578188 0.10204948 1.673582e-01 503
#> -4 -0.23687833 0.1716347 -0.573276220 0.09951955 1.675464e-01 503
#> -3 0.44244127 0.2248854 0.001673899 0.88320864 4.913627e-02 503
#> -2 -0.11927690 0.2150376 -0.540742921 0.30218912 5.791141e-01 508
#> -1 0.17612194 0.1984314 -0.212796424 0.56504030 3.747716e-01 512
#> 0 0.22491529 0.2077587 -0.182284187 0.63211476 2.789950e-01 514
#> 1 0.83813813 0.3030896 0.244093373 1.43218289 5.686795e-03 514
#> 2 1.81274486 0.3270295 1.171778756 2.45371096 2.972281e-08 425
#> 3 1.91941052 0.3874496 1.160023249 2.67879779 7.271751e-07 357
#> 4 1.19567464 0.3374208 0.534342062 1.85700722 3.947464e-04 352
#> 5 0.93520961 0.4367611 0.079173540 1.79124568 3.225483e-02 164
#> 6 2.52517848 0.6401594 1.270489047 3.77986790 7.993038e-05 143
#> 7 2.84129736 0.7979244 1.277394204 4.40520052 3.696419e-04 116
#> 8 2.32288547 0.6932092 0.964220454 3.68155049 8.054483e-04 97
#> 9 2.19105962 0.7773284 0.667524023 3.71459521 4.821776e-03 80
#> 10 1.71382906 0.8242300 0.098367971 3.32929014 3.758893e-02 63
#> 11 1.07651548 1.0376867 -0.957313148 3.11034411 2.995408e-01 50
#> 12 -0.30828629 0.5806224 -1.446285220 0.82971263 5.954476e-01 46
#> 13 0.35079211 2.2190429 -3.998452142 4.70003636 8.743917e-01 11
#> 14 0.81696815 2.6023888 -4.283620131 5.91755643 7.535741e-01 11
#> 15 0.97427612 4.2152885 -7.287537524 9.23608977 8.172145e-01 11
#> 16 -2.49173503 1.9961115 -6.404041662 1.42057159 2.119232e-01 11
#> 17 1.42334105 5.4646049 -9.287087736 12.13376984 7.945047e-01 6
#> 18 -0.01384226 4.8005838 -9.422813627 9.39512910 9.976993e-01 2Borusyak, Jaravel, and Spiess (2021)
also provides a package didimputation to calculate the
average treatment effect in the same way as Liu,
Wang, and Xu (2022).
df.impute <- df.use
df.impute[which(is.na(df.impute$FirstTreat)),"FirstTreat"] <- 0 # replace NA with 0
out.impute <- did_imputation(data = df.impute,
yname = "nat_rate_ord",
gname = "FirstTreat",
tname = "year",
idname = "bfs",
cluster_var = "bfs")
out.impute
#> # A tibble: 1 × 6
#> lhs term estimate std.error conf.low conf.high
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 nat_rate_ord treat 1.51 0.187 1.14 1.87Event 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 didimputation package also supports the event
study estimate, which produces a numerically equivalent estimate for the
dynamic treatment effects post-treatment. However,
didimputation uses a distinct TWFE-regression solely on
untreated observations to obtain the pre-treatment estimates and conduct
a test for parallel trends. In this case, we can specify the option
pretrends = c(-13:-1) to compute the regression
coefficients for the pre-treatment periods ranging from -13 to -1.
model.impute <- did_imputation(data = df.impute,
yname = "nat_rate_ord",
gname = "FirstTreat",
tname = "year",
idname = "bfs",
cluster_var = "bfs",
pretrends = c(-13:-1),
horizon = TRUE)
model.impute$term <- as.numeric(model.impute$term)+1 #set 1 as the first post-treatment period
to_plot <- as.data.frame(model.impute)
esplot(data=to_plot,Period = "term",Estimate = 'estimate',SE = 'std.error',xlim = c(-12,10))
out.impute
#> # A tibble: 1 × 6
#> lhs term estimate std.error conf.low conf.high
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 nat_rate_ord treat 1.51 0.187 1.14 1.87Compare 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.2225988 1.275269 2.14784 1.487699e-14
#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
Sensitivity Analyses and Robust Confidence Set
Rambachan and Roth (2023) introduce a
robust inference method that tolerates violations of the parallel trend
assumption (PTA), provided that deviations in post-treatment periods do
not exceed those observed in the pre-treatment phase. This method is
useful for sensitivity analysis of estimates from FEct
and other models. We demonstrate two approaches for this analysis: one
using estimated coefficients from all 15 pre-treatment periods and
another using coefficients from three specific placebo periods. The
commands createSensitivityResults_relativeMagnitudes and
constructOriginalCS from the Honestdid
package are utilized to construct these robust confidence intervals.
For the first approach, we start with extracting event-study
coefficients of all pre-treatment periods and their variance-covariance
matrix from FEct’s output. By setting a threshold \(\bar{M}\), we allow deviations in
post-treatment periods as long as they are less than \(\bar{M}\) times the maximum difference
between consecutive pre-treatment period coefficients. We generate
confidence intervals for the ATT by varying \(M\) from 0.1 to 0.5, thus adjusting the
allowed level of post-treatment deviations. The l_vec
option dictates the relative weights of post-treatment coefficients in
estimating ATT, proportional to the observation count in each period.
Notably, at \(M=0.3\), the confidence
interval’s lower bound remains slightly above zero, but at \(M=0.4\), the confidence interval includes
zero. This suggests we can still reject the null hypothesis of no
treatment effect when deviations in post-treatment periods do not exceed
\(0.3\) times the maximum difference
between consecutive pre-treatment period coefficients.
# obtain the estimated beta and estimated vcov matrix from the output of FEct
index.use <- which(rownames(out.fect$est.att) %in% as.character(c(-14:10)))
beta.hat <- out.fect$est.att[index.use,1]
vcov.hat <- out.fect$att.vcov[index.use,index.use]
count <- out.fect$count[which(rownames(out.fect$est.att) %in% as.character(c(1:10)))]
# construct the robust confidence intervals under M=0.1,...,0.5.
honest.result <- createSensitivityResults_relativeMagnitudes(betahat = beta.hat,
sigma = vcov.hat,
numPrePeriods = 15,
numPostPeriods = 10,
l_vec = count/sum(count),
Mbarvec = seq(0.1,0.5,by=0.1))
# construct the robust confidence intervals under M=0
originalResults <- constructOriginalCS(betahat = beta.hat,
sigma = vcov.hat,
numPrePeriods = 15,
numPostPeriods = 10,
l_vec = count/sum(count))
createSensitivityPlot_relativeMagnitudes(honest.result, originalResults)
By modifying the l_vec option, we can also obtain
separate confidence intervals for each post-treatment period. In the
figure below, we apply the same command with \(\bar{M}=0.5\) to derive robust confidence
intervals. These intervals allow deviations from the parallel trend
assumption in post-treatment periods, provided they remain smaller than
0.5 times the maximum deviation observed between two consecutive
pre-treatment event study coefficients. These robust confidence
intervals are depicted in dark red in the figure.
T.post <- 10
dte_base <- rep(0,T.post)
dte_output <- cbind.data.frame(lb = c(),
ub = c(),
method = c(),
Delta = c(),
Mbar = c(),
post = c())
for(t.post in c(1:T.post)){
dte_l <- dte_base
dte_l[t.post] <- 1
honest.dte <- createSensitivityResults_relativeMagnitudes(betahat = beta.hat,
sigma = vcov.hat,
numPrePeriods = 15,
numPostPeriods = T.post,
l_vec = dte_l,
Mbarvec = c(0,0.5))
honest.dte <- as.data.frame(honest.dte)
honest.dte$post <- t.post
dte_output <- rbind(dte_output,honest.dte)
}dte_output <- as.data.frame(dte_output)
p.fect.honest <- esplot(fect.output,Period = 'Time',Estimate = 'ATT',
SE = 'S.E.',CI.lower = "CI.lower", main = "FEct",ylab = "Coefficients and 95% CI",
CI.upper = 'CI.upper',xlim = c(-12,10),ylim = c(-6,8),Count = 'count',show.count = T)
p.fect.honest <- p.fect.honest + geom_linerange(aes(x=post,ymin=lb,ymax=ub),data = dte_output[which(dte_output$Mbar==0.5),],color = 'maroon',size=1.2)
p.fect.honest <- p.fect.honest + geom_linerange(aes(x=post,ymin=lb,ymax=ub),data = dte_output[which(dte_output$Mbar==0),],color = 'black',size=1.2)
p.fect.honest
We can also apply the above procedure to event-study results derived
from other estimators, provided their pre-treatment coefficients and
variance-covariance matrices are available. Another concern, highlighted
by Roth (2024), is that some HTE-robust
estimators produce pre-treatment and post-treatment coefficients
asymmetrically, necessitating cautious interpretation by researchers.
For instance, in the CSDID section, omitting the
base_period = "universal" option exposes us to this
risk.
This concern about asymmetry motivates a second approach to obtaining
robust confidence intervals. We remove observations in specified placebo
periods within pre-treatment for model fitting, impute the
counterfactual for these and the post-treatment periods, and calculate
the event-study coefficients. This method ensures symmetric generation
of coefficients for both post-treatment and placebo periods. By setting
a threshold \(\bar{M}\), we allow
deviations in post-treatment periods as long as they are less than \(\bar{M}\) times the maximum difference
between consecutive placebo period coefficients. FEct includes built-in
functionality for this placebo test. We can set
placeboTest = TRUE and
placebo.period = c(-2, 0), which defines the placebo
periods as the three periods immediately preceding treatment onset.
After obtaining the coefficients and variance-covariance matrix for the
placebo and post-treatment periods, we can generate the robust
confidence interval in the same manner as previously described. The
result suggests we can still reject the null hypothesis of no treatment
effect when deviations in post-treatment periods do not exceed \(0.6\) times the maximum difference between
consecutive placebo period coefficients.
out.fect.placebo <- fect(nat_rate_ord~indirect, data = df, index = c("bfs","year"),
method = 'fe', se = TRUE,
placeboTest = TRUE, placebo.period = c(-2,0))index.use <- which(rownames(out.fect.placebo$est.att) %in% as.character(c(-2:T.post)))
beta.hat.p<- out.fect.placebo$est.att[index.use,1]
vcov.hat.p <- out.fect.placebo$att.vcov[index.use,index.use]
count <- out.fect.placebo$count[which(rownames(out.fect.placebo$est.att) %in% as.character(c(1:T.post)))]
honest.result.p <- createSensitivityResults_relativeMagnitudes(betahat = beta.hat.p,
sigma = vcov.hat.p,
numPrePeriods = 3,
numPostPeriods = T.post,
l_vec = count/sum(count),
Mbarvec = seq(0.2,1,by=0.2))
#> Warning in .ARP_computeCI(betahat = betahat, sigma = sigma, numPrePeriods =
#> numPrePeriods, : CI is open at one of the endpoints; CI length may not be
#> accurate
#> Warning in .ARP_computeCI(betahat = betahat, sigma = sigma, numPrePeriods =
#> numPrePeriods, : CI is open at one of the endpoints; CI length may not be
#> accurate
#> Warning in .ARP_computeCI(betahat = betahat, sigma = sigma, numPrePeriods =
#> numPrePeriods, : CI is open at one of the endpoints; CI length may not be
#> accurate
#> Warning in .ARP_computeCI(betahat = betahat, sigma = sigma, numPrePeriods =
#> numPrePeriods, : CI is open at one of the endpoints; CI length may not be
#> accurate
originalResults.p <- constructOriginalCS(betahat = beta.hat.p,
sigma = vcov.hat.p,
numPrePeriods = 3,
numPostPeriods = T.post,
l_vec = count/sum(count))
createSensitivityPlot_relativeMagnitudes(honest.result.p, originalResults.p)
By modifying the l_vec option, we can also obtain
separate confidence intervals for each post-treatment period similarly.
In the figure below, we use the same command with \(\bar{M}=0.5\) to derive robust confidence
intervals using coefficients from the placebo periods. We depict the
placebo periods in blue and the robust confidence intervals in pink.
dte_base <- rep(0,T.post)
dte_output <- cbind.data.frame(lb = c(),
ub = c(),
method = c(),
Delta = c(),
Mbar = c(),
post = c())
for(t.post in c(1:T.post)){
dte_l <- dte_base
dte_l[t.post] <- 1
honest.dte <- createSensitivityResults_relativeMagnitudes(betahat = beta.hat.p,
sigma = vcov.hat.p,
numPrePeriods = 3,
numPostPeriods = T.post,
l_vec = dte_l,
Mbarvec = c(0,0.5))
honest.dte <- as.data.frame(honest.dte)
honest.dte$post <- t.post
dte_output <- rbind(dte_output,honest.dte)
}fect.output.p <- as.data.frame(out.fect.placebo$est.att)
fect.output.p$Time <- as.numeric(rownames(fect.output.p))
p.placebo.honest <- esplot(fect.output,Period = 'Time',Estimate = 'ATT',
SE = 'S.E.',CI.lower = "CI.lower", main = "FEct",ylab = "Coefficients and 95% CI",
CI.upper = 'CI.upper',xlim = c(-12,10),ylim = c(-6,8),Count = 'count',show.count = T)
dte_output <- as.data.frame(dte_output)
p.placebo.honest <- p.placebo.honest + geom_linerange(aes(x=post,ymin=lb,ymax=ub),data = dte_output[which(dte_output$Mbar==0.5),],color = 'maroon1',linewidth = 1)
p.placebo.honest <- p.placebo.honest + geom_linerange(aes(x=post,ymin=lb,ymax=ub),data = dte_output[which(dte_output$Mbar==0),],color = 'black',linewidth = 1)
p.placebo.honest <- p.placebo.honest + geom_linerange(aes(x=Time,ymin=CI.lower,ymax=CI.upper),data = fect.output.p[which(fect.output.p$Time%in%c(-2:0)),],color = 'blue',linewidth = 1)
p.placebo.honest
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 <- gs2020
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,start0 = TRUE)
# 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.02200166 0.07457539 0.1633995Event 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.02234626 0.08929319 0.1768889 2.587186e-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