inst/doc/extensions.R

In did: Treatment Effects with Multiple Periods and Groups

## ----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
#