
#############################################################################################################
############################################# MAIN EFFECTS ##################################################
#############################################################################################################


rm(list = ls()) ## Start from clean workspace

require(Amelia)

load("data/ghana_citizen_imputed_manip.RData")
load("data/uganda_citizen_imputed_manip.RData")
load("data/ghana_citizen_imputed_full.RData")
load("data/uganda_citizen_imputed_full.RData")

numImp <- 5 # This will be used in the imputation loop

## Source a function to make equations (to remove clutter)

make.equation <- function(dvs, ivs, interactions = c()){
  return(paste(dvs, "~", paste(c(ivs, interactions), sep = "", collapse = " + ")))
}

##################### LOOP OVER IMPUTATIONS ####################################

imp.results <- c() # Empty container to store the results
for(m in 0:numImp){ # Start this at 0 if you want to also do analysis over unimputed data
  
  ## Set up data 
  
  data.list <- list(ug.imp.full[ug.imp.full$imp == m, ], 
                    ug.imp.manip[ug.imp.manip$imp == m, ], 
                    gh.imp.full[gh.imp.full$imp == m, ], 
                    gh.imp.manip[gh.imp.manip$imp == m, ])
  manip.list <- c("no", "yes", "no", "yes") # Note, these two naming vectors should properly match the data list in terms of contents
  country.list <- c("uganda", "uganda", "ghana", "ghana")
  var.info <- read.csv("data/mass_analysis_varinfo.csv", stringsAsFactors = F)
  
  ######################## EQUATION PREP ##########################################
  
  treat.vars <- list(c("oil", "aid", "ngo"), # Tax as base
                     c("aid", "ngo", "tax"), # Oil as base
                     c("ngo", "tax", "oil"), # Aid as base
                     c("tax", "oil", "aid"), # NGO as base
                     c("ngo")) # Non-NGO as base
  ref.group <- c("tax", "oil", "aid", "ngo", "non-ngo")
  vars.fe <- c("as.factor(enum_id)")
  
  ## Note: what the treat.vars list object is doing is setting up a set of dummies
  ## to be included in the estimating equation such that the omitted variable
  ## is the base category. The ``ref group'' object is the omitted feature
  ## and is used below to specify the reference group in the output table
  
  ####################### BEGIN LOOP ###############################################
  
  table.results <- c() ## Container for results
  
  for(l in 1:length(treat.vars)){ # For each set of treatment variables (e.g. each element in list)
    
    ## GENERATE EQUATIONS 
    
    treat.use <- treat.vars[[l]] # Grab the treatment variables to include in equation
    treat.ref <- ref.group[l] # Set the reference group
    eqs.fe <- make.equation(dvs = var.info$var_name, ivs = c(treat.use, vars.fe)) # Add FEs (defined above)
    eqs.type <- c("fe") # Name types models to run, will loop over this below
    eqs.list <- list(eqs.fe) # Gather corresponding equations into lists
    
    ## NOTE: CAN REMOVE CERTAIN MODELS BY COMMENTING OUT AND ADJUST eqs.type and eqs.list objects
    
    
    ############################ ESTIMATION ######################
    
    model.obj <- list()
    for(k in 1:length(eqs.type)){
      for(j in 1:length(data.list)){ # This list is defined above, includes full UG + GH datasets and then also with subjects who failed the manipulation check removed
        for(i in 1:length(var.info$var_name)){
          
          tab <- c() # Container for the i loop
          mod <- glm(eqs.list[[k]][i], data = data.list[[j]])
          if(length(treat.use) == 1){ # if it's just a single dummy, have to transpose
            tab <- data.frame(t(as.data.frame(summary(mod)$coefficients[c(treat.use), ])))
            rownames(tab) <- treat.use}
          else{ # if multiple treatment dummies, can take summary frame as-is
            tab <- data.frame(summary(mod)$coefficients[c(treat.use), ])}
        
          ## Create summary object from estimates
          
          colnames(tab) <- c("est", "se", "t", "p")
          tab$iv <- rownames(tab)
          tab$model <- eqs.type[k]
          tab$country <- country.list[j]
          tab$manip <- manip.list[j]
          tab$var_group <- var.info$var_index[i]
          tab$type <- var.info$var_type[i]
          tab$sig <- ifelse(tab$p >= 0.10, "no", "yes")
          tab$dv <- var.info$var_name[i]
          tab$label <- var.info$var_label[i]
          tab$reference_category <- treat.ref
          tab$imputation <- m
          table.results <- rbind(tab, table.results)
        } # Close i loop (over variables)
      } # Close j loop (over datasets)
    } # Close model loop (over different equations)
  } # Close l loop (over different base categories)
  
  useCols <- c("label", "iv", "est", "se", "t", "p", "sig",
               "country", "manip", "var_group", 
               "reference_category", "dv", "model", "imputation")
  lm.results <- table.results[, useCols]
  lm.results <- lm.results[with(lm.results, order(country, reference_category, iv, dv, 
                                                  var_group, model, manip, imputation)), ]
  lm.results$result_id <- 1:nrow(lm.results)
  rownames(lm.results) <- NULL
  
  ## Store the results in a list, where each element is an imputation
  
  imp.results <- rbind(imp.results, lm.results)
  
  
} # Close m loop over imputations

