## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 6, fig.height = 4)

## ----setup--------------------------------------------------------------------
library(panglm)

## ----did-data-----------------------------------------------------------------
set.seed(1)
d <- expand.grid(id = 1:200, year = 2015:2022)
d$first <- ifelse(d$id <= 100, 2019, NA)
a <- rnorm(200)
e <- d$year - d$first
post <- !is.na(e) & e >= 0
d$visits <- 3 + a[d$id] + 0.1 * (d$year - 2015) +
  ifelse(post, 1 + 0.3 * e, 0) + rnorm(nrow(d))
d$admitted <- rbinom(nrow(d), 1, plogis(-1 + a[d$id] +
                                          0.05 * (d$year - 2015) + 0.8 * post))

## ----did----------------------------------------------------------------------
fit <- did(visits ~ 1, data = d, index = c("id", "year"), treat = "first")
coef(fit)
confint(fit)

## ----did-binomial-------------------------------------------------------------
fitb <- did(admitted ~ 1, data = d, index = c("id", "year"),
            treat = "first", family = "binomial")
exp(cbind(OR = coef(fitb), confint(fitb)))

## ----eventstudy, fig.alt = "Event-study estimates with 95% confidence intervals: flat before adoption, positive and rising after"----
es <- eventstudy(visits ~ 1, data = d, index = c("id", "year"),
                 treat = "first", leads = 4, lags = 3)
es
plot(es)

## ----stag-data----------------------------------------------------------------
set.seed(1)
s <- expand.grid(province = 1:120, year = 2012:2021)
s$adoption <- rep(c(2016, 2018, 2020, NA), each = 30)[s$province]
s$urban <- rep(runif(120), 10)[s$province]
s$pop <- rep(round(runif(120, 50, 500)), 10)[s$province]
s$region <- ceiling(s$province / 10)
a <- rnorm(120)
e <- s$year - s$adoption
post <- !is.na(e) & e >= 0
s$visits <- 50 + 5 * a[s$province] + 0.5 * (s$year - 2012) * (1 + s$urban) +
  ifelse(post, 2 + e, 0) + rnorm(nrow(s))
# true ATT averaged over the treated province-years
mean(2 + e[post])
coef(did(visits ~ 1, data = s, index = c("province", "year"),
         treat = "adoption"))

## ----attgt--------------------------------------------------------------------
cs <- attgt(visits ~ 1, data = s, index = c("province", "year"),
            treat = "adoption")
head(cs$estimates)
summary(cs)

## ----attgt-agg----------------------------------------------------------------
summary(cs, type = "dynamic")
summary(cs, type = "group")$estimates

## ----attgt-boot, fig.alt = "Dynamic ATT by event time with pointwise intervals and wider uniform bands"----
csb <- attgt(visits ~ urban, data = s, index = c("province", "year"),
             treat = "adoption", control = "notyet", weights = "pop",
             cluster = "region", cband = TRUE, biters = 199, seed = 1)
summary(csb, type = "dynamic", min_e = -3)
plot(csb)

## ----etwfe--------------------------------------------------------------------
et <- etwfe(visits ~ urban, data = s, index = c("province", "year"),
            treat = "adoption")
summary(et, type = "event")
summary(etwfe(visits ~ urban, data = s, index = c("province", "year"),
              treat = "adoption", ivar = TRUE))$overall

## ----etwfe-pois---------------------------------------------------------------
s$admissions <- rpois(nrow(s), exp(1 + 0.3 * a[s$province] +
                                     0.03 * (s$year - 2012) + 0.2 * post))
s$rural <- as.numeric(s$urban < 0.5)
etp <- etwfe(admissions ~ 1, data = s, index = c("province", "year"),
             treat = "adoption", family = "poisson", xvar = "rural")
summary(etp)