```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
```

## Empirical Application: How the Abolition of Elected Councils Affects Local Public Services in Vietnam

Refer to Tuesday's lecture ([slides](causal3900.github.io/assets/slides/9-1_did.pdf) and recording) for some background on this experiment. Here's a [link](https://www.cambridge.org/core/journals/american-political-science-review/article/impact-of-recentralization-on-public-services-a-differenceindifferences-analysis-of-the-abolition-of-elected-councils-in-vietnam/3477854BAAFE152DC93C594169D64F58) to the paper from the original experiment and a [link](https://www.cambridge.org/core/journals/political-analysis/article/using-multiple-pretreatment-periods-to-improve-differenceindifferences-and-staggered-adoption-designs/747F618FD4AD82A536823521D89310F7) to the paper that re-analyzed this experiment, which we're basing this lab on. 

### Download the data
```{r}
install <- function(package) {
  if (!require(package, quietly = TRUE, character.only = TRUE)) {
    install.packages(package, repos = "http://cran.us.r-project.org", type = "binary")
    library(package, character.only = TRUE)
  }
}

install("cowplot")
install("dplyr")
install("tidyverse")

malesky2014 <- read_csv("https://github.com/causal3900/causal3900.github.io/raw/main/assets/data/malesky2014.csv")
malesky2014 <- drop_na(malesky2014, lnarea, lnpopden, city) # data cleaning
```


### About the data

Variables in the data set:

- `treatment`: binary variable indicating if the record is from the treatment group ($1$) or the control group ($0$)
- `pro4`: educational and cultural program number 4 outcome variable (binary: $1$ if they participated in the program, $0$ if not)
- `tapwater`: tap water outcome variable (binary: $1$ if they drink tap water, $0$ if not)
- `agrext`: agricultural center outcome variable (binary: $1$ if they use the center, $0$ if not)
- `year`: the year the record is from
- `post_treat`: binary variable indicating if the record is from the pre-treatment period ($0$) or the post-treatment period ($0$)

```{r}
head(malesky2014)
```


### Visualizing Trends in the Treatment and Control Groups

In a DID design, the first step is to visualize the trends of the treatment and control groups to check if the parallel trends assumption is credible. Here, we look at three possible outcomes: "Education and Cultural Program," "Tap Water," and "Agricultural Center."

```{r}
# Plot the outcome for `Education and Cultural Program`
ed_and_culture <- malesky2014 %>%
  mutate(treatment = ifelse(treatment == 1, "Treated in 2009","Untreated")) %>%
  group_by(treatment, year) %>%
  summarize(estimate = mean(pro4)) %>%
  ggplot(aes(x = year, y = estimate, color = treatment)) + 
  geom_point() + 
  geom_line() + 
  ggtitle("Ed. and Cultural Program") +
  theme_minimal() +
  theme(legend.position = "none",
        legend.background = element_rect(fill = "white", colour = "white"),
        plot.title = element_text(hjust = 0.5))

# Plot the outcome for `Tap Water`
tap_water <- malesky2014 %>%
  mutate(treatment = ifelse(treatment == 1, "Treated in 2009","Untreated")) %>%
  group_by(treatment, year) %>%
  summarize(estimate = mean(tapwater)) %>%
  ggplot(aes(x = year, y = estimate, color = treatment)) +
  geom_point() +
  geom_line() +
  ggtitle("Tap Water") +
  theme_minimal() +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

# Plot the outcome for `Agricultural Center`
ag_center <- malesky2014 %>%
  mutate(treatment = ifelse(treatment == 1, "Treated in 2009","Untreated")) %>%
  group_by(treatment, year) %>%
  summarize(estimate = mean(agrext)) %>%
  ggplot(aes(x = year, y = estimate, color = treatment)) +
  geom_point() +
  geom_line() +
  ggtitle("Agricultural Center") +
  theme_minimal() + 
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))
```

Let's compare all three of these plots side-by-side!
```{r}
legend <- get_legend(
  ed_and_culture + 
    guides(color = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom")
)
pgrid <- cowplot::plot_grid(ed_and_culture, tap_water, ag_center, nrow = 1)
plot_grid(pgrid, legend, ncol = 1, rel_heights = c(1, .1))
```
Notice how we see parallel trends in the pretreatment period for the outcome "Education and Cultural Program," but not for "Tap Water" or "Agricultural Center."

### Basic Difference in Differences Analysis with One Pretreatment Period

In the data our treatment indicator variable is called `treatment`, and our
post-treatment time period indicator is called `post_treat`. You will need to
create the interaction variable that is `treatment`*`post_treat`. Call this
new variable `did`.

```{r}
# remove the extra pretreatment periods from the dataset
onePerMalesky <- malesky2014 %>%
  filter(year==2008 | year ==2010)

# Now, create the interaction term between treatment and post_treat
onePerMalesky <- onePerMalesky %>%
  mutate(did = ...create your interaction variable here...)
```

Now that you have your `treatment`, `post_treat`, and `did` variables, you are
ready to estimate the causal effect of de-centralizing government on whether or
not people participate in an Educational and Cultural Program. Your outcome
variable in this case is `pro4`. We will estimate the causal effect with an 
OLS regression. Use the formula that we discussed in the slides to create 
your model formula.
```{r}
# regression for Education and Cultural Program on treatment, time period, and interaction term
didReg <- lm(...put your formula here..., data = onePerMalesky)

summary(...put your DiD model object here...)

# estimate
print(didReg$coefficients["did"])
```


