rm(list = ls())

############################## MASS - Randomization Inference #############################


rm(list = ls())
require(stargazer)
require(starpolishr)
require(Hmisc)


## 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 = " + ")))
}


results <- load("results/mass_ri_results_appendix.RData")
load("data/uganda_citizen_imputed_full.RData")

dumMod <- subset(ug.imp.full, imp==0) #Create dummy host model as a placeholder 

varInfo <- read.csv("data/mass_analysis_varinfo.csv", stringsAsFactors = F)
final.results <- ri.results



###### Action Measures, Main Effects Unimputed (RI)
plotVars <- varInfo$var_name[varInfo$include_plot == 1 & varInfo$var_index=="action"]
var_group <- c(plotVars) # Include all variable names you want to report
names(var_group) <- varInfo$var_label[varInfo$include_plot == 1 & varInfo$var_index=="action"]

oilaid <- final.results[final.results$treatment %in% c("oilaid")
                        & final.results$var_group == "action"
                        & final.results$dv %in% plotVars, ]

taxaid <- final.results[final.results$treatment %in% c("taxaid")
                        & final.results$var_group == "action"
                        & final.results$dv %in% plotVars, ]

taxoil <- final.results[final.results$treatment %in% c("taxoil")
                        & final.results$var_group == "action"
                        & final.results$dv %in% plotVars, ]

res.list <- list(oilaid[oilaid$country == "uganda", ], 
                 oilaid[oilaid$country == "ghana", ],
                 taxaid[taxaid$country == "uganda", ],
                 taxaid[taxaid$country == "ghana", ],
                 taxoil[taxoil$country == "uganda", ],
                 taxoil[taxoil$country == "ghana", ])

coef.list <- list()
for(j in 1:length(res.list)){
  coef.loop <- c()
  for(i in 1:length(var_group)){
    coef.loop <- c(coef.loop, res.list[[j]]$est[res.list[[j]]$dv == var_group[i]])
  }
  names(coef.loop) <- var_group
  coef.list[[j]] <- coef.loop
}


p.list <- list()
for(j in 1:length(res.list)){
  p.loop <- c()
  for(i in 1:length(var_group)){
    p.loop <- c(p.loop, res.list[[j]]$two_tailed[res.list[[j]]$dv == var_group[i]])
  }
  names(p.loop) <- var_group
  p.list[[j]] <- p.loop
}

host.model <- lm(make.equation(dvs = "type_high", ivs = var_group), data = dumMod) # -1 to account for intercept in host model

# Table S78. Action Index (RI)
tab.s78 <- stargazer(host.model, host.model, host.model, host.model, host.model, host.model, style = "APSR", 
          colnames = T, omit = "Constant", covariate.labels = c(names(var_group)), 
          coef = coef.list,
          se = p.list, 
          p = p.list,
          intercept.bottom = F,
          column.separate = c(2, 2, 2),
          column.labels = c("Oil Treat-Aid Ctr", "Tax Treat-Aid Ctr", "Tax Treat - Oil Ctr"), 
          dep.var.labels.include = T, 
          dep.var.labels = c("Uganda", "Ghana", "Uganda", "Ghana", "Uganda", "Ghana"), 
          label = "tab:action:ri", 
          omit.stat = c("all"),
          title  = "Action Index, Randomization Inference",
          multicolumn = F, 
          model.numbers = T)

tab.s78 <- star_notes_tex(tab.s78, note.type = "threeparttable", 
                          note ="Estimates are based on randomization inference with two-tailed p-values in parentheses. Columns 1, 3, and 5 report the Ugandan results while Columns 2, 4, and 6 report the Ghanaian estimates.") 

cat(tab.s78, file = paste0(getwd(),"/tables/s78.tex"),sep="\n") 



###### Action Measures, NGO Main Effects Unimputed (RI)
plotVars <- varInfo$var_name[varInfo$include_plot == 1 & varInfo$var_index=="action"]
var_group <- c(plotVars) # Include all variable names you want to report
names(var_group) <- varInfo$var_label[varInfo$include_plot == 1 & varInfo$var_index=="action"]

ngonon <- final.results[final.results$treatment %in% c("ngo")
                        & final.results$var_group == "action"
                        & final.results$dv %in% plotVars, ]


res.list <- list(ngonon[ngonon$country == "uganda", ], 
                 ngonon[ngonon$country == "ghana", ])

