treatment_effects.Rmd

---
title: "Treatment effects"
output:
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
set.seed(1)
source("install_packages.r")
require(knitr)
require(dplyr)
require(plotly)
require(ggplot2)
require(gridExtra)
require(MatchIt)
require(ATE)
```



## Data

```{r}
dta <- read.csv("dta/TableF19-3.csv")
```

LaLonde (1986) Earnings Data, 2,490 Control Observations and 185 Treatment Observations 
Source: LaLonde (1986).

| Variable | Description                              |   |
|----------|------------------------------------------|---|
| t        | Treatment dummy variable,                |   |
| age      | Age in years,                            |   |
| educ     | Education in years,                      |   |
| marr     | Dummy variable for married,              |   |
| black    | Dummy variable for black,                |   |
| hisp     | Dummy variable for Hispanic,             |   |
| nodegree | Dummy for no degree (not used),          |   |
| re74     | Real earnings in 1974, scaled by 10,000, |   |
| re75     | Real earnings in 1975, scaled by 10,000, |   |
| re78     | Real earnings in 1978, scaled by 10,000, |   |
| age2     | Age squared,                             |   |
| educ2    | Educ squared,                            |   |
| re742    | Re74 squared                             |   |
| re752    | Re75 squared                             |   |
| blacku74 | Black times 1(re74                       |   |

## Data exploration 

```{r}
summary(dta)
sum(dta$T)
```

```{r, echo=FALSE}
grid.arrange(
  ggplot(dta) + geom_boxplot(aes(as.factor(T),AGE)),
  ggplot(dta) + geom_boxplot(aes(as.factor(T),EDUC)),
  ncol=2
  )

grid.arrange(
  ggplot(dta[dta$T==0,]) +
    geom_bar(aes(x=factor(1), fill=as.factor(BLACK))) + 
    coord_polar(theta = "y") + 
    guides(fill=guide_legend(title="BLACK")),
  ggplot(dta[dta$T==1,]) +
    geom_bar(aes(x=factor(1), fill=as.factor(BLACK))) +
    coord_polar(theta = "y") +
    guides(fill=guide_legend(title="BLACK")),
  ncol=2
)
grid.arrange(
  ggplot(dta[dta$T==0,]) +
    geom_bar(aes(x=factor(1), fill=as.factor(HISP))) + 
    coord_polar(theta = "y") + 
    guides(fill=guide_legend(title="HISP"))+ylab("untreated"),
  ggplot(dta[dta$T==1,]) +
    geom_bar(aes(x=factor(1), fill=as.factor(HISP))) +
    coord_polar(theta = "y") +
    guides(fill=guide_legend(title="HISP"))+ylab("treated"),
  ncol=2
)
```



```{r, echo=FALSE}
g <- ggplot(dta) + ylim(0,16) + xlab("Treatment (1,0)")
grid.arrange(
  g + geom_boxplot(aes(as.factor(T),RE74)),
  g + geom_boxplot(aes(as.factor(T),RE75)),
  g + geom_boxplot(aes(as.factor(T),RE78)),
  ncol=3
  )
g0 <- ggplot(dta[dta$T==0,]) +ylim(0,16) + xlab("Treatment (1,0)")
g1 <- ggplot(dta[dta$T==1,])+ ylim(0,7)  + xlab("Treatment (1,0)")
grid.arrange(
  g0 + geom_boxplot(aes(as.factor(T),RE74)),
  g0 + geom_boxplot(aes(as.factor(T),RE75)),
  g0 + geom_boxplot(aes(as.factor(T),RE78)),
  ncol=3)
grid.arrange(
  g1 + geom_boxplot(aes(as.factor(T),RE74)),
  g1 + geom_boxplot(aes(as.factor(T),RE75)),
  g1 + geom_boxplot(aes(as.factor(T),RE78)),
  ncol=3
  )
```



```{r}
dta %>% group_by(T) %>% summarise(mean=mean(RE78),sd=sd(RE78))
```
```{r}
t.test(dta$RE78~dta$T)
```


## Difference in means: pre treatment cov

```{r}
dta %>% group_by(T) %>% select(AGE,EDUC,MARR,BLACK,HISP,NODEGREE) %>% summarise_all(funs(mean(.)))
```
```{r}
lapply(c("AGE","EDUC","MARR","BLACK","HISP","NODEGREE"), function(x){t.test(dta[,x]~dta[,"T"])})
```

## Propensity score estimation


```{r}
summary(glm1 <- glm(T~.-RE78,   family = binomial(), data=dta))
```

```{r}
prs_df <- data.frame(pr_score = predict(glm1, type = "response"),
                     T = glm1$model$T)
sample_n(prs_df, size=5)
```

```{r, echo=FALSE}
grid.arrange(
  ggplot(prs_df[prs_df$T==0,])+geom_histogram(aes(x=pr_score))+ggtitle(""),
  ggplot(prs_df[prs_df$T==1,])+geom_histogram(aes(x=pr_score)),
  ggplot(prs_df)+
    geom_density(aes(x=pr_score, 
                     group=as.factor(T),
                     color=as.factor(T), 
                     fill=as.factor(T)),
                 alpha = 0.2)+
    xlab("Probability to be part of treatment group")
)
```