## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width = 8,
  fig.height = 8,
  message = FALSE,
  warning = FALSE
)
library(ggplot2); library(dplyr); library(tidyr)
data.table::setDTthreads(1L)
options(dplyr.summarise.inform = FALSE, scipen = 999, digits = 5)
theme_set(theme_bw(base_size = 12))

## ----draw_data----------------------------------------------------------------
library(childpen)

data <- simulate_data(n_individuals = 2000, treatment_groups = 24:28)
head(data)

## ----estimation---------------------------------------------------------------
res <- multiple_treatment_group_analysis(
  data            = data,
  treatment_groups = 24:25,
  periods_post    = 2,
  periods_pre     = NULL,
  verbose         = FALSE
)

## ----did_plot-----------------------------------------------------------------
res |>
  filter(method %in% c("DID_Female", "DID_Male"),
         d %in% 24:25,
         event_time %in% 0:2) |>
  ggplot(aes(x = event_time, y = est, ymin = ci_l, ymax = ci_h,
             color = method, fill = method)) +
  geom_ribbon(color = NA, alpha = 0.2) +
  geom_point() + geom_line() +
  facet_grid(cols = vars(d), rows = vars(estimand), scales = "free") +
  labs(x = "Event Time", y = "Estimate +/- 95% CI",
       color = "Estimator", fill = "Estimator",
       title = "DID estimates by gender and treatment group") +
  theme(legend.position = "bottom")

## ----td_plot------------------------------------------------------------------
res |>
  filter(method == "TD",
         d %in% 24:25,
         event_time %in% 0:2) |>
  ggplot(aes(x = event_time, y = est, ymin = ci_l, ymax = ci_h)) +
  geom_ribbon(color = NA, alpha = 0.2, fill = "steelblue") +
  geom_point(color = "steelblue") + geom_line(color = "steelblue") +
  facet_grid(cols = vars(d), rows = vars(estimand), scales = "free") +
  labs(x = "Event Time", y = "Estimate +/- 95% CI",
       title = "TD: ATE(female) - ATE(male)") +
  theme(legend.position = "bottom")

## ----ntd_conv_plot------------------------------------------------------------
res |>
  filter(method == "NTD_Conv",
         d %in% 24:25,
         event_time %in% 0:2) |>
  ggplot(aes(x = event_time, y = est, ymin = ci_l, ymax = ci_h)) +
  geom_ribbon(color = NA, alpha = 0.2, fill = "darkorange") +
  geom_point(color = "darkorange") + geom_line(color = "darkorange") +
  facet_grid(cols = vars(d), rows = vars(estimand), scales = "free") +
  labs(x = "Event Time", y = "Estimate +/- 95% CI",
       title = expression(paste("NTD_Conv: ", theta[f], " - ", theta[m]))) +
  theme(legend.position = "bottom")

## ----ntd_new_plot-------------------------------------------------------------
res |>
  filter(method == "NTD_New",
         d %in% 24:25,
         event_time %in% 0:2) |>
  ggplot(aes(x = event_time, y = est, ymin = ci_l, ymax = ci_h)) +
  geom_ribbon(color = NA, alpha = 0.2, fill = "darkgreen") +
  geom_point(color = "darkgreen") + geom_line(color = "darkgreen") +
  facet_grid(cols = vars(d), rows = vars(estimand), scales = "free") +
  labs(x = "Event Time", y = "Estimate +/- 95% CI",
       title = expression(paste("NTD_New: ", Delta, rho,
                                " (effect on gender earnings ratio)"))) +
  theme(legend.position = "bottom")