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)

small <- acs[c("percent.transit.commute","population.density","average.commute.time")]
b1 <- rep(NA, 1000)
b3 <- rep(NA, 1000)
for(i in 1:1000){
  bs.data <- small[sample(1:nrow(small), nrow(small),replace=T),]
  m.bs <- lm(average.commute.time ~ percent.transit.commute*population.density, data=bs.data)
  b1[i] <- coef(m.bs)["percent.transit.commute"]
  b3[i] <- coef(m.bs)["percent.transit.commute:population.density"]
}
cor(b1,b3)

x <- rnorm(100)
y <- 2*x*rnorm(100)

cov(x,2*y) ==2*cov(x,y)

vcov(m3)


var.margin <- vcov(m3)["percent.transit.commute","percent.transit.commute"] + 
              population.density^2*vcov(m3)["percent.transit.commute:population.density","percent.transit.commute:population.density"] + 
              population.density*2*vcov(m3)["percent.transit.commute:population.density","percent.transit.commute"] 
            

se.margin <- sqrt(var.margin)

upper.bounds <- marginal.effect +se.margin*1.96
lower.bounds <- marginal.effect - se.margin*1.96


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)
points(population.density, upper.bounds, type="l", lty=2)
points(population.density, lower.bounds, type="l", lty=2)
rug(acs$population.density)

library(margins)
cplot(m3, x="population.density",dx="percent.transit.commute", what="effect")
abline(h=0)

plot(acs$median.rent, acs$gini)
abline(lm(gini ~ median.rent, data=acs), col="firebrick", lwd=2)


m1 <- lm(gini ~ median.rent, data=acs)
m2 <- lm(gini ~ median.rent + factor(census.region), data=acs)

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

#Let the relationship change in each region
m3 <- lm(gini ~ median.rent*factor(census.region), data=acs)
stargazer(m1, m2,m3, type="text")

median.rent <- seq(200,2500,1)

pred.midwest <- coef(m3)["(Intercept)"] +
                coef(m3)["median.rent"]*median.rent
pred.northeast <- coef(m3)["(Intercept)"] +
                  coef(m3)["factor(census.region)northeast"] +
                  coef(m3)["median.rent"]*median.rent +
                  coef(m3)["median.rent:factor(census.region)northeast"]*median.rent
pred.south <- coef(m3)["(Intercept)"] +
             coef(m3)["factor(census.region)south"] +
             coef(m3)["median.rent"]*median.rent +
             coef(m3)["median.rent:factor(census.region)south"]*median.rent
pred.west <- coef(m3)["(Intercept)"] +
            coef(m3)["factor(census.region)west"] +
            coef(m3)["median.rent"]*median.rent +
            coef(m3)["median.rent:factor(census.region)west"]*median.rent

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

plot(acs$median.rent, acs$gini, col="gray80")
points(median.rent, pred.midwest, lwd=2, col=cols[1], type="l")
points(median.rent, pred.northeast, lwd=2, col=cols[2], type="l")
points(median.rent, pred.south, lwd=2, col=cols[3], type="l")
points(median.rent, pred.west, lwd=2, col=cols[4], type="l")
legend("topright", c("Midwest","Northeast","South","West"), lty=1, col=cols)

x <- runif(100,0,100)
y <- 3 + 5*x + rnorm(100,0,40)

#Scatterplot
plot(x,y)

#Regression
m <- lm(y~x)

#Predicted values
pred <- coef(m)["(Intercept)"] + coef(m)["x"]*x

#Draw lines
plot(x,y)
points(x, pred, type="l", col="firebrick")

coef(m)["(Intercept)"] + coef(m)["x"]*200

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

#Some cleaning because I copied data from wikipedia
names(mwr) <- c("time","date")
mwr <- mwr[-5,]
mwr$seconds <- seconds(hms(mwr$time))
mwr$date <- dmy(mwr$date)
mwr$date[1:30] <- mwr$date[1:30]-years(100)
mwr$year <- year(mwr$date)
mwr$seconds <- as.numeric(mwr$seconds)
mwr$year <- as.numeric(mwr$year)

#How many seconds is 2 hours?
2*3600
#What does data look like?
head(mwr)

#Plotting this relationship:
library(ggplot2)
p <- ggplot(mwr, aes(x=year, y=seconds)) + geom_point() +
  labs(x="Date", y = "Time") +
  geom_smooth(method = lm, formula = y ~ poly(x, 1), se = FALSE)
p

m <- lm(seconds ~ year, data=mwr)
summary(m)

future.years <- seq(2023,2050,1)

prediction <- coef(m)["(Intercept)"] + coef(m)["year"]*future.years
cbind(future.years, prediction/3600)

