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

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

# A_Breadth_of_engagement --------------- ###### 
  
  ### List of index components
  toestimate <- c((CBO_groups_final %>% select(contains("bin_")) %>% colnames))
  
  ## a. Estimation with main outcome var first ####
  
  hyp1_breadth <- rbind(
    lm(A_Breadth_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_Breadth_of_engagement"), #User-written
    
    ## b. Estimations with index components ####
    
    map_dfr(toestimate, function(var){
      set.seed(seed) 
      lm(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. Creating variable description ####
  
  hyp1_breadth$Description <- as.vector(c("Breadth of CBO engagement in muni. governance", c(" Has met: Muni. Council", paste0("  ''             ",   #Listing small description for variables in forest plot
                                                                                                                               c("Cadre de Concer.",
                                                                                                                                 "Secretary General",
                                                                                                                                 "Muni. Admin", 
                                                                                                                                 "Mayor", 
                                                                                                                                 "Councilor", 
                                                                                                                                 "Prefecture", 
                                                                                                                                 "Health auth.", 
                                                                                                                                 "Education auth.", 
                                                                                                                                 "Water auth.", 
                                                                                                                                 "Village chiefs / CVD"))), 
                                          c(" Has collected: data on muni. services", 
                                            paste0("  ''                     ", "opinions")), 
                                          c(" Has interviewed service providers", 
                                            " Has presented at village meetings", 
                                            " Has collaborated with oth. CBOs or equiv.", 
                                            " Has organized: multi-stakeholder meet.", 
                                            " ''                        meet. with CBO or equiv.",
                                            " Has created Whatsapp group or equiv.",
                                            " Has organized information campaign",
                                            " Has reached out to media outlets")
  ))
  
  
  ## d. Adding False Discovery Rate p-value adjustment for all index components ####
  
  hyp1_breadth$adj.p.value <- stats::p.adjust( #Only doing p-value adjustments for index components
    (hyp1_breadth %>% 
       mutate(p.value = ifelse(outcome == "A_Breadth_of_engagement", NA, p.value)))$p.value, 
    method = "fdr", nrow(hyp1_breadth)) #Creating a vector of adjusted p-values
  
  ## e. Creating two sub-groups to distinguish aggregated variables and its components ####
  
  hyp1_breadth %<>% mutate(sub = ifelse(outcome == "A_Breadth_of_engagement", "Main indicator", "Components"))  
  
  ## f. Calculating control group means ####
  covariates <- c("")
  hyp1_breadth$control.group.mean <- 
    sapply(c("A_Breadth_of_engagement", toestimate), function(var){
      round(eval(parse(text = paste0(
        "mean(
            (CBO_groups_final %>% 
              filter(treatment == 0) %>%
              drop_na(",  var, ",", covariates,"))$", var, "
              )"))), digits = 2)
    })
  
  ## g. Creating forest plots: creating a fake metanalysis type table in order to be able to make forest plot ####
  png(file = paste0(path, 'H1_', 'A_Breadth_of_engagement2.png'), res = 200, width = 12, height = 6.5, units = "in") #Preparing the export of the graph
  
  meta::metagen( #Creating meta type object
    TE = (hyp1_breadth %>% filter(outcome != "A_Breadth_of_engagement"))$estimate,
    seTE = (hyp1_breadth %>% filter(outcome != "A_Breadth_of_engagement"))$std.error,
    lower = (hyp1_breadth %>% filter(outcome != "A_Breadth_of_engagement"))$conf.low, #Lower bound of CI
    upper = (hyp1_breadth %>% filter(outcome != "A_Breadth_of_engagement"))$conf.high,
    data = (hyp1_breadth %>% filter(outcome != "A_Breadth_of_engagement")),
    method.ci = "t", #t-test for confidence interval
    studlab = (hyp1_breadth %>% filter(outcome != "A_Breadth_of_engagement"))$Description, #Vector created above
    subgroup = (hyp1_breadth %>% filter(outcome != "A_Breadth_of_engagement"))$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", "Change in Proportion"), 
                 xlab = "Change in proportion & \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()

