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

## REQUIRES USER-WRITTEN FUNCTIONS: 
# execute()
# scale.control()
# scale.group()
# tex()
# stars()
# se()
# obs()

##### ESTIMATION BEGINS HERE ####

## LAGGED DEPENDENT VARIABLE SPECIFICATION WITH RANDOMIZATION BLOCK FIXED EFFECTS

fit.1 <- summary(lm_robust(scale.control(total_points_all, treat=cboincentives) ~ 
                             cboincentives*scale.group(conflict_fatalities_lag2, treat=cboincentives) +
                             cboincentives*scale.group(total_points_all_lag2, treat=cboincentives), 
                           data=SUPERMUN_final %>% filter(year==2020 & 
                                                      municipality_included==1 & 
                                                      is.na(total_points_all)==F & 
                                                      is.na(conflict_fatalities_lag2)==F & 
                                                      is.na(total_points_all_lag2)==F), 
                           fixed_effects=factor(randomization_block_cbo)))
fit.1

fit.2 <- summary(lm_robust(scale.control(total_points_all, treat=cboincentives) ~ 
                             cboincentives*scale.group(conflict_fatalities_lag2, treat=cboincentives) + 
                             cboincentives*scale.group(total_points_all_imp_lag2, treat=cboincentives), 
                           data=SUPERMUN_final %>% filter(year==2020 & 
                                                      municipality_included==1 & 
                                                      is.na(total_points_all)==F & 
                                                      is.na(conflict_fatalities_lag2)==F & 
                                                      is.na(total_points_all_imp_lag2)==F), 
                           fixed_effects=factor(randomization_block_cbo)))
fit.2

fit.3 <- summary(lm_robust(scale.control(total_points_all_imp, treat=cboincentives) ~ 
                             cboincentives*scale.group(conflict_fatalities_lag2, treat=cboincentives) + 
                             cboincentives*scale.group(total_points_all_imp_lag2, treat=cboincentives), 
                           data=SUPERMUN_final %>% filter(year==2020 & 
                                                      municipality_included==1 & 
                                                      is.na(total_points_all_imp)==F & 
                                                      is.na(conflict_fatalities_lag2)==F & 
                                                      is.na(total_points_all_imp_lag2)==F), 
                           fixed_effects=factor(randomization_block_cbo)))
fit.3

fit.sd <- summary(lm_robust(scale.control(total_points_sd, treat=cboincentives) ~ 
                              cboincentives*scale.group(total_points_sd_imp_lag2, treat=cboincentives), 
                            data=SUPERMUN_final %>% filter(year==2020 & 
                                                       municipality_included==1 & 
                                                       is.na(total_points_sd)==F & 
                                                       is.na(total_points_sd_imp_lag2)==F), 
                            fixed_effects=factor(randomization_block_cbo)))
fit.sd

fit.ic <- summary(lm_robust(scale.control(total_points_ic, treat=cboincentives) ~ 
                              cboincentives*scale.group(total_points_ic_imp_lag2, treat=cboincentives), 
                            data=SUPERMUN_final %>% filter(year==2020 & 
                                                       municipality_included==1 & 
                                                       is.na(total_points_ic)==F & 
                                                       is.na(total_points_ic_imp_lag2)==F), 
                            fixed_effects=factor(randomization_block_cbo)))
fit.ic


## Endline-only specification, without adjusting for prior-year values would be least biased, but much less efficient. 
summary(lm_robust(total_points_all ~ cboincentives, data=SUPERMUN_final %>% filter(year==2020 & municipality_included==1), fixed_effects=factor(randomization_block_cbo)))
summary(lm_robust(total_points_all_imp ~ cboincentives, data=SUPERMUN_final %>% filter(year==2020 & municipality_included==1), fixed_effects=factor(randomization_block_cbo)))
# --> QUESTION: Should we report this specification too, at the risk of being under-powered, because we don't include baseline values? 


###### TABLE OUTPUT BEGINS HERE #######

## For each indicator, show estimated treatment effect and 95% CI