coef.list <- list()
for(j in 1:length(res.list)){
  coef.loop <- c()
  for(i in 1:length(var_group)){
    coef.loop <- c(coef.loop, res.list[[j]]$est[res.list[[j]]$dv == var_group[i]])
  }
  names(coef.loop) <- var_group
  coef.list[[j]] <- coef.loop
}


p.list <- list()
for(j in 1:length(res.list)){
  p.loop <- c()
  for(i in 1:length(var_group)){
    p.loop <- c(p.loop, res.list[[j]]$two_tailed[res.list[[j]]$dv == var_group[i]])
  }
  names(p.loop) <- var_group
  p.list[[j]] <- p.loop
}

# Table S79. Action Index, NGO Treatment, non-NGO Control (RI)
host.model <- lm(make.equation(dvs = "type_high", ivs = var_group), data = dumMod) # -1 to account for intercept in host model
tab.s79 <- stargazer(host.model, host.model, style = "APSR", 
          colnames = T, omit = "Constant", covariate.labels = c(names(var_group)), 
          coef = coef.list,
          se = p.list, 
          p = p.list,
          intercept.bottom = F,
          column.separate = c(2),
          column.labels = c("NGO Treatment, Non-NGO Control"), 
          title  = "Action Index, NGO Pooled Results, Randomization Inference",
          dep.var.labels.include = T, 
          dep.var.labels = c("Uganda", "Ghana"), 
          label = "tab:action:ngo:ri", 
          omit.stat = c("all"),
          multicolumn = F, 
          model.numbers = T)

tab.s79 <- star_notes_tex(tab.s79, note.type = "threeparttable", 
  note ="Estimates are based on randomization inference with two-tailed p-values in parentheses. Column 1 reports the Ugandan results while Column 2 reports the Ghanaian estimates.") 

cat(tab.s79, file = paste0(getwd(),"/tables/s79.tex"),sep="\n") 


###### Benefit Measures, Main Effects Unimputed (RI)
plotVars <- varInfo$var_name[varInfo$include_plot == 1 & varInfo$var_index=="benefit"]
var_group <- c(plotVars) # Include all variable names you want to report
names(var_group) <- varInfo$var_label[varInfo$include_plot == 1 & varInfo$var_index=="benefit"]

oilaid <- final.results[final.results$treatment %in% c("oilaid")
                        & final.results$var_group == "benefit"
                        & final.results$dv %in% plotVars, ]

taxaid <- final.results[final.results$treatment %in% c("taxaid")
                        & final.results$var_group == "benefit"
                        & final.results$dv %in% plotVars, ]

taxoil <- final.results[final.results$treatment %in% c("taxoil")
                        & final.results$var_group == "benefit"
                        & final.results$dv %in% plotVars, ]

res.list <- list(oilaid[oilaid$country == "uganda", ], 
                 oilaid[oilaid$country == "ghana", ],
                 taxaid[taxaid$country == "uganda", ],
                 taxaid[taxaid$country == "ghana", ],
                 taxoil[taxoil$country == "uganda", ],
                 taxoil[taxoil$country == "ghana", ])

coef.list <- list()
for(j in 1:length(res.list)){
  coef.loop <- c()
  for(i in 1:length(var_group)){
    coef.loop <- c(coef.loop, res.list[[j]]$est[res.list[[j]]$dv == var_group[i]])
  }
  names(coef.loop) <- var_group
  coef.list[[j]] <- coef.loop
}


p.list <- list()
for(j in 1:length(res.list)){
  p.loop <- c()
  for(i in 1:length(var_group)){
    p.loop <- c(p.loop, res.list[[j]]$two_tailed[res.list[[j]]$dv == var_group[i]])
  }
  names(p.loop) <- var_group
  p.list[[j]] <- p.loop
}

host.model <- lm(make.equation(dvs = "type_high", ivs = var_group), data = dumMod) # -1 to account for intercept in host model

# Table S80. Benefit Index (RI)
tab.s80 <- stargazer(host.model, host.model, host.model, host.model, host.model, host.model, style = "APSR", 
          colnames = T, omit = "Constant", covariate.labels = c(names(var_group)), 
          coef = coef.list,
          se = p.list, 
          p = p.list,
          intercept.bottom = F,
          column.separate = c(2, 2, 2),
          column.labels = c("Oil Treat-Aid Ctr", "Tax Treat-Aid Ctr", "Tax Treat - Oil Ctr"), 
          dep.var.labels.include = T, 
          dep.var.labels = c("Uganda", "Ghana", "Uganda", "Ghana", "Uganda", "Ghana"), 
          label = "tab:benefit:ri", 
          omit.stat = c("all"),
          title  = "Benefit Index, Randomization Inference",
          multicolumn = F, 
          model.numbers = T)

