library(AER)
data("CASchools")
head(CASchools)

#Creating Variables
CASchools$STR <- CASchools$students/CASchools$teachers
summary(CASchools$STR)

CASchools$score <- (CASchools$read + CASchools$math)/2
summary(CASchools$score)

#Changing names
dat <- CASchools

summary(lm(score ~ income, data=dat))

library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 1), se = FALSE)
p

#YOU MUST USE THE RAW=T OPTION WHEN YOU USE THE POLY FUNCTION
#OR IT WILL NOT WORK

m <- lm(score ~ poly(income, 2, raw=T), data=dat )
summary(m)

dat$income2  <- dat$income^2

m <- lm(score ~ income + income2, data=dat )
summary(m)

#Plotting this relationship:
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 2), se = FALSE)
p

summary(m)

coef(m)[2] + 2*coef(m)[3]*20

#and at 30
coef(m)[2] + 2*coef(m)[3]*30

#and at 50
coef(m)[2] + 2*coef(m)[3]*50

p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 2), se = FALSE)
p

#Plotting this relationship:
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 4), se = FALSE)
p

m <- lm(score ~ poly(income, 4, raw=T), data=dat )
coef(m)

#20
coef(m)[2] + 2*coef(m)[3]*20 + 3*coef(m)[4]*20^2 + 4*coef(m)[5]*20^3

#and at 30
coef(m)[2] + 2*coef(m)[3]*30 + 3*coef(m)[4]*30^2 + 4*coef(m)[5]*30^3

#and at 50
coef(m)[2] + 2*coef(m)[3]*50 + 3*coef(m)[4]*50^2 + 4*coef(m)[5]*50^3

#Plotting this relationship:
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(span=.1)
p

#Plotting this relationship:
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 1), se = FALSE)
p

#Plotting this relationship:
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 4), se = FALSE)
p

#Plotting this relationship:
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 2), se = FALSE)
p

library(AER)
data("CASchools")
#Creating Variables
CASchools$STR <- CASchools$students/CASchools$teachers

CASchools$score <- (CASchools$read + CASchools$math)/2
dat <- CASchools

m <- lm(score ~ income, data=dat)

plot(dat$income, dat$score)
abline(m, col="firebrick", lwd=2)

summary(m)

dat$k8 <- NA
dat$k8[dat$grades=="KK-08"] <- 1
dat$k8[dat$grades=="KK-06"] <- 0
table(dat$k8)

plot(dat$income[dat$k8==1], dat$score[dat$k8==1], col="forestgreen", pch=16)
points(dat$income[dat$k8==0], dat$score[dat$k8==0], col="goldenrod", pch=16)
abline(m, col="black", lwd=2)

m <- lm(score ~ income + k8, data=dat)
summary(m)

income <- seq(0,70, 1)
pred1 <- coef(m)["(Intercept)"] + coef(m)["income"]*income + coef(m)["k8"]
pred0 <- coef(m)["(Intercept)"] + coef(m)["income"]*income

plot(dat$income[dat$k8==1], dat$score[dat$k8==1], col="forestgreen", pch=16)
points(dat$income[dat$k8==0], dat$score[dat$k8==0], col="goldenrod", pch=16)
points(income, pred1, col="forestgreen", lwd=2, type="l")
points(income, pred0, col="goldenrod", lwd=2, type="l")

abline(lm(score ~income, data=dat), lwd=2)

m <- lm(score ~ income*k8, data=dat)
summary(m)

coef(m)["income"] + coef(m)["income:k8"]*0

coef(m)["income"] + coef(m)["income:k8"]*1

income <- seq(0,70, 1)
pred1 <- coef(m)["(Intercept)"] +coef(m)["income"]*income + coef(m)["k8"] + coef(m)["k8"] + coef(m)["income:k8"]*income
pred0 <- coef(m)["(Intercept)"] + coef(m)["income"]*income


plot(dat$income[dat$k8==1], dat$score[dat$k8==1], col="forestgreen", pch=16)
points(dat$income[dat$k8==0], dat$score[dat$k8==0], col="goldenrod", pch=16)
points(income, pred1, col="forestgreen", lwd=2, type="l")
points(income, pred0, col="goldenrod", lwd=2, type="l")

library(rio)
acs <- import("https://github.com/marctrussler/IIS-Data/raw/main/ACSCountyData.csv")

acs$median.income <- acs$median.income/1000
acs$gini <- acs$gini*100

plot(acs$median.income, acs$gini, xlab="Median Income (thousands)",
     ylab="Gini Coefficient")
abline(lm(acs$gini ~ acs$median.income), col="firebrick", lwd=2)