tex(code="
  \\begin{tabularx}{\\textwidth}{l X X X }
  \\hline
  & \\multicolumn{3}{c}{\\emph{DV: SUPERMUN municipal performance score (standardized)}}\\\\
  \\\\
  & (1) & (2) & (3)  \\\\
  &  No adjustment & Missing baseline values imputed & Missing baseline and endline values imputed\\\\
  \\hline 
  \\\\
  Standardized ATE  & <<stars(fit.1,'cboincentives',2)>>  & <<stars(fit.2,'cboincentives',2)>>  & <<stars(fit.3,'cboincentives',2)>> \\\\
                    & (<<se(fit.1,'cboincentives',2)>>)   & (<<se(fit.2,'cboincentives',2)>>)   & (<<se(fit.3,'cboincentives',2)>>)\\\\
  \\\\
  Baseline performance  & <<stars(fit.1,'scale.group(total_points_all_lag2, treat = cboincentives)',2)>>  & <<stars(fit.2,'scale.group(total_points_all_imp_lag2, treat = cboincentives)',2)>>  & <<stars(fit.3,'scale.group(total_points_all_imp_lag2, treat = cboincentives)',2)>> \\\\
                    & (<<se(fit.1,'scale.group(total_points_all_lag2, treat = cboincentives)',2)>>)   & (<<se(fit.2,'scale.group(total_points_all_imp_lag2, treat = cboincentives)',2)>>)   & (<<se(fit.3,'scale.group(total_points_all_imp_lag2, treat = cboincentives)',2)>>)\\\\
  \\\\
  Treatment $\\times$ Baseline performance  & <<stars(fit.1,'cboincentives:scale.group(total_points_all_lag2, treat = cboincentives)',2)>>  & <<stars(fit.2,'cboincentives:scale.group(total_points_all_imp_lag2, treat = cboincentives)',2)>>  & <<stars(fit.3,'cboincentives:scale.group(total_points_all_imp_lag2, treat = cboincentives)',2)>> \\\\
                    & (<<se(fit.1,'cboincentives:scale.group(total_points_all_lag2, treat = cboincentives)',2)>>)   & (<<se(fit.2,'cboincentives:scale.group(total_points_all_imp_lag2, treat = cboincentives)',2)>>)   & (<<se(fit.3,'cboincentives:scale.group(total_points_all_imp_lag2, treat = cboincentives)',2)>>)\\\\
  \\\\
  Baseline conflict exposure  & <<stars(fit.1,'scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>  & <<stars(fit.2,'scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>  & <<stars(fit.3,'scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>> \\\\
                    & (<<se(fit.1,'scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>)   & (<<se(fit.2,'scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>)   & (<<se(fit.3,'scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>)\\\\
  \\\\
  Treatment $\\times$ Baseline conflict exposure  & <<stars(fit.1,'cboincentives:scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>  & <<stars(fit.2,'cboincentives:scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>  & <<stars(fit.3,'cboincentives:scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>> \\\\
                    & (<<se(fit.1,'cboincentives:scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>)   & (<<se(fit.2,'cboincentives:scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>)   & (<<se(fit.3,'cboincentives:scale.group(conflict_fatalities_lag2, treat = cboincentives)',2)>>)\\\\
  \\\\
  \\hline
  \\\\
  Observations (Municipalities)     & <<fit.1$nobs>>                                   & <<fit.2$nobs>>                                   & <<fit.3$nobs>>\\\\
  \\hline\\hline
  \\end{tabularx}
  \\footnotesize{All effects are estimated via OLS regressions with fixed effects for randomization blocks. In columns (2) and (3), missing values on individual performance indicators have been replaced with prior-year values and used together with the other, available indicators to calculate the aggregate performance score. 
  HC2 standard errors in parentheses. $^{*} p<0.1$, $^{**} p<0.05$, $^{***} p<0.01$ (two-sided). }", 
    file=paste0(here("Output", "For Publication", "Tables","T_H5.tex")))

    