###############################################
#                                             #
#          Gabriel Moran  Thesis              #
#                                             #
###############################################
# R program to run models for thesis research #
# April 13, 2020                              #
# Gabriel Moran (with Eric Reinhardt)         #
###############################################


###############################################
############### R initialization & setup ######
###############################################
# No need to change any of this unless errors #

# ipak function: install and load multiple R packages.
# check to see if packages are installed. Install them if they are not, then load them into the R session.

ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

# usage
packages <- c("readstata13","stargazer","survival","survminer","tidyverse",
              "dplyr","simPH","foreign","MASS","ggthemes","ggplot2","sandwich",
              "cowplot","psych","interplot","gridExtra","lsmeans")
ipak(packages)
update.packages() # not strictly necessary, disable or say 'no' if source of errors

options(scipen=15) #display quantities in float format, not scientific notation

# enable use of function 'list' for 'ls' ?
rm(list=ls())

# define working directory # GABRIEL <- change as needed
# to find out the local pathname for a file on a Mac...
# ...open the finder to that folder & file...
# ...right click on the file name, then hold down the Option key...
# ...and you'll see "Copy '<filename>' as Pathname" - left-click that
# ...and paste into the quotes below, taking out the final file name,
# ...leaving just the foldername in the quotes
mypath <- "/Users/gabeb/OneDrive/Desktop/Honors Thesis/" # modify accordingly
setwd(mypath)

###############################################
########### import & process dataset     ######
###############################################

"thesis" <- read.csv("ThesisData.csv", header = TRUE, fileEncoding = "UTF-8-BOM") # <-- put in the correct *.csv file name here
data <- thesis
attach(data)
names(data)

########################################
# Summarise data
########################################

#summary statistics for the whole file
# entire dataset
stargazer(data,type="text",summary.stat = c("n","mean","min","max","sd"),summary.logical = FALSE)

# count of available observations per country in the full dataset file
# this assumed that 'statenumber' is the country code variable in your file
as.data.frame(table(data$statenumber))

# descriptive statistics on each variable by country in the full dataset file
describeBy(data,group=data$statenumber,na.rm=TRUE,skew=FALSE,ranges=FALSE,fast=TRUE)

########################################
# Bivariate graphing
########################################

# Establish colors to be used for each country
colors <- c("#f3bb32", "firebrick","springgreen4")
# Establish the R group identifying variable discriminating the countries
grps <- as.factor(data$statenumber)

# DShare against IncumShare1
plot(IncumShare1,DShare,ylab="Incumbent Vote Share Increase",col=colors[grps],cex.lab=1.3,cex=2.1)
      abline(lm(DShare ~ IncumShare1),col="blue",lty="dotted")
      abline(h=0)
      text(DShare ~IncumShare1, labels=countynumber,data=data, cex=0.65)
      legend("topleft", legend=paste(c("Seychelles","Kenya","Nigeria")), pch=16, col=colors)


# DShare against lnChina
plot(lnChina,DShare,ylab="Incumbent Vote Share Increase",
      col = colors[grps],cex.lab=1.3,cex=2.1) + # change cex=1.3 if you want the dots to be smaller
      abline(lm(DShare ~ lnChina),col="#f3bb32",lty="dotted")
      abline(h=0)
      lines(lowess(lnChina, DShare, f=0.33), col="black", lwd=2)
      text(DShare ~lnChina, labels=countynumber,data=data, cex=0.65)
      legend("bottom", legend=paste(c("Seychelles","Kenya","Nigeria")), pch=16, col=colors)

# DShare against lnChinapc
plot(lnChinapc,DShare,ylab="Incumbent Vote Share Increase",
      col = colors[grps],cex.lab=1.3,cex=2.1) + # change cex=1.3 if you want the dots to be smaller
      abline(lm(DShare ~ lnChinapc),col="#f3bb32",lty="dotted")
      abline(h=0)
      lines(lowess(lnChinapc, DShare, f=0.33), col="black", lwd=2)
      text(DShare ~lnChinapc, labels=countynumber,data=data, cex=0.65)
      legend("bottom", legend=paste(c("Seychelles","Kenya","Nigeria")), pch=16, col=colors)


# DShare against lnWB
plot(lnWB,DShare,ylab="Incumbent Vote Share Increase",
      col = colors[grps],cex.lab=1.3,cex=2.1) + # change cex=1.3 if you want the dots to be smaller
      abline(lm(DShare ~ lnWB),col="blue",lty="dotted")
      abline(h=0)
      lines(lowess(lnWB, DShare, f=0.5), col="black", lwd=2)
      text(DShare ~lnWB, labels=countynumber,data=data, cex=0.65)
      legend("bottom", legend=paste(c("Seychelles","Kenya","Nigeria")), pch=16, col=colors)

