Nothing
## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ---- echo=FALSE, results="hide", warning=FALSE, message=FALSE----------------
library(did)
# Source the currently version of the did package (based on our Dropbox)
# fldr <- here::here("R/")
# sapply(paste0(fldr,list.files(fldr)), source)
# Source simulation designs
source(here::here("vignettes/setup_sims.R"))
## ----echo=FALSE---------------------------------------------------------------
time.periods <- 5
reset.sim()
te <- 0
te.e <- c(-1,rep(1,time.periods-1))
bett <- betu <- rep(0,time.periods)
set.seed(12345)
# generate data set with these parameters
dta <- build_sim_dataset()
dta$G <- ifelse(dta$G==0, 0, dta$G+1)
dta$G <- ifelse(dta$G==6, 0, dta$G)
## -----------------------------------------------------------------------------
nrow(dta)
head(dta)
## ---- fig.width=8,fig.height=5, fig.align='center', out.width="90%", dpi = 200----
# estimate group-time average treatment effects using att_gt method
# (and ignoring pre-treatment "dip")
attgt.ignoredip <- att_gt(yname = "Y",
tname = "period",
idname = "id",
gname = "G",
xformla = ~1,
data = dta,
)
# summarize the results
summary(attgt.ignoredip)
# make dynamic effects plot
p <- ggdid(aggte(attgt.ignoredip, "dynamic"))
# add actual treatment effects to the plot
library(ggplot2)
truth <- cbind.data.frame(e = seq(-3,2), att.e = c(0,0,-1,1,1,1))
p <- p + geom_line(data = truth, aes(x = e, y = att.e), inherit.aes = FALSE, color = "blue")
p <- p + geom_point(data = truth, aes(x = e, y = att.e), inherit.aes = FALSE, color = "blue")
p
## -----------------------------------------------------------------------------
compute.attgt <- function(dta) {
# pick up all groups
groups <- unique(dta$G)
# pick up all time periods
time.periods <- unique(dta$period)
# sort the groups and drop the untreated group
groups <- sort(groups)[-1]
# sort the time periods and drop the first two
# (can't compute treatment effects for these two
# periods with one-period anticipation -- w/o anticipation
# we would just drop one period here)
time.periods <- sort(time.periods)[-c(1,2)]
# drop last time period (because we don't know if
# these units are affected by anticipation or not
# and we are being very careful)
# (if you were worried about more than one anticipation
# period here, would need to drop more time periods
# from the end)
time.periods <- time.periods[-length(time.periods)]
# list to store all group-time average treatment effects
# that we calculate
attgt.list <- list()
counter <- 1
# loop over all groups
for (g in groups) {
# get the correct "base" period for this group
# (subtract 2 to avoid anticipation)
main.base.period <- g - 2
# loop over all time periods
for (tp in time.periods) {
#----------------------------------------------------
# if it's a pre-treatment time period (used for the
# pre-test, we need to adjust the base period)
# group not treated yet
if (tp < g) {
# move two periods before
base.period <- tp - 2
} else {
# this is a post-treatment period
base.period <- main.base.period
}
#----------------------------------------------------
#----------------------------------------------------
# now, all we need to do is collect the right subset
# of the data and estimate a 2x2 DiD
# get group g and untreated group
this.data <- subset(dta, G==g | G==0)
# get current period and base period data
this.data <- subset(this.data, period==tp | period==base.period)
# set up to compute 2x2 estimator
Ypost <- subset(this.data, period==tp)$Y
Ypre <- subset(this.data, period==base.period)$Y
# dummy variable being in group g
G <- 1*(subset(this.data, period==tp)$G == g)
# compute 2x2 estimator using DRDID package
# (in this unconditional case, it would be straightforward
# to calculate the 2x2 att just using averages, but we
# like the DRDID package as it will work for more complicated
# cases as well)
attgt <- DRDID::reg_did_panel(Ypost, Ypre, G, covariates=NULL)$ATT
# save results
attgt.list[[counter]] <- list(att=attgt, group=g, time.period=tp)
counter <- counter+1
#----------------------------------------------------
}
}
#-----------------------------------------------------------------------------
# aggregate into dynamic effects
# turn results into a data.frame
attgt.results <- do.call("rbind.data.frame", attgt.list)
# add event time to the results
attgt.results$e <- attgt.results$time.period - attgt.results$group
# calculate relative sizes of each group
# (will be used as weights)
n.group <- sapply(groups, function(gg) nrow(subset(dta, G==gg)))
# merge in group sizes
ngroup.mat <- cbind(groups, n.group)
attgt.results <- merge(attgt.results, ngroup.mat, by.x = "group", by.y = "groups")
# event times to calculate dynamic effects
eseq <- unique(attgt.results$e)
eseq <- sort(eseq)
# calculate average effects by event time
att.e <- c()
counter <- 1
for (this.e in eseq) {
# get subset of results at this event time
res.e <- subset(attgt.results, e==this.e)
# calculate weights by group size
res.e$weight <- res.e$n.group / sum(res.e$n.group)
# calculate dynamic effect as weighted average
att.e[counter] <- sum(res.e$att * res.e$weight)
# on to the next one
counter <- counter+1
}
# store dynamic effects results
dyn.results <- cbind.data.frame(e = eseq, att.e = att.e)
# return group-time average treatment effects and dynamic effects
return(list(attgt.results=attgt.results[,c("group","att","time.period")],
dyn.results=dyn.results))
}
## -----------------------------------------------------------------------------
anticipation.results <- compute.attgt(dta)
anticipation.results
## ----bootstrap, eval=FALSE----------------------------------------------------
# # the number of bootstrap iterations
# biters <- 100
#
# # list to store bootstrap results
# boot.res <- list()
#
# # loop for each nonparametric bootstrap iteration
# for (b in 1:biters) {
# # draw a bootstrap sample; here, we'll call an outside function
# bdata <- BMisc::blockBootSample(dta, "id")
#
# # call our function for estimating dynamic effects on the
# # bootstrapped data
# boot.res[[b]] <- compute.attgt(bdata)$dyn.results$att.e
# }
#
# # use the bootstrapped results to compute standard errors
# boot.res <- t(simplify2array(boot.res))
# boot.se <- apply(boot.res, 2, sd)
#
# # add the standard errors to the main results
# dyn.results <- anticipation.results$dyn.results
# dyn.results$att.se <- boot.se
## ---- eval=FALSE--------------------------------------------------------------
# p <- ggplot(data=dyn.results, aes(x=e, y=att.e)) +
# geom_line() +
# geom_point() +
# geom_errorbar(aes(ymin=(att.e-1.96*att.se), ymax=(att.e+1.96*att.se)), width=0.1) +
# theme_bw()
#
# p <- p + geom_line(data=truth, aes(x=e, y=att.e), inherit.aes=FALSE, color="blue")
# p <- p + geom_point(data=truth, aes(x=e, y=att.e), inherit.aes=FALSE, color="blue")
# p
#
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.