## Subset to the imputed data only and combine them

res.id <- unique(imp.results$result_id)
mi.results <- c()
for(i in 1:length(res.id)){
  df.tmp <- imp.results[imp.results$result_id == i & imp.results$imputation != 0, ]
  est.mi <- unlist(mi.meld(q = cbind(df.tmp[, "est"], df.tmp[, "est"]), 
                           se = cbind(df.tmp[, "se"], df.tmp[, "se"])))[c(1, 3)] # element 1-2 is beta, 3-4 is se (because of repetition)
  # Note, need to duplicate here because mi.meld function won't handle vectors 
  all.mi <- c(est.mi[1], est.mi[2], est.mi[1]/est.mi[2], 2*(1 - pnorm(abs(est.mi[1]/est.mi[2])))) 
  mi.results <- rbind(mi.results, all.mi)}
final.results <- imp.results[imp.results$imp == 0, ]
final.results[, paste("mi_", c("est", "se", "t", "p"), sep = "")] <- mi.results
final.results$mi_sig <- ifelse(final.results$p <= 0.05, 1, 0)

# Round the numeric columns
numeric.vars <- c("est", "se", "t", "p", 
                  paste("mi_", c("est", "se", "t", "p"), sep = ""))
final.results[, numeric.vars] <- round(final.results[, numeric.vars], digits = 5)

# Save object
save(final.results, file = "results/mass_maineffects_results_appendix.RData")


#############################################################################################################
############################################# SUBGROUP EFFECTS ##############################################
#############################################################################################################


rm(list = ls()) ## Start from clean workspace


load("data/ghana_citizen_imputed_manip.RData")
load("data/uganda_citizen_imputed_manip.RData")
numImp <- 5 # This will be used in the imputation loop

## Source a function to make equations (to remove clutter)

make.equation <- function(dvs, ivs, interactions = c()){
  return(paste(dvs, "~", paste(c(ivs, interactions), sep = "", collapse = " + ")))
}

##################### LOOP OVER IMPUTATIONS ####################################