var <- vcov(m)["(Intercept)","(Intercept)"] + future.years^2 * vcov(m)["year","year"] + 2*future.years * vcov(m)["year","(Intercept)"]
se.pred <- sqrt(var)
upper <- (prediction + se.pred*1.96)/3600
lower <- (prediction - se.pred*1.96)/3600

plot(future.years, prediction/3600, type="l", ylim=c())
points(future.years, lower, type="l", lty=2)
points(future.years, upper, type="l", lty=2)

mwr$prediction <-  predict(m)
head(mwr)

prediction2 <- predict(m, newdata=data.frame(year=future.years))
cbind(prediction, prediction2)

prediction2 <- predict(m, newdata=data.frame(year=future.years), interval="confidence")
prediction2


#Plotting this relationship:
library(ggplot2)
p <- ggplot(mwr, aes(x=year, y=seconds)) + geom_point() +
  labs(x="Date", y = "Time") +
  geom_smooth(method = lm, formula = y ~ poly(x, 3), se = FALSE)
p


m2 <- lm(seconds ~ poly(year, 3), data=mwr)
summary(m2)

prediction.poly <- predict(m2, newdata=data.frame(year=future.years), interval="confidence")
plot(future.years, prediction.poly[,1], type="l", ylim=c(6000,8000))
points(future.years, prediction.poly[,2], type="l")
points(future.years, prediction.poly[,3], type="l")
abline(h=7200, lty=3)

prediction.poly <- predict(m2, newdata=data.frame(year=future.years), interval="confidence")
plot(future.years, prediction.poly[,1], type="l", ylim=c(6000,11000), xlim=c(1900,2050))
points(future.years, prediction.poly[,2], type="l")
points(future.years, prediction.poly[,3], type="l")
points(mwr$year, mwr$seconds, pch=16)
abline(h=7200, lty=3)

head(mwr)
prediction.poly <- predict(m2, newdata=data.frame(year=future.years), interval="confidence")
plot(future.years, prediction.poly[,1], type="l", ylim=c(6000,11000), xlim=c(2000,2050))
points(future.years, prediction.poly[,2], type="l")
points(future.years, prediction.poly[,3], type="l")
points(mwr$year, mwr$seconds, pch=16)
points(2023,7235, pch=16, col="firebrick" )
abline(h=7200, lty=3)

library(rio)

cces <- import("https://github.com/marctrussler/IIS-Data/raw/main/CCES.csv")
head(cces)
table(cces$year)


table(cces$vote.intent)
cces$vote.intent <- factor(cces$vote.intent)
cces$vote.intent <- relevel(cces$vote.intent, ref="No")

cces.16 <- cces[cces$year==2016,]
m <- lm(vv.turnout ~ vote.intent, data=cces.16)
summary(m)

predict(m, newdata= data.frame(vote.intent=c("No", "Maybe", "Yes", "Already voted")))

cces.training <- cces[cces$year %in% c(2016,2018),]
cces.test <- cces[cces$year==2020,]
m1 <- lm(vv.turnout ~ vote.intent, data=cces.training)
m2 <- lm(vv.turnout ~  vote.intent +prim.turnout + educ:white, data=cces.training)
m3 <- lm(vv.turnout ~ female + vote.intent +prim.turnout + educ:white:agegrp + pid3, data=cces.training)

summary(m1)
predict(m1, newdata= data.frame(vote.intent=c("No", "Maybe", "Yes", "Already voted")))

cces.test$predict1 <- predict(m1, newdata=cces.test)
head(cces.test)
table(cces.test$predict1)

#Classification Error
mean((cces.test$predict1>.5 & cces.test$vv.turnout==0) | (cces.test$predict1<.5 & cces.test$vv.turnout==1),na.rm=T)

cces.test$predict2 <- predict(m2, newdata=cces.test)
cces.test$predict3 <- predict(m3, newdata=cces.test)

#Classification Error Model 2
mean((cces.test$predict2>.5 & cces.test$vv.turnout==0) | (cces.test$predict2<.5 & cces.test$vv.turnout==1),na.rm=T)

#Classification Error Model 3
mean((cces.test$predict3>.5 & cces.test$vv.turnout==0) | (cces.test$predict3<.5 & cces.test$vv.turnout==1),na.rm=T)



#Stupid model
m4 <- lm(vv.turnout ~ pid3:vote.intent:prim.turnout + female:educ:white:agegrp, data=cces.training)
cces.test$predict4<- predict(m4, newdata=cces.test)

#Classification Error Model 4
mean((cces.test$predict4>.5 & cces.test$vv.turnout==0) | (cces.test$predict4<.5 & cces.test$vv.turnout==1),na.rm=T)