###########################################################################
# PUBLG100: Introduction to Quantitative Methods
#
# Week 7 Seminar: Interactions
#
#

# Set your working directory
#
# CAUTION: Make sure the directory you specify here matches the working directory on your computer.
#          We're using N:/PUBLG100 only for illustration purposes and it would only work if you're 
#          using a UCL dekstop. If you're using your own laptop, then replace N:/PUBLG100 with the 
#          appropriate directory (or folder)
setwd("N:/PUBLG100")


# Verify that your working directory is set correctly
getwd()


## ----message=FALSE-------------------------------------------------------
library(foreign)
library(Zelig)
library(texreg)
library(dplyr)

## ------------------------------------------------------------------------
rm(list = ls())

## ------------------------------------------------------------------------
world_data <- read.dta("QoG2012.dta")

names(world_data)

## ------------------------------------------------------------------------
world_data <- select(world_data, 
                     hdi = undp_hdi,
                     corruption_control = wbgi_cce,
                     gdp = wdi_gdpc,
                     former_col)

head(world_data)

## ------------------------------------------------------------------------
summary(world_data) 

## ------------------------------------------------------------------------
table(world_data$former_col) 

## ------------------------------------------------------------------------
plot(hdi ~ corruption_control, 
     data = world_data, 
     col = as.factor(former_col), 
     pch = 20)

# add a legend
legend("bottomright", 
       c("former_col = 0", "former_col = 1"), 
       col=c("black","red"), 
       pch = 20)

## ------------------------------------------------------------------------
model1 <- lm(hdi ~ corruption_control + former_col, data = world_data)

screenreg(model1)

## ---- echo = FALSE-------------------------------------------------------
# ignore the following code, it's used for generating the content on the webpage
get_intercept_eq <- function(model, i, dummy) {
  coefs <- coefficients(model)
  sprintf("`Intercept` + (`%s` * %d)<br>= %.2f + (%.2f * %d)<br>= %.2f", names(coefs)[i], dummy, coefs[1], coefs[i], dummy, coefs[1] + (coefs[i] * dummy))
}

get_slope_eq <- function(model, i) {
  coefs <- coefficients(model)
  sprintf("`%s`<br> = %.2f", names(coefs)[i], coefs[i])
}


## ------------------------------------------------------------------------
model_coef <- coefficients(model1)
model_coef

## ------------------------------------------------------------------------
plot(hdi ~ corruption_control, data = world_data,
     col = as.factor(former_col), 
     pch = 20,
     xlab = "Corruption Control",
     ylab = "Human Development Index",
     main = "Effect of a Dummy Variable")

# the regression line when former_col = 0
abline(model_coef[1], model_coef[2], col = "black")

# the regression line when former_col = 1
abline(model_coef[1] + model_coef[3], model_coef[2], col = "red")

# add a legend to the plot
legend("bottomright", 
       c("former_col = 0", "former_col = 1"), 
       col=c("black","red"), 
       lwd = 1,   # line width = 1 for adding a line to the legend
       pch = 20)

## ------------------------------------------------------------------------
model2 <- lm(hdi ~ corruption_control + former_col + corruption_control:former_col, 
             data = world_data)

screenreg(list(model1, model2))
anova(model1, model2)


## ------------------------------------------------------------------------
z.out <- zelig(hdi ~ corruption_control + former_col + corruption_control:former_col, 
               data = world_data, 
               model = "ls")

## ------------------------------------------------------------------------
range(world_data$corruption_control, na.rm = TRUE)

## ------------------------------------------------------------------------
# set covariates for countries that weren't colonised
x.out1 <- setx(z.out, former_col = 0, corruption_control = -2:3)

# set covariates for colonised countries
x.out2 <- setx(z.out, former_col = 1, corruption_control = -2:3)

# simulate 
s.out <- sim(z.out, x = x.out1, x1 = x.out2)

## ------------------------------------------------------------------------
ci.plot(s.out, 
        ci = 95,
        xlab = "Control of Corruption",
        ylab = "Expected Value of Human Development Index",
        main = "Effect of Control of Corruption by Colonial Past"
        ) 

# add text to Zelig plot. 
# NOTE: x and y are just plot coordinates so you can change them to control where 
#       the text is displayed
text(x = -1, y = 0.8, labels = "former_col = 0")
text(x = 1, y = 0.7, labels = "former_col = 1")

## ------------------------------------------------------------------------
plot(hdi ~ gdp, data = world_data)

## ------------------------------------------------------------------------
model1 <- lm(hdi ~ gdp, data = world_data)
summary(model1)

plot(hdi ~ gdp, data = world_data)
abline(model1)

## ------------------------------------------------------------------------
model2 <- lm(hdi ~ gdp + I(gdp^2), data = world_data)
summary(model2)

## ------------------------------------------------------------------------
anova(model1, model2)

## ------------------------------------------------------------------------
z.out <- zelig(hdi ~ gdp + I(gdp^2), data = world_data, model = "ls")

# setting covariates; GDP/captia is a sequence from 0 to 60000 in steps of 1000
x.out <- setx(z.out, gdp = seq(0, 60000, 1000))

# simulate using our model and our covariates
s.out <- sim(z.out, x = x.out)

# plot the results
ci.plot(s.out, ci = 95)