imp.results <- c() # Empty container to store the results
for(m in 0:numImp){ # Start this at 0 if you want to also do analysis over unimputed data
  
  ## Set up data 
  
  data.list <- list(ug.imp.manip[ug.imp.manip$imp == m, ], 
                    gh.imp.manip[gh.imp.manip$imp == m, ])
  manip.list <- c("yes", "yes") # Note, these two naming vectors should properly match the data list in terms of contents
  country.list <- c("uganda", "ghana")
  var.info <- read.csv("data/mass_analysis_varinfo.csv", stringsAsFactors = F)
   
  
  ######################## EQUATION PREP ##########################################
  
  treat.vars <- list(c("oil", "aid", "ngo"), # Tax as base
                     c("aid", "ngo", "tax"), # Oil as base
                     c("ngo", "tax", "oil"), # Aid as base
                     c("tax", "oil", "aid"), # NGO as base
                     c("ngo")) # Non-NGO as base 
  ref.group <- c("tax", "oil", "aid", "ngo", "non-ngo")
  interaction.vars <- c("type_high", "gov_supp", "oil_region", "female", "trust_low", 
                        "noncoethnic", "perf_low", "corr_high", "poverty_high",
                        "poverty_govsupp", "ngo_hightrust")
  vars.fe <- c("as.factor(enum_id)")
  
  ## Note: what the treat.vars list object is doing is setting up a set of dummies
  ## to be included in the estimating equation such that the omitted variable
  ## is the base category. The ``ref group'' object is the omitted feature
  ## and is used below to specify the reference group in the output table
  
  ####################### BEGIN LOOP ###############################################
  
  table.results <- c() ## Container for results
  
  for(l in 1:length(treat.vars)){ # For each set of treatment variables (e.g. each element in list)
    for(u in 1:length(interaction.vars)){ # For each interaction term 
      ## GENERATE EQUATIONS 
      
      treat.use <- treat.vars[[l]] # Grab the treatment variables to include in equation
      treat.ref <- ref.group[l] # Set the reference group
      interact.terms <-  unlist(lapply(treat.use, FUN = function(c){
        paste(c, ":", interaction.vars[u], sep = "")})) # Generate interaction terms, leave base term out (b/c of randomization)
      eqs.interact <- make.equation(dvs = var.info$var_name, ivs = c(treat.use, interaction.vars[u], vars.fe), 
                                    interactions = interact.terms) # Add interaction term (keep FE)
      eqs.type <- c("interaction") # Name types models to run, will loop over this below
      eqs.list <- list(eqs.interact) # Gather corresponding equations into lists
      
      ## NOTE: CAN REMOVE CERTAIN MODELS BY COMMENTING OUT AND ADJUST eqs.type and eqs.list objects
      
      
      ############################ ESTIMATION ######################
      
      model.obj <- list()
      for(k in 1:length(eqs.type)){
        for(j in 1:length(data.list)){ # This list is defined above, includes full UG + GH datasets and then also with subjects who failed the manipulation check removed
          for(i in 1:length(var.info$var_name)){
            
            ## Run models and store objects depending on variable type
            
            tab <- c() # Build container to hold model estimates and identifying information
            mod <- glm(eqs.list[[k]][i], data = data.list[[j]])
            model.obj[[i]] <- mod
            interact.rows <- rownames(summary(mod)$coefficients)[grep(":", rownames(summary(mod)$coefficients))]
            tab <- data.frame(summary(mod)$coefficients[c(treat.use, interact.rows), ])
            
            
            ## Create summary object from estimates
            colnames(tab) <- c("est", "se", "t", "p")
            tab$iv <- rownames(tab)
            tab$model <- eqs.type[k]
            tab$country <- country.list[j]
            tab$manip <- manip.list[j]
            tab$var_group <- var.info$var_index[i]
            tab$type <- var.info$var_type[i]
            tab$sig <- ifelse(tab$p >= 0.10, "no", "yes")
            tab$dv <- var.info$var_name[i]
            tab$label <- var.info$var_label[i]
            tab$reference_category <- treat.ref
            tab$subgroup <- interaction.vars[u]
            tab$imputation <- m
            table.results <- rbind(tab, table.results)
          } #Close u loop (over subgroups)
        } # Close i loop (over variables)
      } # Close j loop (over datasets)
    } # Close model loop (over different equations)
  } # Close l loop (over different base categories)
  
  useCols <- c("label", "iv", "est", "se", "t", "p", "sig",
               "country", "manip", "var_group", 
               "reference_category", "dv", "model", "imputation")
  
  lm.results <- table.results[, useCols]
  lm.results <- lm.results[with(lm.results, order(country, reference_category, iv, dv, 
                                                  var_group, model, manip, imputation)), ]
  lm.results$result_id <- 1:nrow(lm.results)
  rownames(lm.results) <- NULL
  
  ## Store the results in a list, where each element is an imputation
  
  imp.results <- rbind(imp.results, lm.results)
  
  
} # Close m loop over imputations

## Subset to the imputed data only and combine them

res.id <- unique(imp.results$result_id)
mi.results <- c()
for(i in 1:length(res.id)){
  df.tmp <- imp.results[imp.results$result_id == i & imp.results$imputation != 0, ]
  est.mi <- unlist(mi.meld(q = cbind(df.tmp[, "est"], df.tmp[, "est"]), 
                           se = cbind(df.tmp[, "se"], df.tmp[, "se"])))[c(1, 3)] # element 1-2 is beta, 3-4 is se (beta and ses should be identical)
  # Note, need to duplicate estimate and se here because mi.meld function won't handle vectors, 
  all.mi <- c(est.mi[1], est.mi[2], est.mi[1]/est.mi[2], 2*(1 - pnorm(abs(est.mi[1]/est.mi[2])))) # Note that est.mi is of length 2
  mi.results <- rbind(mi.results, all.mi)} # Close i loop

