############### ###############
## F_Est_H5
## Project: CBO
## PI: Malte Lierl
## Written by: Malte Lierl
## Date of creation: 2022-05-22
############### ###############

## REQUIRES USER-WRITTEN FUNCTIONS: 
# scale.control()

## OUTPUT PATH
path <- here("Output/For Publication/Figures/F_Est_")

## LOOP OVER PERFORMANCE INDICATORS AND ESTIMATE STANDARDIZED AVERAGE TREATMENT EFFECT ON EACH

results <- NULL
vars <- NULL

for (i in indicators$indicator[1:18]) {
  fit.i <- lm_robust(scale.control(eval(parse(text=paste0(i))), treat=cboincentives) ~ cboincentives + 
                       cboincentives*scale.group(eval(parse(text=paste0(i,"_lag2"))), treat=cboincentives), 
                     data=SUPERMUN_final %>% filter(year==2020  & 
                                                municipality_included==1 & 
                                                is.na(eval(parse(text=paste0(i))))==F &
                                                is.na(eval(parse(text=paste0(i,"_lag2"))))==F), fixed_effects=randomization_block_cbo)
  out.i <- summary(fit.i)$coefficients["cboincentives",]
  results <- rbind(results, out.i)
  vars.i <- c(indicators$caption_english[indicators$indicator==i],i)
  vars <- rbind(vars, vars.i)
}

results <- data.frame(vars, results, row.names=NULL)
names(results) <- c("caption", "varname", "estimate", "se", "tstat", "pval", "lower95", "upper95", "df")
results

results.reordered <- results[results$varname%in%c("total_points_ic", "total_points_sd")==FALSE,]
results.reordered <- results.reordered[order(results.reordered$estimate, decreasing=TRUE, na.last=NA),]
results.reordered <- rbind(results[results$varname%in%c("total_points_ic", "total_points_sd"),], results.reordered)
results.reordered$y <- nrow(results.reordered) - 1:nrow(results.reordered) + 1
results.reordered$caption.ordered <- ordered(results.reordered$caption, levels=results.reordered$y, labels=results.reordered$caption)


## VISUALIZE AVERAGE TREATMENT EFFECTS ON DIFFERENT PERFORMANCE INDICATORS
results.viz <- ggplot(data=results.reordered) + 
  geom_vline(aes(xintercept=0), lwd=0.1, lty=2, alpha=0.5) +
  geom_errorbarh(aes(xmin=lower95, xmax=upper95, y=y), height=0, position=position_nudge(y=-0), data=results.reordered) + 
  #geom_point(aes(x=estimate, y=caption), color=gray50, alpha=0.5)+
  geom_point(aes(x=estimate, y=y), position=position_nudge(y=-0), data=results.reordered)+
  ylab("") + xlab("Average Treatment Effect (Standardized, 95% CI)") + 
  coord_cartesian(xlim=c(-1,1))+
  scale_y_continuous(element_blank(), breaks=results.reordered$y, labels=results.reordered$caption, minor_breaks= results.reordered$y)+ 
  theme_light()
results.viz
ggsave(file=paste0(path,"H5.png"), plot=results.viz, width=7, height=5, dpi=200)

