###########################################################################
# PUBLG100: Introduction to Quantitative Methods
#
# Week 4 Seminar: Bivariate linear regression models
#
#

# Change your working directory
setwd("N:/PUBLG100")

# Check your working directory
getwd()


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

## ------------------------------------------------------------------------
install.packages("https://cran.r-project.org/src/contrib/Archive/Zelig/Zelig_4.2-1.tar.gz", 
                 repos=NULL,
                 type="source")

install.packages("texreg")

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

## ------------------------------------------------------------------------
# load the communities datasets
communities <- read.csv("http://uclspp.github.io/PUBLG100/data/communities.csv")
communities_employment <- read.csv("http://uclspp.github.io/PUBLG100/data/communities_employment.csv")

## ------------------------------------------------------------------------
names(communities)
str(communities)

## ------------------------------------------------------------------------
str(communities_employment)

## ------------------------------------------------------------------------
# merge the two datasets
communities <- merge(communities, communities_employment, by = c("state", "communityname"))

# explore dataset
names(communities)

## ----eval = FALSE--------------------------------------------------------
## View(communities)

## ------------------------------------------------------------------------
communities <- select(communities, 
                      state, 
                      Community = communityname, 
                      UnemploymentRate = PctUnemployed, 
                      NoHighSchool = PctNotHSGrad,
                      White = racePctWhite)

## ------------------------------------------------------------------------
plot(communities$NoHighSchool, communities$UnemploymentRate,
     xlab = "Adults without High School education (%)",
     ylab = "Unemployment Rate")

## ------------------------------------------------------------------------
model1 <- lm(UnemploymentRate ~ NoHighSchool, data = communities)

## ------------------------------------------------------------------------
summary(model1)

## ------------------------------------------------------------------------
plot(communities$NoHighSchool, communities$UnemploymentRate,
     xlab = "Adults without High School education (%)",
     ylab = "Unemployment Rate")
abline(model1, col = "red")

## ------------------------------------------------------------------------
screenreg(model1)

## ------------------------------------------------------------------------
communities$Minority <- 1 - communities$White

## ------------------------------------------------------------------------
model2 <- lm(UnemploymentRate ~ Minority, data = communities)
summary(model2)

## ------------------------------------------------------------------------
plot(communities$Minority, communities$UnemploymentRate,
     xlab = "Minority population (%)",
     ylab = "Unemployment Rate")
abline(model2, col = "blue")

## ------------------------------------------------------------------------
screenreg(list(model1, model2))

## ------------------------------------------------------------------------
htmlreg(list(model1, model2), file = "lab4_model_comparison.doc")

## ------------------------------------------------------------------------
z.out <- zelig(UnemploymentRate ~ Minority, data = communities, model = "ls")

## ------------------------------------------------------------------------
summary(z.out)

## ------------------------------------------------------------------------
x.out <- setx(z.out, Minority = seq(0, 1, 0.1))
s.out <- sim(z.out, x = x.out)

summary(s.out)
plot(s.out, qi = "pv", xlab = "Minority population (%)")