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

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

# Hypothesis 2 summary --------------- ###### 

  hyp2_sumresults <- bind_rows(
    
    ## a. Estimation: Public goods game - proportion of endowment contributed ####
    # Tobit regression with cluster-robust standard errors 
    # at the CBO level
    AER::tobit(A_game_group_contrib ~ treatment + factor(randomization_block), #Region and security situation
               data = CBO_groups_final_IND, #Individual level dataset derived from CBO level dataset
               weights = IND_weight, # Weights to compensate for th
               # unequal number of interviewee in each CBO
               na.action(na.omit),
               dist = "gaussian", 
               left = 0,
               right = 1) %>%
      coeftest(., 
               vcov = vcovCL(., 
                             #Clustering at CBO level
                             cluster = factor((CBO_groups_final_IND %>% 
                                                 drop_na(A_game_group_contrib))$appcode), #CBO id
                             type = "HC0")) %>% cbo_regression3(outcome.var = "A_game_group_contrib") %>%
      mutate(Description = "Public goods game: proportion of endowment contributed") %>%
      mutate(control.group.mean = mean((CBO_groups_final_IND %>%
                                          drop_na(A_game_group_contrib) %>%
                                          filter(treatment == 0))$A_game_group_contrib, na.rm = T)), 
    
    ## b. Estimation: the number of general membership meetings: ####
    lm(meetings_general ~ treatment + factor(randomization_block), #Region and security situation
       data = CBO_groups_final, #CBO level dataset
       na.action(na.omit)) %>%
      coeftest(., 
               vcov = vcovHC(., 
                             type = "HC2")) %>% cbo_regression3(outcome.var = "meetings_general") %>% 
      mutate(Description = "General membership meetings held by the CBO in 2020")  %>%
      mutate(control.group.mean = mean((CBO_groups_final %>%
                                          drop_na(meetings_general) %>%
                                          filter(treatment == 0))$meetings_general, na.rm = T)),
    
    ## c. Estimation: the number of executive meetings: ####
    lm(meetings_executive ~ treatment + factor(randomization_block), #Region and security situation
       data = CBO_groups_final, #CBO level dataset
       na.action(na.omit)) %>%
      coeftest(., 
               vcov = vcovHC(., 
                             type = "HC2")) %>% cbo_regression3(outcome.var = "meetings_executive") %>%
      mutate(Description = "Executive board meetings held by the CBO in 2020") %>%
      mutate(control.group.mean = mean((CBO_groups_final %>%
                                          drop_na(meetings_executive) %>%
                                          filter(treatment == 0))$meetings_executive, na.rm = T)))
  
  ## d. Estimation: cbo had a budget in 2020 ####
  set.seed(seed) #Need to set seed for bootstrapping
  
  hyp2_sumresults <- bind_rows(hyp2_sumresults,
                               
                               lm(as.numeric(as.character(budget_any)) ~ treatment + factor(randomization_block), #Region and security situation
                                  data = CBO_groups_final, #CBO level dataset
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovBS(., 
                                                        R = 250)) %>% cbo_regression3(outcome.var = "budget_any") %>%
                                 mutate(Description = "CBO had a formal budget in 2020") %>%
                                 mutate(control.group.mean = mean(as.numeric(as.character((CBO_groups_final %>%
                                                                                             drop_na(budget_any) %>%
                                                                                             filter(treatment == 0))$budget_any)), na.rm = T)),
                               
                               ## e. Estimation: net membership change ####
                               lm(A_log_members_total_after ~ treatment + factor(randomization_block) + #Region and security situation
                                    A_log_members_total_before, 
                                  data = CBO_groups_final, #CBO level dataset
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovHC(., 
                                                        type = "HC2")) %>% cbo_regression3(outcome.var = "A_log_members_total_after") %>%
                                 mutate(Description = "Net membership change") %>%
                                 mutate(estimate = exp(estimate) - 1) %>% #De-loging estimates
                                 mutate(std.error = exp(std.error) - 1) %>% #If need de-logged
                                 mutate(conf.low = exp(conf.low) - 1) %>% #If need de-logged
                                 mutate(conf.high = exp(conf.high) - 1) %>% #If need de-logged
                                 mutate(control.group.mean = exp(mean((CBO_groups_final %>%
                                                                         drop_na(A_log_members_total_after, A_log_members_total_before) %>%
                                                                         filter(treatment == 0))$A_log_members_total_after, na.rm = T))))
  
  ## f. Estimation: CBO has adopted COVID-19 policy ####
  set.seed(seed) #Need to set seed for bootstrapping
  
  hyp2_sumresults <- bind_rows(hyp2_sumresults,
                               
                               lm(as.numeric(as.character(covid_policy)) ~ treatment + factor(randomization_block), #Region and security situation
                                  data = CBO_groups_final, #CBO level dataset
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovBS(., 
                                                        R = 250)) %>% cbo_regression3(outcome.var = "covid_policy") %>%
                                 mutate(Description = "CBO has adopted COVID-19 policy") %>%
                                 mutate(control.group.mean = mean(as.numeric(as.character((CBO_groups_final %>%
                                                                                             drop_na(covid_policy) %>%
                                                                                             filter(treatment == 0))$covid_policy)), na.rm = T)),
                               
                               ## g. Estimation: CBO has seized COVID-related opportunities ####
                               lm(A_covid_adaptation ~ treatment + factor(randomization_block), #Region and security situation
                                  data = CBO_groups_final, #CBO level dataset
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovHC(., 
                                                        type = "HC2")) %>% cbo_regression3(outcome.var = "A_covid_adaptation") %>% 
                                 mutate(Description = "CBO has seized COVID-related opportunities") %>%
                                 mutate(control.group.mean = mean((CBO_groups_final %>%
                                                                     drop_na(A_covid_adaptation) %>%
                                                                     filter(treatment == 0))$A_covid_adaptation, na.rm = T)))
  
  
  ## h. Estimation: Perceived internal conflict ####
  set.seed(seed) #Need to set seed for bootstrapping
  
  hyp2_sumresults <- bind_rows(hyp2_sumresults,
                               
                               lm(A_cbo_conflict ~ treatment + factor(randomization_block), #Region and security situation
                                  data = CBO_individuals_final, #CBO level dataset
                                  weights = IND_weight,
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovBS(., 
                                                        cluster = factor(CBO_individuals_final$appcode),
                                                        R = 250)) %>% cbo_regression3(outcome.var = "A_cbo_conflict") %>%
                                 mutate(Description = "Perceived internal conflict") %>%
                                 mutate(control.group.mean = mean((CBO_individuals_final %>%
                                                                     drop_na(A_cbo_conflict) %>%
                                                                     filter(treatment == 0))$A_cbo_conflict, na.rm = T)), #Control group mean
                               
                               ## i. Estimation: Satisfaction with the CBO leadership ####
                               lm(as.numeric(as.character(cbo_satisfaction)) ~ treatment + factor(randomization_block), #Region and security situation
                                  data = CBO_individuals_final, #CBO level dataset
                                  weights = IND_weight,
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovCL(., 
                                                        cluster = factor(CBO_individuals_final$appcode),
                                                        type = "HC2")) %>% cbo_regression3(outcome.var = "cbo_satisfaction") %>%
                                 mutate(Description = "Satisfaction with the CBO leadership") %>%
                                 mutate(control.group.mean = mean(as.numeric(as.character((CBO_individuals_final %>%
                                                                                             drop_na(cbo_satisfaction) %>%
                                                                                             filter(treatment == 0))$cbo_satisfaction)), na.rm = T)),
                               
                               ## j. Estimation: Enthusiasm about being a member ####
                               lm(as.numeric(as.character(cbo_enthusiasm)) ~ treatment + factor(randomization_block), #Region and security situation
                                  data = CBO_individuals_final, #CBO level dataset
                                  weights = IND_weight,
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovCL(., 
                                                        cluster = factor(CBO_individuals_final$appcode),
                                                        type = "HC2")) %>% cbo_regression3(outcome.var = "cbo_enthusiasm") %>%
                                 mutate(Description = "Enthusiasm about being a member") %>%
                                 mutate(control.group.mean = mean(as.numeric(as.character((CBO_individuals_final %>%
                                                                                             drop_na(cbo_enthusiasm) %>%
                                                                                             filter(treatment == 0))$cbo_enthusiasm)), na.rm = T)))
  
  
  ## k. Estimation: Perceived capture of CBO’s activities ####
  set.seed(seed) #Need to set seed for bootstrapping
  
  hyp2_sumresults <- bind_rows(hyp2_sumresults,
                               
                               lm(as.numeric(as.character(cbo_capture)) ~ treatment + factor(randomization_block), #Region and security situation
                                  data = CBO_individuals_final, #CBO level dataset
                                  weights = IND_weight,
                                  na.action(na.omit)) %>%
                                 coeftest(., 
                                          vcov = vcovBS(., 
                                                        cluster = factor(CBO_individuals_final$appcode),
                                                        R = 250)) %>% cbo_regression3(outcome.var = "cbo_capture") %>%
                                 mutate(Description = "Perceived capture of CBO’s activities") %>%
                                 mutate(control.group.mean = mean(as.numeric(as.character((CBO_individuals_final %>%
                                                                                             drop_na(cbo_capture) %>%
                                                                                             filter(treatment == 0))$cbo_capture)), na.rm = T)))
  
  ## l. Adding False Discovery Rate p-value adjustment ####
  hyp2_sumresults$adj.p.value <- stats::p.adjust( 
    hyp2_sumresults$p.value, method = "fdr", nrow(hyp2_sumresults)) #Creating a vector of adjusted p-values
  
  ## m. Creating sub-groups ####
  hyp2_sumresults %<>% mutate(sub = ifelse(row_number() < 2, "Collective Action Capacity", 
                                           ifelse(row_number() >= 2 & row_number() <= 7, "Organizational Capacity", 
                                                  ifelse(row_number() >= 8 & row_number() <= 10, "Internal cohesion", 
                                                         ifelse(row_number() == 11, "Organizational capture", "")))))
  
  ## n. Hyp2 forest plot ####
  png(file = paste0(path, 'H2_', 'summary.png'), res = 200, width = 12, height = 6.5, units = "in") #Preparing the export of the graph
  
  meta::metagen( #Creating meta type object
    TE = hyp2_sumresults$estimate,
    seTE = hyp2_sumresults$std.error,
    lower = hyp2_sumresults$conf.low, #Lower bound of CI
    upper = hyp2_sumresults$conf.high,
    data = hyp2_sumresults,
    method.ci = "t", #t-test for confidence interval
    studlab = hyp2_sumresults$Description, #Vector created above
    subgroup = hyp2_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", "Average Treatment \n Effect"), 
                 xlab = "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()