###########################################################################
# 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)
= %.2f + (%.2f * %d)
= %.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`
= %.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)