summary(lm(gini ~ median.income, data=acs))

m1 <- lm(gini ~ median.income, data=acs)

m2 <- lm(gini ~ median.income + percent.white, data=acs)

#Use stargazer to view two models together
library(stargazer)

stargazer(m1, m2, type="text")

acs$majority.white[acs$percent.white>=50] <- 1
acs$majority.white[acs$percent.white<50] <- 0
m3 <- lm(gini ~ median.income*majority.white, data=acs)

stargazer(m1,m2,m3, type="text")

m4 <- lm(gini ~ median.income*percent.white, data=acs)

stargazer(m1,m2,m3,m4, type="text")
summary(m4)

median.income <- seq(0,150,.01)

pred0 <- coef(m4)["(Intercept)"] +
  coef(m4)["median.income"]*median.income +
  coef(m4)["percent.white"]*0 +
  coef(m4)["median.income:percent.white"]*0*median.income
pred25 <- coef(m4)["(Intercept)"] +
  coef(m4)["median.income"]*median.income +
  coef(m4)["percent.white"]*25 +
  coef(m4)["median.income:percent.white"]*25*median.income
pred50 <- coef(m4)["(Intercept)"] +
  coef(m4)["median.income"]*median.income +
  coef(m4)["percent.white"]*50 +
  coef(m4)["median.income:percent.white"]*50*median.income
pred75 <- coef(m4)["(Intercept)"] +
  coef(m4)["median.income"]*median.income +
  coef(m4)["percent.white"]*75 +
  coef(m4)["median.income:percent.white"]*75*median.income
pred100 <- coef(m4)["(Intercept)"] +
  coef(m4)["median.income"]*median.income +
  coef(m4)["percent.white"]*100 +
  coef(m4)["median.income:percent.white"]*100*median.income

library(RColorBrewer)

cols <- brewer.pal(7, "Set2")

plot(acs$median.income, acs$gini, col="gray80")
points(median.income, pred0,lwd=2, type="l", col=cols[1])
points(median.income, pred25,lwd=2, type="l", col=cols[2])
points(median.income, pred50,lwd=2, type="l",  col=cols[3])
points(median.income, pred75,lwd=2, type="l",  col=cols[4])
points(median.income, pred100,lwd=2, type="l",  col=cols[5])
legend("topleft", c("0% White",
                    "25% White",
                    "50% White",
                    "75% White",
                    "100% White"),
       lty=rep(1,6), col=cols)

percent.white <- seq(0,100,.01)

margin <- coef(m4)["median.income"] + coef(m4)["median.income:percent.white"]*percent.white

plot(percent.white, margin, type="l", ylim=c(-0.15,0))
abline(h=0, lty=2)
rug(acs$percent.white)

median.income <- seq(0,140,1)

margin <- coef(m4)["percent.white"] + coef(m4)["median.income:percent.white"]*median.income

plot(median.income, margin, type="l", ylim=c(-0.15,0))
abline(h=0, lty=2)
rug(acs$median.income)

library(margins)
cplot(m4, x="median.income",dx="percent.white", what="effect")
cplot(m4, dx="median.income",x="percent.white", what="effect")

margins(m4, at=list(percent.white=c(0,25,50,75,100)))
margins(m4, at=list(median.income=c(0,50,150)))

m5 <- lm(gini ~ median.income*percent.white + population.density +
           percent.adult.poverty + percent.car.commute, data=acs)

summary(m5)

percent.white <- seq(0,100,.01)

margin <- coef(m5)["median.income"] + coef(m5)["median.income:percent.white"]*percent.white

plot(percent.white, margin, type="l", ylim=c(-0.15,0.15))
abline(h=0, lty=2)
rug(acs$percent.white)

cplot(m5, dx="median.income",x="percent.white", what="effect", ylim=c(-.15,.15))
abline(h=0, lty=2)

library(stargazer)
m <- lm(average.commute.time ~ percent.transit.commute, data=acs)
stargazer(m, type="text")

m2 <- lm(average.commute.time ~ percent.transit.commute + population.density, data=acs)
stargazer(m,m2, type="text")


m3 <- lm(average.commute.time ~ percent.transit.commute*population.density, data=acs)

stargazer(m,m2,m3, type="text")

summary(acs$population.density)

population.density <- seq(0,73000,1)
marginal.effect <- coef(m3)["percent.transit.commute"] + population.density*coef(m3)["percent.transit.commute:population.density"]

plot(population.density, marginal.effect, type="l", xlab="Population Density", ylab="Effect of % Transit Commute on Commute Time")
abline(h=0, lty=2)
rug(acs$population.density)