## Put these results alongside the unimputed results

het.results <- imp.results[imp.results$imp == 0, ]
het.results[, paste("mi_", c("est", "se", "t", "p"), sep = "")] <- mi.results

## Get rid of some decimal points in the numeric columns

numeric.vars <- c("est", "se", "t", "p", 
                  paste("mi_", c("est", "se", "t", "p"), sep = ""))
het.results[, numeric.vars] <- round(het.results[, numeric.vars], digits = 5)

## Save object to produce tables

save(het.results, file = "results/mass_heteffects_results_appendix.RData")



#########################################################################
########################### MP SUBGROUP #################################
#########################################################################

rm(list = ls()) ## Start from clean workspace

## Load MP data
load("data/ghana_mp_full.RData")
load("data/uganda_mp_full.RData")


## Source a function to make equations (to remove clutter)

make.equation <- function(dvs, ivs, interactions = c()){
  return(paste(dvs, "~", paste(c(ivs, interactions), sep = "", collapse = " + ")))
}


country.list <- c("uganda", "ghana")
data.list <- list(d.ug.mp, d.gh.mp)
names(data.list) <- c("uganda", "ghana")

var.info <- read.csv("data/mp_analysis_varinfo.csv", stringsAsFactors = F)


######################## EQUATION PREP ##########################################



treat.vars <- list(c("oil", "aid", "ngo"), # Tax as base
                   c("aid", "ngo", "tax"), # Oil as base
                   c("ngo", "tax", "oil"), # Aid as base
                   c("tax", "oil", "aid"), # NGO as base
                   c("ngo")) # Non-NGO as base 
ref.group <- c("tax", "oil", "aid", "ngo", "non-ngo")
interaction.vars <- c( "mp_ruling", "services_poor", "corr_high")
vars.fe <- c("as.factor(enum_id)")



table.results <- c() ## Container for results
interact.keep <- c() ## This stores the names of the interaction terms, will be used to subset out the unnecessary comparisons before saving the final results in an RData object
for(l in 1:length(treat.vars)){ # For each set of treatment variables (e.g. each element in list)
  for(u in 1:length(interaction.vars)){ # For each interaction term 
    ## GENERATE EQUATIONS 
    
    treat.use <- treat.vars[[l]] # Grab the treatment variables to include in equation
    treat.ref <- ref.group[l] # Set the reference group
    interact.terms <-  unlist(lapply(treat.use, FUN = function(c){
      paste(c, ":", interaction.vars[u], sep = "")})) # Generate interaction terms, leave base term out (b/c of randomization)
    eqs.interact <- make.equation(dvs = var.info$var_name, ivs = c(treat.use, interaction.vars[u], vars.fe), 
                                  interactions = interact.terms) # Add interaction term (keep FE)
    eqs.type <- c("interaction") # Name types models to run, will loop over this below
    eqs.list <- list(eqs.interact) # Gather corresponding equations into lists
    interact.keep <- c(interact.terms, interact.keep)
    
    ## NOTE: CAN REMOVE CERTAIN MODELS BY COMMENTING OUT AND ADJUST eqs.type and eqs.list objects
    
    
    ############################ ESTIMATION ######################
    
    model.obj <- list()
    for(k in 1:length(eqs.type)){
      for(j in 1:length(data.list)){ # This list is defined above, includes full UG + GH datasets and then also with subjects who failed the manipulation check removed
        for(i in 1:length(var.info$var_name)){
          
          ## Run models and store objects depending on variable type
          
          tab <- c()
          if(var.info$var_type[i] %in% c("continuous", "ordered", "binary")){
            mod <- glm(eqs.list[[k]][i], data = data.list[[j]])
            model.obj[[i]] <- mod
            interact.rows <- rownames(summary(mod)$coefficients)[grep(":", rownames(summary(mod)$coefficients))]
            
            if(length(treat.use) == 1 & eqs.type[k] != "interaction"){ # if it's just a single dummy, have to transpose
              tab <- data.frame(t(as.data.frame(summary(mod)$coefficients[c(treat.use, interact.rows), ])))
              rownames(tab) <- treat.use}
            else{ # if multiple treatment dummies, can take summary frame as-is
              tab <- data.frame(summary(mod)$coefficients[c(treat.use, interact.rows), ])}}
          
          
          ## Create summary object from estimates
          colnames(tab) <- c("est", "se", "t", "p")
          tab$iv <- rownames(tab)
          tab$country <- country.list[j]
          tab$var_group <- var.info$var_index[i]
          tab$type <- var.info$var_type[i]
          tab$sig <- ifelse(tab$p >= 0.10, "no", "yes")
          tab$dv <- var.info$var_name[i]
          tab$label <- var.info$var_label[i]
          tab$reference_category <- treat.ref
          tab$subgroup <- interaction.vars[u]
          table.results <- rbind(tab, table.results)
        } #Close u loop (over subgroups)
      } # Close i loop (over variables)
    } # Close j loop (over datasets)
  } # Close model loop (over different equations)
} # Close l loop (over different base categories)