tab.s80 <- star_notes_tex(tab.s80, note.type = "threeparttable", 
                          note ="Estimates are based on randomization inference with two-tailed p-values in parentheses. Columns 1, 3, and 5 report the Ugandan results while Columns 2, 4, and 6 report the Ghanaian estimates.") 

cat(tab.s80, file = paste0(getwd(),"/tables/s80.tex"),sep="\n") 




###### Benefit Measures, NGO Main Effects Unimputed (RI)
plotVars <- varInfo$var_name[varInfo$include_plot == 1 & varInfo$var_index=="benefit"]
var_group <- c(plotVars) # Include all variable names you want to report
names(var_group) <- varInfo$var_label[varInfo$include_plot == 1 & varInfo$var_index=="benefit"]

ngonon <- final.results[final.results$treatment %in% c("ngo")
                        & final.results$var_group == "benefit"
                        & final.results$dv %in% plotVars, ]


res.list <- list(ngonon[ngonon$country == "uganda", ], 
                 ngonon[ngonon$country == "ghana", ])

coef.list <- list()
for(j in 1:length(res.list)){
  coef.loop <- c()
  for(i in 1:length(var_group)){
    coef.loop <- c(coef.loop, res.list[[j]]$est[res.list[[j]]$dv == var_group[i]])
  }
  names(coef.loop) <- var_group
  coef.list[[j]] <- coef.loop
}


p.list <- list()
for(j in 1:length(res.list)){
  p.loop <- c()
  for(i in 1:length(var_group)){
    p.loop <- c(p.loop, res.list[[j]]$two_tailed[res.list[[j]]$dv == var_group[i]])
  }
  names(p.loop) <- var_group
  p.list[[j]] <- p.loop
}

# Table S81. Benefit Index, NGO Treatment, non-NGO Control (RI)
host.model <- lm(make.equation(dvs = "type_high", ivs = var_group), data = dumMod) # -1 to account for intercept in host model
tab.s81 <- stargazer(host.model, host.model, style = "APSR", 
          colnames = T, omit = "Constant", covariate.labels = c(names(var_group)), 
          coef = coef.list,
          se = p.list, 
          p = p.list,
          intercept.bottom = F,
          column.separate = c(2),
          column.labels = c("NGO Treatment, Non-NGO Control"), 
          title  = "Benefit Index, NGO Pooled Results, Randomization Inference",
          dep.var.labels.include = T, 
          dep.var.labels = c("Uganda", "Ghana"), 
          label = "tab:benefit:ngo:ri", 
          omit.stat = c("all"),
          multicolumn = F, 
          model.numbers = T)

tab.s81 <- star_notes_tex(tab.s81, note.type = "threeparttable", 
                          note ="Estimates are based on randomization inference with two-tailed p-values in parentheses. Column 1 reports the Ugandan results while Column 2 reports the Ghanaian estimates.") 

cat(tab.s81, file = paste0(getwd(),"/tables/s81.tex"),sep="\n") 


####################################### Misappropriation   ################################################

###### Misappropriation Measures, Main Effects Unimputed (RI)
plotVars <- varInfo$var_name[varInfo$include_plot == 1 & varInfo$var_index=="misappropriation"]
var_group <- c(plotVars) # Include all variable names you want to report
names(var_group) <- varInfo$var_label[varInfo$include_plot == 1 & varInfo$var_index=="misappropriation"]

oilaid <- final.results[final.results$treatment %in% c("oilaid")
                        & final.results$var_group == "misappropriation"
                        & final.results$dv %in% plotVars, ]

taxaid <- final.results[final.results$treatment %in% c("taxaid")
                        & final.results$var_group == "misappropriation"
                        & final.results$dv %in% plotVars, ]

taxoil <- final.results[final.results$treatment %in% c("taxoil")
                        & final.results$var_group == "misappropriation"
                        & final.results$dv %in% plotVars, ]

res.list <- list(oilaid[oilaid$country == "uganda", ], 
                 oilaid[oilaid$country == "ghana", ],
                 taxaid[taxaid$country == "uganda", ],
                 taxaid[taxaid$country == "ghana", ],
                 taxoil[taxoil$country == "uganda", ],
                 taxoil[taxoil$country == "ghana", ])