# graph <- grid.arrange(p1, p2, p3, nrow=1,
#             top="Incumbent Vote Share Change against Covariates")
# ggsave("descriptives.pdf",plot=graph,path=mypath, dpi=300,scale=1) # reactivate this to replace the saved file

# Which countynumbers got the Chinese aid in the 1st place? endogeneity exploration... #

plot(IncumShare1, lnChina,ylab="Ln China Aid",
      col = colors[grps],cex.lab=1.3,cex=2.1) + # change cex=1.3 if you want the dots to be smaller
      abline(lm(lnChina ~ IncumShare1 + IncumShare1*IncumShare1),col="#f3bb32",lty="dotted")
#      abline(h=0)
#      lines(lowess(IncumShare1, lnChina, f=0.33), col="black", lwd=2)
      text(lnChina ~IncumShare1, labels=countynumber,data=data, cex=0.65)
      legend("bottom", legend=paste(c("Seychelles","Kenya","Nigeria")), pch=16, col=colors)

plot(IncumShare1, lnChinapc,ylab="Ln China Aid",
      col = colors[grps],cex.lab=1.3,cex=2.1) + # change cex=1.3 if you want the dots to be smaller
      abline(lm(lnChinapc ~ IncumShare1 + IncumShare1*IncumShare1),col="#f3bb32",lty="dotted")
#      abline(h=0)
#      lines(lowess(IncumShare1, lnChinapc, f=0.33), col="black", lwd=2)
      text(lnChinapc ~IncumShare1, labels=countynumber,data=data, cex=0.65)
      legend("bottom", legend=paste(c("Seychelles","Kenya","Nigeria")), pch=16, col=colors)

###############################################
########### ESTIMATE MAIN MODELS         ######
###############################################

### Model 1: dep var = DShare, cross-country fixed effects in 'statenumber'

# base model
model1 <- lm(DShare ~ state
               + IncumShare1
               + lnChina + lnWB,
               data=data, na.action=na.exclude)
stargazer(model1,type="text")

# test of equality of coefficients
library(car) # initialize relevant R library
linearHypothesis(model1, c("lnChina = lnWB")) # joint test of null hypotheses listed under c()

# base model per capita
model1pc <- lm(DShare ~ state
                 + IncumShare1
                 + lnChinapc + lnWBpc,
                 data=data, na.action=na.exclude)
stargazer(model1pc,type="text")

# test of equality of coefficients
# library(car) # initialize relevant R library
linearHypothesis(model1pc, c("lnChinapc = lnWBpc")) # joint test of null hypotheses listed under c()

# base model per capita
model2pc <- lm(DShare ~ state
                 + IncumShare1
                 + China + WB
                 + lnChinapc + lnWBpc,
                 data=data, na.action=na.exclude)
stargazer(model2pc,type="text")

# test of equality of coefficients
# library(car) # initialize relevant R library
linearHypothesis(model2pc, c("lnChinapc = lnWBpc")) # joint test of null hypotheses listed under c()


# does lnChina work differently for the three countries?
model2 <- lm(DShare ~ state
                 + IncumShare1
                 + lnChina + lnWB
                 + state*lnChina + state*lnWB,
                 data=data, na.action=na.exclude)
stargazer(model2,type="text")
# A: it appears to be somewhat positive for Sey & Nig but zero-ish for Kenya


# try dropping observations... e.g., Kenya cases
noK <- subset(data,statenumber<2 | statenumber>2)
summary(noK)
model3 <- lm(DShare ~ state
                 + IncumShare1
                 + lnChina + lnWB,
                 data=noK, na.action=na.exclude)
stargazer(model3,type="text")

# test of equality of coefficients
linearHypothesis(model3, c("lnChina = lnWB")) # joint test of null hypotheses listed under c()

##Robustness Test Control for Capital
#create subset with no capital districts
nocap <-subset(data,Capital<1 | Capital>1)
summary(nocap)

##run both estimate models using nocap
model4 <- lm(DShare ~ state
             + IncumShare1
             + lnChina + lnWB,
             data=nocap, na.action=na.exclude)
stargazer(model4,type="text")

# test of equality of coefficients
linearHypothesis(model4, c("lnChina = lnWB")) # joint test of null hypotheses listed under c()

model4pc <- lm(DShare ~ state
             + IncumShare1
             + lnChinapc + lnWBpc,
             data=nocap, na.action=na.exclude)
stargazer(model4pc,type="text")


# test of equality of coefficients
linearHypothesis(model4pc, c("lnChinapc = lnWBpc")) # joint test of null hypotheses listed under c()

##Robustness Test Country specific Slopes Model
#write out linear model as follows "m.interaction <-lm(DShare ~ lnChinapc*statenumber, data=data)"

m.interaction <-lm(DShare ~ lnChinapc*statenumber, data=data)
#pull up ANOVA table
anova(m.interaction)

#Obtain slopes
m.interaction$coefficients
m.lst <- lstrends(m.interaction, "statenumber", var="lnChinapc")