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

## ----message=FALSE-------------------------------------------------------
library(Zelig) # to predict from regression output and present uncertainty
library(texreg) # to plot regression output to screen or file

## ----eval=FALSE----------------------------------------------------------
## # clear workspace
## rm(list = ls() )
## 
## # set working directory
## setwd("N:/PUBLG100")

## ------------------------------------------------------------------------
# load csv data set
cor.data <- read.csv("http://uclspp.github.io/PUBLG100/data/corruption.csv")

## ------------------------------------------------------------------------
summary(cor.data)

## ------------------------------------------------------------------------
# estimate linear multivariate regression model
m1 <- lm(gdp ~ ti.cpi + region, data = cor.data)

# print model output to word file
htmlreg(m1, file = "model1.doc")

# print output to screen
screenreg(m1)

## ------------------------------------------------------------------------
# estimate linear model excluding regions
m.noregions <- lm(gdp ~ ti.cpi, data = cor.data)

# joint hypothesis testing
anova(m.noregions, m1)

## ------------------------------------------------------------------------
# estimate the linar model including regions
linear.model <- zelig( gdp ~ ti.cpi + region, data = cor.data, model = "ls")

# look at the ti.cpi variable to see the range
summary( cor.data$ti.cpi)

# set the covariates
covariates <- setx( linear.model, ti.cpi = 1:10)

# simulate
predictions <- sim( linear.model, covariates)

# plot results
plot(predictions)

## ----fig.width=13, fig.height=11-----------------------------------------
# estimate the linear model including regions
linear.model <- zelig( gdp ~ ti.cpi + region, data = cor.data, model = "ls")

# look at the labels of the factor variable regions
table(cor.data$region)

# set the covariates
europa <- setx( linear.model, region = "Europe")
amerika <- setx( linear.model, region = "Americas")

# simulate
predictions <- sim( linear.model, x = europa, x1 = amerika)

# look at the first difference
summary(predictions)

# you can do this visually, too
plot(predictions)