useCols <- c("label", "iv", "est", "se", "t", "p", "sig",
             "country", "var_group", "subgroup",
             "reference_category", "dv")

lm.results <- table.results[, useCols]
lm.results <- lm.results[with(lm.results, order(country, reference_category, iv, dv, 
                                                var_group)), ]
lm.results$result_id <- 1:nrow(lm.results)
rownames(lm.results) <- NULL

numeric.vars <- c("est", "se", "t", "p")
lm.results[, numeric.vars] <- round(lm.results[, numeric.vars], digits = 5)

## Keep only the comparisons of interest

lm.results$comparison <- paste(lm.results$iv, lm.results$reference_category, sep = "-")
lm.results <- lm.results[lm.results$iv %in% interact.keep, ]

save(lm.results, file = "results/mp_heteffects_appendix.RData")


###################################################################################
################################ MASS RI ANALYSIS #################################
###################################################################################

rm(list = ls()) ## Start from clean workspace

set.seed(1234)
# Function to make equations easier
make.equation <- function(dvs, ivs, interactions = c()){
  return(paste(dvs, "~", paste(c(ivs, interactions), sep = "", collapse = " + ")))
}

packages <- c("ri2","randomizr", "dplyr")
lapply(packages, require, character.only = T)

load("data/ghana_citizen_imputed_full.RData")
load("data/uganda_citizen_imputed_full.RData")

## Data and setup

unimp.ug <- ug.imp.full[ug.imp.full$imp == 0, ] ## Uganda - Full Sample, unimputed
unimp.gh <- gh.imp.full[gh.imp.full$imp == 0, ] ## Ghana - Full Sample, unimputed
var.info <- read.csv("data/mass_analysis_varinfo.csv", stringsAsFactors = F)

## Create binary variables that will automatically subset out all except 2 treatment conditions
## Note: 1 = oil, 2 = aid, 3 = tax, 4 = ngo

# Ghana
unimp.gh$oilaid <- dplyr::recode(unimp.gh$treatment, `1` = "1", `2`= "0", `3`="NA", `4`="NA") # Note that function will spit a warning because "NA" will be treated as NAs by R
unimp.gh$oilaid <- as.numeric(unimp.gh$oilaid)
unimp.gh$taxaid <- dplyr::recode(unimp.gh$treatment, `1` = "NA", `2`= "0", `3`="1", `4`="NA")
unimp.gh$taxaid <- as.numeric(unimp.gh$taxaid)
unimp.gh$taxoil <- dplyr::recode(unimp.gh$treatment, `1` = "0", `2` = "NA", `3` = "1", `4` = "NA")
unimp.gh$taxoil <- as.numeric(unimp.gh$taxoil) 

#Uganda
unimp.ug$oilaid <- dplyr::recode(unimp.ug$treatment, `1` = "1", `2`= "0", `3`="NA", `4`="NA")
unimp.ug$oilaid <- as.numeric(unimp.ug$oilaid)
unimp.ug$taxaid <- dplyr::recode(unimp.ug$treatment, `1` = "NA", `2`= "0", `3`="1", `4`="NA")
unimp.ug$taxaid <- as.numeric(unimp.ug$taxaid)
unimp.ug$taxoil <- dplyr::recode(unimp.ug$treatment, `1` = "0", `2` = "NA", `3` = "1", `4` = "NA")
unimp.ug$taxoil <- as.numeric(unimp.ug$taxoil)

