###########################################################################
# PUBLG100: Introduction to Quantitative Methods
#
# Week 5 Seminar: Multiple linear regression models
#
#

# 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())

## ------------------------------------------------------------------------
# load dataset in Stata format
world_data <- read.dta("qog_std_cs_jan15.dta")

# check the dimensions of the dataset
dim(world_data)

## ------------------------------------------------------------------------
world_data <- select(world_data,
                     country = cname, 
                     political_stability = wbgi_pse,
                     latitude = lp_lat_abst, 
                     globalization = dr_ig, 
                     inst_quality = ti_cpi)

## ------------------------------------------------------------------------
head(world_data)

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

## ------------------------------------------------------------------------
world_data <- filter(world_data, 
                     !is.na(latitude), 
                     !is.na(globalization), 
                     !is.na(inst_quality))

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

## ------------------------------------------------------------------------
# transform latitude variable
world_data$latitude <- world_data$latitude * 90

## ------------------------------------------------------------------------
plot(political_stability ~ latitude, data = world_data)

## ------------------------------------------------------------------------
latitude_model <- lm(political_stability ~ latitude, data = world_data)

# add the line
plot(political_stability ~ latitude, data = world_data)
abline(latitude_model)

# regression output
screenreg(latitude_model)

## ------------------------------------------------------------------------
# model with more explanatory variables
inst_model <- lm(political_stability ~ latitude + globalization + inst_quality, 
                 data = world_data)

screenreg(list(latitude_model, inst_model))

## ------------------------------------------------------------------------
anova(latitude_model, inst_model)

## ------------------------------------------------------------------------
# multivariate regression model
z.out <- zelig(political_stability ~ latitude + globalization + inst_quality, 
               data = world_data, 
               model = "ls")

# setting covariates
x.out <- setx(z.out, inst_quality = seq(0, 10, 1))

# simulation
s.out <- sim(z.out, x = x.out)

# plot results
ci.plot(s.out, ci = 95, xlab = "Quality of Institutions")

## ----fig.width=13, fig.height=11-----------------------------------------
# 25th percentile of quality of institutions
x.out.low <- setx(z.out, inst_quality = quantile(world_data$inst_quality, 0.25))

# 75th percentile of quality of institutions
x.out.high <- setx(z.out, inst_quality = quantile(world_data$inst_quality, 0.75))

# simulation
s.out <- sim(z.out, x = x.out.low, x1 = x.out.high)

# summary output
summary(s.out)