coef.list <- list()
for(j in 1:length(res.list)){
  coef.loop <- c()
  for(i in 1:length(var_group)){
    coef.loop <- c(coef.loop, res.list[[j]]$est[res.list[[j]]$dv == var_group[i]])
  }
  names(coef.loop) <- var_group
  coef.list[[j]] <- coef.loop
}


p.list <- list()
for(j in 1:length(res.list)){
  p.loop <- c()
  for(i in 1:length(var_group)){
    p.loop <- c(p.loop, res.list[[j]]$two_tailed[res.list[[j]]$dv == var_group[i]])
  }
  names(p.loop) <- var_group
  p.list[[j]] <- p.loop
}


host.model <- lm(make.equation(dvs = "type_high", ivs = var_group), data = dumMod) # -1 to account for intercept in host model

#Table S82. Misappropriation Index (RI)
tab.s82 <- stargazer(host.model, host.model, host.model, host.model, host.model, host.model, style = "APSR", 
          colnames = T, omit = "Constant", covariate.labels = c(names(var_group)), 
          coef = coef.list,
          se = p.list, 
          p = p.list,
          intercept.bottom = F,
          column.separate = c(2, 2, 2),
          column.labels = c("Oil Treat-Aid Ctr", "Tax Treat-Aid Ctr", "Tax Treat - Oil Ctr"), 
          dep.var.labels.include = T, 
          dep.var.labels = c("Uganda", "Ghana", "Uganda", "Ghana", "Uganda", "Ghana"), 
          label = "tab:misapp:ri", 
          omit.stat = c("all"),
          title  = "Misappropriation Index, Randomization Inference",
          multicolumn = F, 
          model.numbers = T)


tab.s82 <- star_notes_tex(tab.s82, note.type = "threeparttable", 
                          note ="Estimates are based on randomization inference with two-tailed p-values in parentheses. Columns 1, 3, and 5 report the Ugandan results while Columns 2, 4, and 6 report the Ghanaian estimates.") 

cat(tab.s82, file = paste0(getwd(),"/tables/s82.tex"),sep="\n") 


###### misappropriation Measures, NGO Main Effects Unimputed (RI)
plotVars <- varInfo$var_name[varInfo$include_plot == 1 & varInfo$var_index=="misappropriation"]
var_group <- c(plotVars) # Include all variable names you want to report
names(var_group) <- varInfo$var_label[varInfo$include_plot == 1 & varInfo$var_index=="misappropriation"]

ngonon <- final.results[final.results$treatment %in% c("ngo")
                        & final.results$var_group == "misappropriation"
                        & final.results$dv %in% plotVars, ]


res.list <- list(ngonon[ngonon$country == "uganda", ], 
                 ngonon[ngonon$country == "ghana", ])

coef.list <- list()
for(j in 1:length(res.list)){
  coef.loop <- c()
  for(i in 1:length(var_group)){
    coef.loop <- c(coef.loop, res.list[[j]]$est[res.list[[j]]$dv == var_group[i]])
  }
  names(coef.loop) <- var_group
  coef.list[[j]] <- coef.loop
}


p.list <- list()
for(j in 1:length(res.list)){
  p.loop <- c()
  for(i in 1:length(var_group)){
    p.loop <- c(p.loop, res.list[[j]]$two_tailed[res.list[[j]]$dv == var_group[i]])
  }
  names(p.loop) <- var_group
  p.list[[j]] <- p.loop
}


host.model <- lm(make.equation(dvs = "type_high", ivs = var_group), data = dumMod) # -1 to account for intercept in host model

# Table S83. Misappropriation Index, NGO Treatment, non-NGO Control (RI)
tab.s83 <- stargazer(host.model, host.model, style = "APSR", 
          colnames = T, omit = "Constant", covariate.labels = c(names(var_group)), 
          coef = coef.list,
          se = p.list, 
          p = p.list,
          intercept.bottom = F,
          column.separate = c(2),
          column.labels = c("NGO Treatment, Non-NGO Control"), 
          title  = "Misappropriation Index, NGO Pooled Results, Randomization Inference",
          dep.var.labels.include = T, 
          dep.var.labels = c("Uganda", "Ghana"), 
          label = "tab:misapp:ngo:ri", 
          omit.stat = c("all"),
          multicolumn = F, 
          model.numbers = T)

tab.s83 <- star_notes_tex(tab.s83, note.type = "threeparttable", 
                          note ="Estimates are based on randomization inference with two-tailed p-values in parentheses. Column 1 reports the Ugandan results while Column 2 reports the Ghanaian estimates.") 

cat(tab.s83, file = paste0(getwd(),"/tables/s83.tex"),sep="\n") 