## Specify which subsets and treatment variables you want

treat.vars <- list(c("oilaid"), # Note that these variables are constructed such that the first condition is treatment and second is control, all subjects not in either of these 2 treatments will be NA
                   c("taxaid"),
                   c("taxoil"),
                   c("ngo"))
ref.group <- c("aid", "aid", "oil", "non-ngo")
indiv.treats <- list(c("oil", "aid"), 
                     c("tax", "aid"),
                     c("tax", "oil"),
                     c("oil", "tax", "aid", "ngo"))


data.list <- list(unimp.ug, unimp.gh) 

table.results <- c()
country.list <- c("uganda", "ghana")

for(j in 1:length(data.list)){ # This list is defined above, includes full UG + GH unimputed datasets
  for(i in 1:length(var.info$var_name)){ # For each variable in var.info
    for(l in 1:length(treat.vars)){ # For each set of treatment variables (e.g. each element in list)
      dat <- data.list[[j]] ## Set dataset
      which.treats <- indiv.treats[[l]]  ## Set possible treatments
      dv.use <- var.info$var_name[i] ## Pick DV
      treat.use <- treat.vars[[l]]# Grab the treatment variables to include in equation
      ri.eq <- make.equation(dvs = dv.use, ivs = treat.use)  ## Formula for RI
      dat.using <- na.omit(dat[tolower(dat$treat) %in% which.treats, c(treat.use, dv.use)])  ## Exclude NAs from dataset
      declar.eq <- declare_ra(N = nrow(dat.using)) # NOTE: You cannot use the num_arms argument here or conduct_ri will throw an error. 
      declar.eq
      
      ri_out <- conduct_ri(
        formula =  ri.eq,
        assignment = treat.use,
        declaration = declar.eq,
        sharp_hypothesis = 0,
        data = dat.using)
      
      tab <- data.frame(treatment = summary(ri_out)$term, # coefficient name
                        estimate = summary(ri_out)$estimate, # point estimate
                        two_tailed = summary(ri_out)$two_tailed_p_value, # two-tailed p
                        sig = ifelse(summary(ri_out)$two_tailed_p_value < 0.10, "yes", "no"), # signif dummy
                        dv = var.info$var_name[i], var_group = var.info$var_index[i], # dv and var family
                        country = country.list[j], reference_category = ref.group[l]) # country
      table.results <- rbind(tab, table.results)
    }
  }
}


ri.results <- table.results
ri.results <- ri.results[, c("treatment", "estimate", "two_tailed", 
                             "sig", "dv", "reference_category", "var_group", "country")]
rownames(ri.results) <- NULL
ri.results$estimate <- as.numeric(ri.results$estimate)
ri.results$two_tailed <- as.numeric(ri.results$two_tailed)
ri.results[, c("estimate", "two_tailed")] <- round(ri.results[, c("estimate", "two_tailed")], digits = 3)

save(ri.results, file = "results/mass_ri_results_appendix.RData")

#################################################################################
################################ MP RI ANALYSIS #################################
#################################################################################

rm(list = ls()) ## Start from clean workspace

set.seed(1234)
# Function to make equations easier
make.equation <- function(dvs, ivs, interactions = c()){
  return(paste(dvs, "~", paste(c(ivs, interactions), sep = "", collapse = " + ")))
}

packages <- c("ri2","randomizr", "dplyr")
lapply(packages, require, character.only = T)

load("data/ghana_mp_full.RData")
load("data/uganda_mp_full.RData")

## Data and setup

var.info <- read.csv("data/mp_analysis_varinfo.csv", stringsAsFactors = F)

## Create binary variables that will automatically subset out all except 2 treatment conditions
## Note: 1 = oil, 2 = aid, 3 = tax, 4 = ngo


# Ghana
d.gh.mp$oilaid <- dplyr::recode(d.gh.mp$treatment, `1` = "1", `2`= "0", `3`="NA", `4`="NA") # Note that function will spit a warning because "NA" will be treated as NAs by R
d.gh.mp$oilaid <- as.numeric(d.gh.mp$oilaid)
d.gh.mp$taxaid <- dplyr::recode(d.gh.mp$treatment, `1` = "NA", `2`= "0", `3`="1", `4`="NA")
d.gh.mp$taxaid <- as.numeric(d.gh.mp$taxaid)
d.gh.mp$taxoil <- dplyr::recode(d.gh.mp$treatment, `1` = "0", `2` = "NA", `3` = "1", `4` = "NA")
d.gh.mp$taxoil <- as.numeric(d.gh.mp$taxoil) 
d.gh.mp$treat <- dplyr::recode(d.gh.mp$treatment, `1` = "oil", `2`= "aid", `3`="tax", `4`="ngo")

