Remember to answer all the questions in R markdown and produce a PDF. Email your completed homework (R markdown file and PDF) to \href{jonnyphillips@gmail.com}{jonnyphillips@gmail.com} by midnight the night before class. Remember to refer to the example code from this week and the last couple of weeks for coding guidance.

Load the Boas and Hidalgo dataset `combined_data.RDS` using `readRDS()`. There are many tools and settings that can be adjusted when using matching - we will try and repicate the Boas and Hidalgo results using more user-friendly tools than in their own code, so we will not recover identical values. 

```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
#rm(list=ls())
knitr::opts_chunk$set(include=TRUE, echo=FALSE, warning=FALSE, message=FALSE)
library(foreign)
library(tidyverse)
library(knitr)
library(broom)
library(stargazer)
library(xtable)
library(Zelig)
library(ZeligChoice)
```

```{r}
#Pre-process for sharable dataset
#load("replication.RData")
#combined_data <- cbind(match.data[,c("treat","pctVV","elected","male")],genmatch.covar)
#saveRDS(combined_data,"combined_data.RDS")
d <- readRDS("combined_data.RDS")
```

**1. The treatment variable in Boas and Hidalgo is whether a councillor that applied for a media licence received approval before the 2004 election (`treat`). The outcome variable is the councillor's vote share in the 2004 elections (`pctVV`). Conduct and interpret a basic linear regression of the outcome on treatment with no controls.**

```{r, results='asis'}
Q1 <- d %>% zelig(pctVV ~ treat, data=., model="ls",cite=FALSE)

Q1 %>% from_zelig_model() %>%
  stargazer(header=FALSE,title="Q1")
```

Councillors that receive approved media licences before the 2004 elections are associated with a `r Q1 %>% from_zelig_model() %>% tidy() %>% filter(term=="treat") %>% select(estimate) %>% round(3)` % points increase in vote share in the 2004 election. However, this effect is not causal.

