library(MASS)
set.seed(19104)
sigma<-rbind(c(1,.2), c(.2,1))
sigma
mu <- c(0,0)
d <- as.data.frame(mvrnorm(n=100, mu=mu, Sigma=sigma))
names(d) <- c("x","y")
plot(d$x, d$y)
abline(lm(d$y ~ d$x), lwd=2, col="firebrick")
summary(lm(d$y ~ d$x))


par(mfrow=c(2,3))
for(i in 1:6){
d <- as.data.frame(mvrnorm(n=100, mu=mu, Sigma=sigma))
names(d) <- c("x","y")
m <- lm(d$y ~ d$x)

plot(d$x, d$y, main=paste("Slope = ", round(m$coefficients["d$x"],3)))
abline(lm(d$y ~ d$x), lwd=2, col="firebrick")
}

par(mfrow=c(1,2))
plot(d$x, d$y, type="n", xlim=c(-3,3), ylim=c(-3,3))
betas <- rep(NA, 1000)

for(i in 1:1000){
d <- as.data.frame(mvrnorm(n=100, mu=mu, Sigma=sigma))
names(d) <- c("x","y")
m <- lm(d$y ~ d$x)
abline(lm(d$y ~ d$x),  lwd=2, col="gray80")
betas[i] <- m$coefficients["d$x"]
}
abline(h=0, lty=2)
abline(v=0, lty=2)

plot(density(betas), main="Sampling Distribution of Slope coefficient")
legend("topright", paste("se=", round(sd(betas), 3)))

eval <- seq(-.2, .6, .001)
plot(density(betas), main="Sampling Distribution of Slope coefficient")
points(eval, dnorm(eval, mean=mean(betas), sd=sd(betas)), type="l", col="firebrick")

par(mfrow=c(1,2))
plot(d$x, d$y, type="n", xlim=c(-3,3), ylim=c(-3,3))
betas <- rep(NA, 1000)

for(i in 1:1000){
d <- as.data.frame(mvrnorm(n=500, mu=mu, Sigma=sigma))
names(d) <- c("x","y")
m <- lm(d$y ~ d$x)
abline(lm(d$y ~ d$x),  lwd=2, col="gray80")
betas[i] <- m$coefficients["d$x"]
}
abline(h=0, lty=2)
abline(v=0, lty=2)

plot(density(betas), main="Sampling Distribution of Slope coefficient", xlim=c(-0.2,0.6))
legend("topright", paste("se=", round(sd(betas), 3)))


sigma<-rbind(c(1,.2), c(.2,.5))
sigma
par(mfrow=c(1,2))
plot(d$x, d$y, type="n", xlim=c(-3,3), ylim=c(-3,3))
betas <- rep(NA, 1000)

for(i in 1:1000){
d <- as.data.frame(mvrnorm(n=100, mu=mu, Sigma=sigma))
names(d) <- c("x","y")
m <- lm(d$y ~ d$x)
abline(lm(d$y ~ d$x),  lwd=2, col="gray80")
betas[i] <- m$coefficients["d$x"]
}
abline(h=0, lty=2)
abline(v=0, lty=2)

plot(density(betas), main="Sampling Distribution of Slope coefficient")
legend("topright", paste("se=", round(sd(betas), 3)))

sigma<-rbind(c(.5,.2), c(.2,1))

par(mfrow=c(1,2))
plot(d$x, d$y, type="n", xlim=c(-3,3), ylim=c(-3,3))
betas <- rep(NA, 1000)

for(i in 1:1000){
d <- as.data.frame(mvrnorm(n=100, mu=mu, Sigma=sigma))
names(d) <- c("x","y")
m <- lm(d$y ~ d$x)
abline(lm(d$y ~ d$x),  lwd=2, col="gray80")
betas[i] <- m$coefficients["d$x"]
}
abline(h=0, lty=2)
abline(v=0, lty=2)

plot(density(betas), main="Sampling Distribution of Slope coefficient")
legend("topright", paste("se=", round(sd(betas), 3)))


betas <- rep(NA, 1000)

par(mfrow=c(1,3))
for(i in 1:3){
sigma<-rbind(c(5,.2), c(.2,1))
d <- as.data.frame(mvrnorm(n=100, mu=mu, Sigma=sigma))
names(d) <- c("x","y")
plot(d$x, d$y, type="p", xlim=c(-7,7), ylim=c(-3,3))
abline(lm(d$y ~ d$x),  lwd=2, col="firebrick")
}



betas <- rep(NA, 1000)

par(mfrow=c(1,3))
for(i in 1:3){
sigma<-rbind(c(.5,.2), c(.2,1))
d <- as.data.frame(mvrnorm(n=100, mu=mu, Sigma=sigma))
names(d) <- c("x","y")
plot(d$x, d$y, type="p", xlim=c(-7,7), ylim=c(-3,3))
abline(lm(d$y ~ d$x),  lwd=2, col="firebrick")
}


summary(m)
sqrt(mean(m$residuals^2)/sum((d$x-mean(d$x))^2))

eval <- seq(-3,3,.001)
plot(eval, dt(eval, df=98), type="l", main="t df=98")

t.stat <- m$coefficients["d$x"]/.12728
eval <- seq(-3,3,.001)
plot(eval, dt(eval, df=98), type="l", main="t df=98")
abline(v=c(-t.stat, t.stat), lty=2)

pt(-t.stat, df=98)*2

#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")

#install.packages("rgl")
#library(rgl)
#plot3d(dat$STR, dat$english,dat$score, xlab="Student:Teacher Ratio", ylab="% ESL Students", zlab="Test Scores")

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)