#install.packages("AER")
library(AER)

data("CASchools")

head(CASchools)

CASchools$STR <- CASchools$students/CASchools$teachers
summary(CASchools$STR)

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

dat <- CASchools

plot(dat$STR, dat$score, xlab="Student:Teacher Ratio", ylab="Test Score", pch=16)

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

plot(dat$STR, dat$score, xlab="Student:Teacher Ratio", ylab="Test Score", pch=16)
abline(lm(dat$score ~ dat$STR), lwd=2, col="firebrick")

dat$high.esl <-NA 
dat$high.esl[dat$english>=median(dat$english)]<- 1
dat$high.esl[dat$english<median(dat$english)]<- 0

plot(dat$STR, dat$score, xlab="Student:Teacher Ratio", ylab="Test Score", pch=16)
abline(lm(dat$score ~ dat$STR), lwd=2, col="firebrick")

plot(dat$STR, dat$score, xlab="Student:Teacher Ratio", ylab="Test Score", pch=16, type="n")
points(dat$STR[dat$high.esl==0], dat$score[dat$high.esl==0], pch=16, col="darkblue")
points(dat$STR[dat$high.esl==1], dat$score[dat$high.esl==1], pch=16, col="firebrick")
legend("topright", c("Low ESL", "High ESL"), pch=c(16,16), col=c("darkblue", "firebrick"))

plot(dat$STR, dat$score, xlab="Student:Teacher Ratio", ylab="Test Score", pch=16, type="n")
points(dat$STR[dat$high.esl==0], dat$score[dat$high.esl==0], pch=16, col="darkblue")
points(dat$STR[dat$high.esl==1], dat$score[dat$high.esl==1], pch=16, col="firebrick")
abline(a=691.32, b=-1.3963, col="darkblue", lwd=3)
abline(a=691.32-19.49, b=-1.3963, col="firebrick", lwd=3)
legend("topright", c("Low ESL", "High ESL"), pch=c(16,16), col=c("darkblue", "firebrick"))

plot(dat$STR, dat$score, xlab="Student:Teacher Ratio", ylab="Test Score", pch=16, type="n", xlim=c())
points(dat$STR[dat$high.esl==0], dat$score[dat$high.esl==0], pch=16, col="darkblue")
points(dat$STR[dat$high.esl==1], dat$score[dat$high.esl==1], pch=16, col="firebrick")
abline(a=691.32, b=-1.3963, col="darkblue", lwd=3)
abline(a=691.32-19.49, b=-1.3963, col="firebrick", lwd=3)
abline(lm(dat$score ~ dat$STR), lwd=3)
legend("topright", c("Low ESL", "High ESL"), pch=c(16,16), col=c("darkblue", "firebrick"))
abline(v=0, lty=2)

m <- lm(score ~ STR + high.esl, data=dat)
summary(m)

plot(dat$STR, dat$score, xlab="Student:Teacher Ratio", ylab="Test Score", pch=16)
abline(lm(dat$score ~ dat$STR), lwd=2, col="firebrick")

plot(dat$english, dat$score, xlab="% ESL Students", ylab="Test Score", pch=16)
abline(lm(dat$score ~ dat$english), lwd=2, col="firebrick")

library(plotly)
fig <- plot_ly(dat, x = ~STR, y = ~english, z = ~score)
fig <- fig %>% add_markers()
fig

m <- lm(score ~ STR + english, data=dat)
summary(m)

m <- lm(score ~ STR + english + lunch + income, data=dat)
summary(m)

m <- lm(score ~english + income + STR +  lunch, data=dat)
summary(m)

summary(dat$income)

dat$income2 <- dat$income*1000

m <- lm(score ~english + income2 + STR +  lunch, data=dat)
summary(m)

m <- lm(score ~english + income + STR +  lunch, data=dat)
summary(m)

sd(dat$STR)

-.56*1.89

#We can do this a bit more programatically to determine the effect of a one sd shift for each of our variables

vars <- c("english","income","STR","lunch")

#Get sd of each
sds <- rep(NA, length(vars))
for(i in 1:length(vars)){
  sds[i] <- sd( dat[,vars[i]] )
}

#The coef() command will pull the coefficients from a regression
c <- coef(m)[vars]

#Standardized impacts:
summary(m)
c*sds


#Regression with income in real dollars
m <- lm(score ~english + income2 + STR +  lunch, data=dat)
summary(m)

