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

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

## ----eval = FALSE--------------------------------------------------------
## # Change your working directory
## setwd("N:/PUBLG100")
## 
## # Check your working directory
## getwd()
## 
## # clear the environment
## rm(list = ls())

## ------------------------------------------------------------------------
CA_benefits <- read.csv("CABenefits.csv")
CA_students <- read.csv("CAStudents.csv")

head(CA_benefits)
head(CA_students)

## ------------------------------------------------------------------------
CA_admin <- read.csv("CASchoolAdmin.csv")

CA_schools <- merge(CA_admin, CA_students, by = c("county", "school"))
CA_schools <- merge(CA_schools, CA_benefits, by = "district")

head(CA_schools)

## ------------------------------------------------------------------------
CA_schools <- mutate(CA_schools, score = (math + read) / 2)

## ------------------------------------------------------------------------
calworks_model <- lm(score ~ calworks, data = CA_schools)

summary(calworks_model)

screenreg(calworks_model)

## ------------------------------------------------------------------------
plot(score ~ calworks, 
     data = CA_schools,
     xlab = "Percent qualifying for CalWorks income assistance",
     ylab = "Test Scores")

abline(calworks_model, col = "red")

## ------------------------------------------------------------------------
CA_schools <- mutate(CA_schools, STR = students/teachers)

## ------------------------------------------------------------------------
STR_model <- lm(score ~ STR, data = CA_schools)

income_model <- lm(score ~ income, data = CA_schools)

screenreg(list(STR_model, income_model, calworks_model))

## ------------------------------------------------------------------------
htmlreg(list(STR_model, income_model, calworks_model), 
        file = "solutions4_model_comparison.doc")

## ------------------------------------------------------------------------
range(CA_schools$calworks)

## ------------------------------------------------------------------------
z.out <- zelig(score ~ calworks, data = CA_schools, model = "ls")
x.out <- setx(z.out, calworks = seq(0, 100, 0.5))
s.out <- sim(z.out, x = x.out)

ci.plot(s.out, ci = 95, xlab = "Percent qualifying for CalWorks income assistance")

## ----eval = FALSE--------------------------------------------------------
## png(filename = "solutions4_calworks_model.png")
## 
## plot(score ~ calworks,
##      data = CA_schools,
##      xlab = "Percent qualifying for CalWorks income assistance",
##      ylab = "Test Scores")
## abline(calworks_model, col = "red")
## 
## dev.off()

## ----eval = FALSE--------------------------------------------------------
## png(filename = "solutions4_calworks_model_zelig.png")
## 
## ci.plot(s.out,
##         ci = 95,
##         xlab = "Percent qualifying for CalWorks income assistance",
##         ylab = "Expected Value of Test Score")
## 
## dev.off()