#Uganda
d.ug.mp$oilaid <- dplyr::recode(d.ug.mp$treatment, `1` = "1", `2`= "0", `3`="NA", `4`="NA")
d.ug.mp$oilaid <- as.numeric(d.ug.mp$oilaid)
d.ug.mp$taxaid <- dplyr::recode(d.ug.mp$treatment, `1` = "NA", `2`= "0", `3`="1", `4`="NA")
d.ug.mp$taxaid <- as.numeric(d.ug.mp$taxaid)
d.ug.mp$taxoil <- dplyr::recode(d.ug.mp$treatment, `1` = "0", `2` = "NA", `3` = "1", `4` = "NA")
d.ug.mp$taxoil <- as.numeric(d.ug.mp$taxoil)
d.ug.mp$treat <- dplyr::recode(d.ug.mp$treatment, `1` = "oil", `2`= "aid", `3`="tax", `4`="ngo")

## Specify which subsets and treatment variables you want

treat.vars <- list(c("oilaid"), # Note that these variables are constructed such that the first condition is treatment and second is control, all subjects not in either of these 2 treatments will be NA
                   c("taxaid"),
                   c("taxoil"),
                   c("ngo"))
ref.group <- c("oil", "aid", "oil", "non-ngo")
indiv.treats <- list(c("oil", "aid"), # These should be ordered such that first is the "treatment" and second is "control"
                     c("tax", "aid"),
                     c("tax", "oil"),
                     c("oil", "tax", "aid", "ngo"))

data.list <- list(d.ug.mp, d.gh.mp) 

table.results <- c()
country.list <- c("uganda", "ghana")

for(j in 1:length(data.list)){ # This list is defined above, includes full UG + GH unimputed datasets
  for(i in 1:length(var.info$var_name)){ # For each variable in var.info
    for(l in 1:length(treat.vars)){ # For each set of treatment variables (e.g. each element in list)
      dat <- data.list[[j]] ## Set dataset
      which.treats <- indiv.treats[[l]]  ## Set possible treatments (this determines the subset)
      dv.use <- var.info$var_name[i] ## Pick DV
      treat.use <- treat.vars[[l]]# Grab the treatment variables to include in equation
      ri.eq <- make.equation(dvs = dv.use, ivs = treat.use)  ## Formula for RI
      dat.using <- na.omit(dat[tolower(dat$treat) %in% which.treats, c(treat.use, dv.use)])  ## Exclude NAs from dataset
      declar.eq <- declare_ra(N = nrow(dat.using)) # NOTE: You cannot use the num_arms argument here or conduct_ri will throw an error. 
      declar.eq
      
      ri_out <- conduct_ri(
        formula =  ri.eq,
        assignment = treat.use,
        declaration = declar.eq,
        sharp_hypothesis = 0,
        data = dat.using)
      
      tab <- data.frame(treatment = summary(ri_out)$term, # coefficient name
                        estimate = summary(ri_out)$estimate, # point estimate
                        two_tailed = summary(ri_out)$two_tailed_p_value, # two-tailed p
                        sig = ifelse(summary(ri_out)$two_tailed_p_value < 0.10, "yes", "no"), # signif dummy
                        dv = var.info$var_name[i], var_group = var.info$var_index[i], # dv and var family
                        country = country.list[j], reference_category = ref.group[l]) # country
      table.results <- rbind(tab, table.results)
    }
  }
}


ri.results <- table.results
ri.results <- ri.results[, c("treatment", "estimate", "two_tailed", 
                             "sig", "dv", "reference_category", 
                             "var_group", "country")]
rownames(ri.results) <- NULL
ri.results$estimate <- as.numeric(ri.results$estimate)
ri.results$two_tailed <- as.numeric(ri.results$two_tailed)
ri.results[, c("estimate", "two_tailed")] <- round(ri.results[, c("estimate", "two_tailed")], digits = 3)



save(ri.results, file = "results/mp_ri_results_appendix.RData")