#We can do this a bit more programatically to determine the effect of a one sd shift for each of our variables

vars <- c("english","income2","STR","lunch")

#Get sd of each
sds <- rep(NA, length(vars))
for(i in 1:length(vars)){
  sds[i] <- sd( dat[,vars[i]] )
}

#The coef() command will pull the coefficients from a regression
c <- coef(m)[vars]

#Standardized impacts:
summary(m)
c*sds


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

cor(acs$population.density, acs$percent.transit.commute, use="pairwise.complete")

acs$density[acs$population.density<50] <- 1
acs$density[acs$population.density>=50 & acs$population.density<1000] <- 2
acs$density[acs$population.density>=1000] <- 3

table(acs$density)

plot(acs$density, acs$percent.transit.commute)

boxplot(acs$percent.transit.commute ~ acs$density)

boxplot(acs$percent.transit.commute ~ acs$density, outline=F)

m <- lm(percent.transit.commute ~ density, data=acs)
summary(m)

plot(acs$density, acs$percent.transit.commute, xlim=c(0,4))
abline(m, col="firebrick", lwd=2)

acs$rural <- acs$density==1
acs$suburban <- acs$density==2
acs$urban <- acs$density==3

table(acs$rural, acs$density)
table(acs$suburban, acs$density)
table(acs$urban, acs$density)


m <- lm(percent.transit.commute ~ rural + suburban + urban, data=acs)
summary(m)

m <- lm(percent.transit.commute ~ rural + suburban, data=acs)
summary(m)

mean(acs$percent.transit.commute[acs$urban],na.rm=T)

mean(acs$percent.transit.commute[acs$rural],na.rm=T) - mean(acs$percent.transit.commute[acs$urban],na.rm=T) 

m <- lm(percent.transit.commute ~ rural + urban, data=acs)
summary(m)

m1 <- lm(median.income ~ percent.transit.commute, data=acs)
summary(m1)
plot(acs$percent.transit.commute, acs$median.income)
abline(lm(median.income ~ percent.transit.commute, data=acs), lwd=2, col="firebrick")

table(acs$census.region)

acs$midwest[acs$census.region=="midwest"]<- 1
acs$midwest[acs$census.region!="midwest"]<- 0
acs$northeast[acs$census.region=="northeast"]<- 1
acs$northeast[acs$census.region!="northeast"]<- 0
acs$south[acs$census.region=="south"]<- 1
acs$south[acs$census.region!="south"]<- 0
acs$west[acs$census.region=="west"]<- 1
acs$west[acs$census.region!="west"]<- 0

#Check if this worked:
table(acs$census.region, acs$midwest)
table(acs$census.region, acs$northeast)
table(acs$census.region, acs$south)
table(acs$census.region, acs$west)

m <- lm(median.income ~ midwest + northeast + south + west, data=acs)
summary(m)

m <- lm(median.income ~  northeast + south + west, data=acs)
summary(m)

mean(acs$median.income[acs$census.region=="midwest"])

coef(m)["(Intercept)"] + coef(m)["northeast"]
mean(acs$median.income[acs$census.region=="northeast"])

m <- lm(median.income ~  midwest + south + west, data=acs)
summary(m)

m1 <- lm(median.income ~ percent.transit.commute, data=acs)
summary(m1)

m2 <- lm(median.income ~ percent.transit.commute + northeast + south + west, data=acs)
summary(m2)

acs$census.region <- as.factor(acs$census.region)

m <- lm(median.income ~ percent.transit.commute + census.region, data=acs)
summary(m)

acs$census.region <- relevel(acs$census.region, ref="northeast")


m <- lm(median.income ~ percent.transit.commute + census.region, data=acs)
summary(m)


m1 <- lm(median.income ~ percent.transit.commute, data=acs)
summary(m1)

m2 <- lm(median.income ~ percent.transit.commute + northeast + south + west, data=acs)
summary(m2)

#install.packages("coefplot")
library(coefplot)

coefplot(m2)

library(stargazer)

#The "text" type will output into R in a nice looking format:
stargazer(m1, m2, type="text")


stargazer(m1, m2, type="html", style="apsr",
          dep.var.labels = "Median Income",
          covariate.labels = c("% Transit Commuters",
                               "Region: NorthEast",
                               "Region: South",
                               "Region: West"))