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

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

# Hypothesis 1 summary --------------- ###### 

  ## a. Estimations with standardized outcome variables ####
  
  hyp1summary <- bind_rows( 
    
    lm(scale(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
      mutate(Description = "Breadth of CBO engagement in muni. governance"),
    
    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
      mutate(Description = "Intensity of CBO engagement in muni. governance"),
    
    lm(scale(A_DM_know_and_inter) ~ 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_know_and_inter") %>% #User-written
      mutate(Description = "Decision makers' knowledge of and interaction with CBOs")
  )
  
  ## b. Adding False Discovery Rate p-value adjustment for all main outcome vars ####
  hyp1summary$adj.p.value <- stats::p.adjust(hyp1summary$p.value, 
                                             method = "fdr", nrow(hyp1summary)) #Creating a vector of adjusted p-values
  
  ## c. Calculating control group means ####
  hyp1summary$control.group.mean <- 
    c(
      round(
        mean((CBO_groups_final %>% 
                drop_na(A_Breadth_of_engagement) %>%
                mutate(A_Breadth_of_engagement = scale(A_Breadth_of_engagement)) %>%
                filter(treatment == 0))$A_Breadth_of_engagement), 
        digits = 2), 
      round(
        mean((CBO_groups_final %>% 
                drop_na(A_Intensity_of_engagement) %>%
                mutate(A_Intensity_of_engagement = scale(A_Intensity_of_engagement)) %>%
                filter(treatment == 0))$A_Intensity_of_engagement), 
        digits = 2), 
      round(
        mean((DM_finalx2 %>% 
                drop_na(A_DM_know_and_inter, respondent_type) %>%
                mutate(A_DM_know_and_inter = scale(A_DM_know_and_inter)) %>%
                filter(treatmentcommune == 0))$A_DM_know_and_inter, na.rm = T), 
        digits = 2)
    )
  
  ### Creating forest plot: creating a fake metanalysis type table in order to be able to make forest plot
  png(file = paste0(path, 'H1_', 'summary.png'), res = 200, width = 12, height = 6.5, units = "in") #Preparing the export of the graph
  
  meta::metagen( #Creating meta type object
    TE = hyp1summary$estimate,
    seTE = hyp1summary$std.error,
    lower = hyp1summary$conf.low, #Lower bound of CI
    upper = hyp1summary$conf.high,
    data = hyp1summary,
    method.ci = "t", #t-test for confidence interval
    studlab = hyp1summary$Description, #Vector created above
    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", "Average Treatment Effect \n (standardized)"), 
                 xlab = "Standardized ATE & \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()