**2. One potential confounding variable is gender (this could affect the chances of an application being approved if there is bias in the Ministry, and the candidate's vote share if there is bias among voters). Is there balance across control and treatment groups on the `male` variable?**

```{r}
Q2_diff <- d %>% group_by(treat) %>% dplyr::summarise(male=mean(male,na.rm=TRUE)) %>% dplyr::summarise(diff=diff(male)) %>% `*` (100) %>% round(3)

Q2_p.value <- d %>% summarise_at(vars(male),funs(t.test(.[treat==0], .[treat==1])$p.value)) %>% round(3)

```

No. The treatment group has `r Q2_diff`% points more men than the control group, which is a statistically significant difference with a p-value of `r Q2_p.value`.

**3. One 'control for confouding' approach is to use a regression. Conduct the regression of the oucome (`pctVV`) on treatment (`treat`), controlling for `male` and compare the results to your answer to Q1.**

```{r, results='asis'}
d %>% zelig(pctVV ~ treat + male, data=., model="ls",cite=FALSE) %>%
  from_zelig_model() %>%
  stargazer(header=FALSE,title="Q3")                  
```

The results slightly reduce the coefficient estimate on the treatment variable, but it remains significant. Being male is positively associated with the outcome, but the difference is not statistically significant.

**4. An alternative approach to overcome confounding is matching: removing units from the dataset until we have balance on gender. Let's try to implement one-to-one exact matching MANUALLY to illustrate the process:**  
a. Create separate datasets for treated and control units.
b. Create a new column in the treated dataset (eg. matched_control_unit_pctVV) where we are going to save the outcome for the matched control unit.
b. Create a `for` loop so we can analyze each *treated* unit (each row of your treated dataset) separately.  
c. Within your for loop, for each treated unit, identify its gender, and filter the control units dataset to only include that same gender.  
d. There are more control units than treated units, so we can pick any same-gender control unit as our matched unit. *Hint:* Use `sample_n(1)` to randomly pick one row of your control unit data.  
e. Now select only the outcome variable for this matched control unit and save that value to the new column you made in the *treated* dataset (eg. matched_control_unit_pctVV). Make sure you save it only to the row of the current treated unit you are analysing in your `for` loop.

```{r}
d_treat <- d %>% filter(treat==1)
d_control <- d %>% filter(treat==0)
d_treat[,"matched_control_unit_pctVV"] <- NA

for (i in 1:dim(d_treat)[1]){
  gender_treated <- d_treat[i,] %>% select(male) %>% as.numeric()
  matched_unit <- d_control %>% filter(male==gender_treated) %>% sample_n(1)
  d_treat[i,"matched_control_unit_pctVV"] <- matched_unit %>% select(pctVV) %>% as.numeric()
}

```

**5. With this matched dataset, we now have balance on gender by definition. To estimate the treatment effect using a difference-in-means calculate the difference in outcomes between your treated and matched control units in the treated units dataset (i.e. between the `pctVV` and `matched_control_unit_pctVV` variables). Then average these differences to get an average treatment effect and interpret the result. Also conduct a t-test to get a p-value.**

```{r}
Q5_diff <- d_treat %>% mutate(diff_outcome=pctVV-matched_control_unit_pctVV) %>% summarise(diff_outcome=mean(diff_outcome,na.rm=TRUE)) %>%
  round(3)

Q5_p.value <- t.test(d_treat$pctVV,d_treat$matched_control_unit_pctVV)$p.value %>% round(3)

```

The estimated average treatment effect for receiving an approved media licence before the election is a `r Q5_diff`% points increase in vote share, which is statistically significant with a p-value of `r Q5_p.value`.

**6. We can also use post-matching regression to do the analysis stage. To do so, we need to reorganize our dataset to include rows for both the treated and matched control units, with their corresponding outcomes.**  
a. Select only the 'pctVV' and 'matched_control_unit_pctVV' columns.
b. Transform the data from wide to long format using `gather`. Your dataset should now have 622 rows.
c. Use `mutate` with `ifelse` to rename the treatment column values from 'pctVV' and 'matched_control_unit_pctVV' to 1 and 0 respectively.
d. Conduct the simple linear regression of the outcome on treatment for this dataset.  
e. Compare the results to your answer in Q5.  

```{r, results='asis'}
Q6 <- d_treat %>% select(pctVV,matched_control_unit_pctVV) %>% gather(key="Treatment",value="pctVV_all") %>% mutate(Treatment=ifelse(Treatment=="pctVV",1,0))

Q6_reg <- Q6 %>% zelig(pctVV_all ~ Treatment,data=.,model="ls",cite=FALSE) 

Q6_reg %>%
  from_zelig_model() %>%
  stargazer(header=FALSE,title="Q6")
```

The estimate treatment effect is `r Q6_reg %>% from_zelig_model() %>% tidy() %>% filter(term=="Treatment") %>% select(estimate) %>% round(3)`, with a p-value of `r Q6_reg %>% from_zelig_model() %>% tidy() %>% filter(term=="Treatment") %>% select(p.value) %>% round(3)`. This is exactly the same as in Q5.

**7. To match on continuous variables, we can't match 'exactly' so we do 'nearest' neighbour matching. Another potential confounder is  the size of the electorate (`log.valid.votes`). Use `matchit` to construct a datset of treated and control units matched on (only) the size of the electorate with nearest neighbour matching. In your matched data, how many treated and control units are there? *Hint:* Access the units summary using `output$nn`, where `output` is the result of `matchit`.**

```{r, results='asis'}
library(MatchIt)
matched_data1 <- matchit(treat~log.valid.votes, data=d,method="nearest")
matched_data1$nn %>% kable(caption="Q7")

```

**8. Check that this matching has improved balance by comparing the average difference in `log.valid.votes` in the original dataset and in the matched dataset. *Hint:* Remember you can use `match.data(output_of_matchit)` to get the matched dataset.**

```{r, echo=FALSE, warning=FALSE, message=FALSE, results='asis'}
Q8_full <- d %>% group_by(treat) %>% 
  dplyr::summarize(log.valid.votes=mean(log.valid.votes))
  
Q8_matched <- match.data(matched_data1) %>% group_by(treat) %>% dplyr::summarize(log.valid.votes=mean(log.valid.votes))

Q8_full %>% left_join(Q8_matched,by="treat") %>% dplyr::rename(Before_Matching=log.valid.votes.x,After_Matching=log.valid.votes.y) %>% xtable() %>% print(comment=FALSE)
```

Balance improves considerably after matching.

\newpage

**9. To see which units have been matched, plot two layers on a graph of the confounder `log.valid.votes` against the outcome `pctVV`. One layer should be the full dataset - make this a specific size, shape and alpha (transparency). The second layer should be the matched dataset from Q8 - make this a different size, shape and alpha so that the matched units stand out. Group and colour each layer by treatment status.**

```{r}
d %>% ggplot() + geom_point(aes(x=log.valid.votes,y=pctVV,group=factor(treat),color=factor(treat)),shape=1,alpha=0.4) +
  geom_point(data=match.data(matched_data1),aes(x=log.valid.votes,y=pctVV,group=factor(treat),color=factor(treat)),shape=2)
```


**10. Use a simple least squares regression to estimate the effect of the treatment ('treat') on the outcome ('pctVV') in the matched dataset. Interpret the result.**

```{r, results='asis'}
Q10 <- match.data(matched_data1) %>% zelig(pctVV~treat, data=., model="ls",cite=FALSE)

Q10 %>%
  from_zelig_model() %>%
  stargazer(header=FALSE,title="Q10")
```

Receiving an approved media licence application before the election is associated with a `r Q10 %>% from_zelig_model() %>% tidy() %>% filter(term=="treat") %>% dplyr::select(estimate) %>% round(3)` % points increase in vote share, but this is not statistically significant.

**11. Conduct the same regression as in Q8 but add a control for `log.valid.votes`, the variable we already matched on. Does this change the estimate of the treatment effect compared to your answer in Q8? Why/why not?**

```{r,results='asis'}
Q11 <- match.data(matched_data1) %>% zelig(pctVV~treat + log.valid.votes, data=., model="ls",cite=FALSE)
#Check if need match.data here...
Q11 %>%
  from_zelig_model() %>%
  stargazer(header=FALSE,title="Q11")
```

The coefficient on treatment is essentially unchanged because we have already gone a long way towards achieving balance through the matching procedure, and the regression only controls for any residual imbalance. The p-value is slightly smaller because we have absorbed some of the variation with the control, allowing greater precision.

**12. Now lets include all of the matching variables that Boas and Hidalgo use, and use nearest neighbour matching to construct a matched dataset. Use the list of matching variables provided below for matching. **

```{r, echo=TRUE,tidy=TRUE}
covars <- c("occBlue.collar", "occEducation", "occGovernment", "occMedia", "occNone", "occOther", "occPolitician", "occWhite.collar", "lat", "long", "ran.prior", "incumbent", "log.valid.votes", "party.prior.pctVV", "prior.pctVV", "elec.year", "match.partyPCB", "match.partyPC.do.B", "match.partyPDT", "match.partyPFL", "match.partyPL", "match.partyPMDB", "match.partyPMN", "match.partyPP", "match.partyPPS", "match.partyPSB", "match.partyPSC", "match.partyPSDB", "match.partyPSDC", "match.partyPSL", "match.partyPT", "match.partyPTB", "match.partyPV", "uf.rs", "uf.sp", "yob", "eduMore.than.Primary..Less.than.Superior", "eduSome.Superior.or.More", "log.total.assets", "pt_pres_1998", "psdb_2000", "hdi_2000", "income_2000", "log.num.apps")

covars_formula <- paste0(covars,collapse=" + ")
```

```{r}
matched_data1 <- matchit(as.formula(paste0("treat~",covars_formula)), data=d,method="nearest")

#match.data(matched_data1)

#Include matchit plots?
```

**13. For your matched dataset from Q12, report a before-matching and after-matching balance table similar to the first four columns of Boas and Hidalgo's Table 3. The standardized difference in means they use allows us to compare everything on the same scale and requires calculating the difference in means and then dividing by the standard deviation of the treated group:**  
a. Calculate the standard deviation of the treated variable for the covariates. *Hint:* `summarise_at` is useful to apply the same summary to the list of covariates provided above.
b. Summarise the average of each covariate separately for the treated and control groups, and then calculate the difference.   
c. Divide this difference in means by the standard deviation of the treated group for each variable, and multiply by 100.  
d. Calulate the p-value of the difference in means for each covariate between treated and control groups. *Hint:* Try something like `d %>%  summarise_at(covars,funs(t.test(.[treat==0], .[treat==1])$p.value))`.   
e. Now combine the columns for standardized differences in means and p-values for the full and matched datasets into a single table. *Hint:* Use `gather()` to turn wide rows into long columns, and `left_join` to combine the columns into a single table.  

```{r, results='asis'}
treat_sds <- d %>% filter(treat==1) %>% summarise_at(covars,sd)

fulldata_mean_diff <- d %>% group_by(treat) %>% summarise_at(covars,mean) %>% summarise_all(diff) %>% dplyr::select(-treat)
fulldata_mean_diff_std <- (fulldata_mean_diff/treat_sds)*100 
fulldata_mean_diff_std <- fulldata_mean_diff_std %>% gather()

fulldata_p <- d %>%  summarise_at(covars,funs(t.test(.[treat==0], .[treat==1])$p.value)) %>% gather()

Q13_data <- match.data(matched_data1)

treat_sds_matched <- Q13_data %>% filter(treat==1) %>% summarise_at(covars,sd)

matcheddata_mean_diff <- Q13_data %>% group_by(treat) %>% summarise_at(covars,mean) %>% summarise_all(diff) %>% dplyr::select(-treat)
matcheddata_mean_diff_std <- (matcheddata_mean_diff/treat_sds_matched)*100
matcheddata_mean_diff_std <- matcheddata_mean_diff_std %>% gather()

matched_p <- Q13_data %>% summarise_at(covars,funs(t.test(.[treat==0], .[treat==1])$p.value)) %>% gather()

fulldata_mean_diff_std %>% left_join(matcheddata_mean_diff_std, by=c("key")) %>% 
  left_join(fulldata_p,by=c("key")) %>%
  left_join(matched_p,by=c("key")) %>%
  dplyr::rename(Before_Mean=value.x, After_Mean=value.y, Before_p=value.x.x, After_p=value.y.y) %>% 
  mutate_at(vars(Before_Mean,After_Mean,Before_p,After_p),funs(round(.,3))) %>%
  xtable(caption="Q13") %>% print(comment=FALSE)
```

**14. Implement the regression of the outcome on treatment (excluding controls) on the matched dataset and compare your answer to that of Boas and Hidalgo in Table 4.**

```{r, results='asis'}
Q14 <- Q13_data %>% zelig(pctVV ~ treat, data=., model="ls", cite=FALSE) 

Q14 %>%
  from_zelig_model() %>%
  stargazer(header=FALSE,title="Q14")
```

The estimate is considerably smaller than Boas and Hidalgo's at `r Q14 %>% from_zelig_model() %>% tidy() %>% filter(term=="treat") %>% dplyr::select(estimate) %>% round(3)`, and is statistically insignificant with a p-value of `r Q14 %>% from_zelig_model() %>% tidy() %>% filter(term=="treat") %>% dplyr::select(p.value) %>% round(3)`

**15. One risk with the nearest-neighbour method is that the control unit can still be far away from the treated unit if there are no good matches, so imbalance can still be a problem. We can put a 'caliper' on the matching method so that it only includes 'close' matches within a specific number of standard deviations. Re-run the matching process from Q12 but with a caliper of 0.01 standard deviations, and then conduct the regression from Q14. Note how the number of units changes and interpret if the result is any different.**

```{r, echo=FALSE, warning=FALSE, message=FALSE, results='asis'}
Q15_data <- matchit(as.formula(paste0("treat~",covars_formula)), data=d,method="nearest", caliper = 0.01)
Q15 <- match.data(Q15_data) %>% zelig(pctVV ~ treat, data=., model="ls", cite=FALSE) 

Q15 %>% 
  from_zelig_model() %>% 
  stargazer(header=FALSE,title="Q15")
```

The coefficient estimate is higher at `r Q15 %>% from_zelig_model() %>% tidy() %>% filter(term=="treat") %>% dplyr::select(estimate) %>% round(3)`, and has a smaller p-value of `r Q15 %>% from_zelig_model() %>% tidy() %>% filter(term=="treat") %>% dplyr::select(p.value) %>% round(3)`. The numer of units reduces from 622 to `r Q15_data$nn %>% as.data.frame() %>% rownames_to_column() %>% filter(rowname=="Matched") %>% mutate(N=Control+Treated) %>% dplyr::select(N)` as treated units with no matching control units within the caliper distance are discarded.

**16. Boas and Hidalgo use genetic matching, which is a complex automated process. We can implement genetic matching easily by loading the package `rgenoud` and choosing `method="genetic"` in the matching process in `matchit`. Then run the same regression again and compare the results to those in Q14 and in Boas and Hidalgo. *Hint:* Genetic matching might take 10-20 minutes, so in your R chunk options set `cache=TRUE` so that it doesn't re-run every time you knit to PDF.**

```{r,cache=TRUE, echo=FALSE, warning=FALSE, message=FALSE, results='hide', cache=TRUE}
library(rgenoud)
Q16_data <- matchit(as.formula(paste0("treat~",covars_formula)), data=d,method="genetic")


```

```{r, echo=FALSE, warning=FALSE, message=FALSE, results='asis'}
Q16 <- match.data(Q16_data) %>% zelig(pctVV ~ treat, data=., model="ls", cite=FALSE) 

Q16 %>% 
  from_zelig_model() %>% 
  stargazer(header=FALSE,title="Q16")
```

The coefficient estimate is higher at `r Q16 %>% from_zelig_model() %>% tidy() %>% filter(term=="treat") %>% dplyr::select(estimate) %>% round(3)`, and has a smaller p-value of `r Q16 %>% from_zelig_model() %>% tidy() %>% filter(term=="treat") %>% dplyr::select(p.value) %>% round(3)`. However, the esimate and p-value are still a little distant from the Boas and Hidalgo results. 

**17. As with regression, these estimates only have a causal interpretation if matching is successful at producing balance on ALL of the confounding covariates. Provide an example of a confounding variable which would bias Boas and Hidalgo's estimates despite their matching procedure.**

A confounding variable could be anything that is not included in the matching procedure (or anything in the matching procedure which remains imbalanced). For example, they do not have a perfect measure of the councillors' personal wealth; if bribes or status affect the chance of getting an application approved and also boost electoral success, the results are not causal estimates.