############### ###############
## F_Est_H3_summary
## Project: CBO
## Author: Kamil Kouhen
## Date of creation: 23/06/2022
############### ###############

### Creating uniform path
path <- here("Output/For Publication/Figures/F_Est_")

# Hypothesis 3 summary --------------- ###### 
  
  ## a. Estimation: A_bin_decideurs_important ####
  set.seed(seed) #Need to set seed for bootstrapping
  
  hyp3_sumresults <- 
    lm(A_bin_decideurs_important ~ treatmentcommune + factor(situation) + factor(region) + #Region and security situation
         factor(respondent_type), #Controlling for respondent type
       data = DM_final, #CBO level dataset
       weights = DM_weight,
       na.action(na.omit)) %>%
    coeftest(., 
             vcov = vcovBS(., 
                           cluster = factor(DM_final$commune),
                           R = 250)) %>% cbo_regression3(outcome.var = "A_bin_decideurs_important") %>%
    mutate(Description = "...most important actors to keep happy in order to maintain influence") %>%
    mutate(control.group.mean = mean((DM_final %>%
                                        drop_na(A_bin_decideurs_important) %>%
                                        filter(treatmentcommune == 0))$A_bin_decideurs_important, na.rm = T))
  
  ## b. Estimation: A_bin_decideurs_justify ####
  set.seed(seed) #Need to set seed for bootstrapping
  
  hyp3_sumresults <- bind_rows(hyp3_sumresults, 
                               lm(A_bin_decideurs_justify ~ treatmentcommune + factor(randomization_block) + #Region and security situation
                                    factor(respondent_type), #Controlling for respondent type
                                  data = DM_final, #CBO level dataset
                                  weights = DM_weight,
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovBS(., 
                                                        cluster = factor(DM_final$commune),
                                                        R = 250)) %>% cbo_regression3(outcome.var = "A_bin_decideurs_justify") %>%
                                 mutate(Description = "... actors to whom municipal decision makers spend the most effort justifying their decisions") %>%
                                 mutate(control.group.mean = mean((DM_final %>%
                                                                     drop_na(A_bin_decideurs_justify) %>%
                                                                     filter(treatmentcommune == 0))$A_bin_decideurs_justify, na.rm = T)),                          
                               
                               ## c. Estimation: pressure_civilsociety ####
                               lm(as.numeric(as.character(pressure_civilsociety)) ~ treatmentcommune + factor(randomization_block) + #Region and security situation
                                    factor(respondent_type), #Controlling for respondent type
                                  data = DM_final, #CBO level dataset
                                  weights = DM_weight,
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovHC(., 
                                                        type = "HC2")) %>% cbo_regression3(outcome.var = "pressure_civilsociety") %>%
                                 mutate(Description = "Perceived likelihood of civil society pressure if their municipality fails to deliver on its
  responsibilities without good justification") %>%
                                 mutate(control.group.mean = mean(as.numeric(as.character((DM_final %>%
                                                                                             drop_na(pressure_civilsociety) %>%
                                                                                             filter(treatmentcommune == 0))$pressure_civilsociety)), na.rm = T))                             
  )
  
  #Adding False Discovery Rate p-value adjustment
  hyp3_sumresults$adj.p.value <- stats::p.adjust( 
    hyp3_sumresults$p.value, method = "fdr", nrow(hyp3_sumresults)) #Creating a vector of adjusted p-values
  
  #Creating two sub-groups to distinguish aggregated variables and its components
  hyp3_sumresults %<>% mutate(sub = ifelse(row_number() < 3, "Community groups and civil society organizations are among the three...", 
                                           ifelse(row_number() == 3, "", ""))) 
  
  ##Creating forest plot: creating a fake metanalysis type table in order to be able to make forest plot
  png(file = paste0(path, 'H3_', 'summary.png'), res = 200, width = 17, height = 6.5, units = "in") #Preparing the export of the graph
  
  meta::metagen( #Creating meta type object
    TE = hyp3_sumresults$estimate,
    lower = hyp3_sumresults$conf.low, #Lower bound of CI
    upper = hyp3_sumresults$conf.high,
    data = hyp3_sumresults,
    method.ci = "t", #t-test for confidence interval
    studlab = hyp3_sumresults$Description, #Vector created above
    subgroup = hyp3_sumresults$sub, 
    level = 0.95
  ) %>% 
    meta::forest(overall = F, #Creating plot
                 rightcols = c("p.value", "adj.p.value", "observations"),
                 rightlabs = c("P-value", "Adj. P-value", "Observations"),
                 leftcols = c("studlab", "control.group.mean", "estimate"),
                 leftlabs = c("Variable", "Control Group \n Mean", "ATE \n (change in proportion/probablity)"), 
                 subgroup = T,
                 subgroup.hetstat = F,
                 xlab = "ATE (change in proportion/probablity) & \n 95 percent confidence interval",
                 fixed = F, #Removing indicator for FE
                 random = F, #Removing indicator for RE
                 subgroup.name = NULL,
                 fontfamily = gpar(fontfamily = "serif"),
                 fontsize = 11,
                 print.I2 = F, 
                 print.tau2 = F, 
                 print.pval.Q = F, 
                 addrows.below.overall = 1)
  
  dev.off()


