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

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

# A_Intensity_of_engagement --------------- ######

  ### List of index components
  toestimate <- c("meet_mayor_num", "meet_councilor_num", "meet_cm", "meet_cdc", "meet_sg_num", "meet_admin_num", "meet_health_total", "meet_education_total", "meet_water_total",  "organize_stakeholder_num", "organize_collab_num", "organize_villagemeeting_num", "pressure_campaign_num")
  
  ## a. Estimation with main outcome var first ####
  
  hyp1_intensity <- rbind(
    lm(A_Intensity_of_engagement ~ treatment + factor(randomization_block), #Region and security situation
       data = CBO_groups_final,
       na.action(na.omit)) %>%
      coeftest(., 
               vcov = vcovHC(., 
                             type = "HC2")) %>% cbo_regression3(outcome.var = "A_Intensity_of_engagement"), #User-written
    
    ## b. Estimations with index components ####
    
    map_dfr(toestimate, function(var){
      set.seed(seed) 
      #Standardized outcome variable
      lm(scale(eval(parse(text = paste0("CBO_groups_final$", var)))) ~ treatment + factor(randomization_block), 
         data = CBO_groups_final,
         na.action(na.omit)) %>% 
        coeftest(., 
                 vcov = vcovBS(., 
                               R = 250)) %>% cbo_regression3(outcome.var = var)}),
    
    ## c. Number of yearly DM/CBO interactions ####
    
    lm(scale(A_DM_freqinter) ~ treatmentcommune + factor(randomization_block) + #Region and security situation 
         factor(respondent_type), #Controlling for respondent type
       data = DM_finalx2,
       weights = DM_weight,
       na.action(na.omit)) %>%
      coeftest(., 
               vcov = vcovCL(., 
                             #Clustering at CBO level
                             cluster = factor(DM_finalx2$appcode), #CBO id 
                             type = "HC2")) %>% cbo_regression3(outcome.var = "A_DM_freqinter")
  )
  
  ## d. Need to create vector with small variable description for plot ####
  hyp1_intensity$Description <- as.vector(c("Average: intensity of engagement", c("Nb. of meetings with: elected - Mayor", 
                                                                                  "''                                    Councilor", 
                                                                                  "''                                    Muni. Council", 
                                                                                  "''                                    Cadre de Concert.", 
                                                                                  "''                          unelected - Secretary General",
                                                                                  "''                          unelected - Muni. Admin",
                                                                                  "''                          service providers - health",
                                                                                  "''                          service providers - education",
                                                                                  "''                          service providers - water",
                                                                                  "Nb. of organized multi-stake. meet.",
                                                                                  "Nb. of CBOs or equiv. collaborated with", 
                                                                                  "Nb. of villages in which CBO presented", 
                                                                                  "Nb. of organized information campaign"), 
                                            "Nb. of yearly DM/CBO interactions"))
  
  ## e. Adding False Discovery Rate p-value adjustment for all index components ####
  hyp1_intensity$adj.p.value <- stats::p.adjust( #Only doing p-value adjustments for index components
    (hyp1_intensity %>% 
       mutate(p.value = ifelse(outcome == "A_Intensity_of_engagement" | outcome == "A_DM_freqinter", NA, p.value)))$p.value, 
    method = "fdr", nrow(hyp1_intensity)) #Creating a vector of adjusted p-values
  
  ## f. Creating three sub-groups to distinguish aggregated variables and its components, and A_DM_freqinter ####
  hyp1_intensity %<>% mutate(sub = ifelse(outcome != "A_DM_freqinter", "According to CBOs", "According to Decision Makers"))
  
  ## g. Calculating control group means ####
  covariates <- c("")
  hyp1_intensity$control.group.mean <- 
    c(
      round(
        mean((CBO_groups_final %>% 
                filter(treatment == 0) %>%
                drop_na(A_Intensity_of_engagement))$A_Intensity_of_engagement), 
        digits = 2),
      
      sapply(c(toestimate), function(var){
        round(eval(parse(text = paste0(
          "mean(
            (CBO_groups_final %>% 
              drop_na(",  var, ",", covariates,") %>%
              mutate(", var, "= scale(", var, ")) %>%
              filter(treatment == 0))$", var, ")"))), 
          digits = 2)
      }), 
      
      round(
        mean((DM_finalx2 %>% 
                drop_na(A_DM_freqinter, respondent_type) %>%
                mutate(A_DM_freqinter = scale(A_DM_freqinter)) %>%
                filter(treatmentcommune == 0))$A_DM_freqinter), 
        digits = 2)
    )
  
  ## h. Creating forest plot: creating a fake metanalysis type table in order to be able to make forest plot ####
  
  ###Preparing the export of the graph
  png(file = paste0(path, 'H1_', 'A_Intensity_of_engagement.png'), res = 200, width = 12, height = 6.5, units = "in") 
  
  meta::metagen( #Creating meta type object
    TE = hyp1_intensity$estimate,
    seTE = hyp1_intensity$std.error,
    lower = hyp1_intensity$conf.low, #Lower bound of CI
    upper = hyp1_intensity$conf.high,
    data = hyp1_intensity,
    method.ci = "t", #t-test for confidence interval
    studlab = hyp1_intensity$Description, #Vector created above
    subgroup = hyp1_intensity$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", "Standardized Average \n Treatment Effect"), 
                 xlab = "Standardized ATE & \n 95 percent confidence interval",
                 subgroup = T,
                 subgroup.hetstat = F,
                 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()

