dat3 <-read.table("lecture4ex3.csv", sep=",", header=TRUE)
dat3$aveage <- ave(dat3$age, dat3$tableid)
dat3$diffage <- dat3$age - dat3$aveage

library(lme4)

out3m1 <- lmer(sat ~ 1 + age + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)

out3m2 <- lmer(sat ~ 1 + diffage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m2)

out3m3 <- lmer(sat ~ 1 + age + aveage + (1|tableid), 
               data=dat3, REML=FALSE)
summary(out3m3)

out3m4 <- lmer(sat ~ 1 + diffage + aveage + (1|tableid), 
               data=dat3, REML=FALSE)
summary(out3m4)

dat3$sage <- dat3$age / 10
dat3$avesage <- ave(dat3$sage, dat3$tableid)
dat3$diffsage <- dat3$sage - dat3$avesage

out3m1 <- lmer(sat ~ 1 + sage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)

out3m2 <- lmer(sat ~ 1 + diffsage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m2)

out3m3 <- lmer(sat ~ 1 + sage + avesage + (1|tableid), 
               data=dat3, REML=FALSE)
summary(out3m3)

out3m4 <- lmer(sat ~ 1 + diffsage + avesage + (1|tableid), 
               data=dat3, REML=FALSE)
summary(out3m4)

out3m5 <- lmer(sat ~ 1 + sage + avesage + sage:avesage 
               + (1 + sage|tableid), data=dat3, REML=FALSE)
summary(out3m5)

out3m6 <- lmer(sat ~ 1 + diffsage + avesage + diffsage:avesage 
               + (1 + diffsage|tableid), data=dat3, REML=FALSE)
summary(out3m6)



dat4 <-read.table("lecture6ex1.csv", sep=",", header=TRUE)
dat4$aveach <- ave(dat4$ach, dat4$schoolid)
dat4$diffach <- dat4$ach - dat4$aveach
out4m1 <- lmer(efficacy ~ 1 + diffach + I(aveach - 50) + private
               + diffach:I(aveach - 50)
               + (1 + diffach|schoolid), data=dat4, REML=FALSE)
summary(out4m1)


out4m0 <- lmer(ach ~ 1 + (1|schoolid), data=dat4, REML=FALSE)
summary(out4m0)

plot(dat4$diffach, dat4$efficacy, type="n", xlab="Math Achievement Deviation within School", ylab="Self Efficacy")
abline(a = 35.35, b = 0.01, col="red")
abline(a = 35.85, b = 0.11, col="green")
abline(a = 34.85, b = 0.21, col="blue")
legend(-40, 50, c("Low", "Medium", "High"), col=c("red", "green", "blue"), 
       lty=1, title="School Math Achievement")


plot(dat4$ach, dat4$efficacy, type="n", xlab="School Math Achievement", ylab="Self Efficacy")
abline(a = 35 - (-0.15*50), b = -0.15, col="red")
abline(a = 35.6 - (-0.05*50), b = -0.05, col="green")
abline(a = 36.2 - (0.05*50), b = 0.05, col="blue")
legend(20, 30, c("Low", "Medium", "High"), col=c("red", "green", "blue"), 
       lty=1, title="Relative Position within School")

sumout4 <- summary(out4m1)
psi <- coef(sumout4)[,1]

vpsi <- vcov(sumout4)
cvec1 <- c(1, 0, -5, 0, 0) 
cvec2 <- c(1, 0, 5, 0, 0) 
cvec3 <- c(1, 0, 15, 0, 0) 
cvec4 <- c(0, 1, 0, 0, -5) 
cvec5 <- c(0, 1, 0, 0, 5) 
cvec6 <- c(0, 1, 0, 0, 15) 
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6)
est <- t(cmat) %*% psi
se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat))
zval <- est/se
pval <- round(pnorm(-abs(zval)) * 2, 6)
cbind(est, se, zval, pval)

cvec1 <- c(1, -10, 0, 0, 0) 
cvec2 <- c(1, 0, 0, 0, 0) 
cvec3 <- c(1, 10, 0, 0, 0) 
cvec4 <- c(0, 0, 1, 0, -10) 
cvec5 <- c(0, 0, 1, 0, 0) 
cvec6 <- c(0, 0, 1, 0, 10) 
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6)
est <- t(cmat) %*% psi
se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat))
zval <- est/se
pval <- round(pnorm(-abs(zval)) * 2, 6)
cbind(est, se, zval, pval)





#######################################

dat3$avefemale <- ave(dat3$female, dat3$tableid)
out3g1 <- lmer(sat ~ 1 + female + avefemale 
               + (1|tableid), data=dat3, REML=FALSE)
out3g2 <- lmer(sat ~ 1 + female + avefemale 
               + (1 + female|tableid), data=dat3, REML=FALSE)
anova(out3g1, out3g2)

out3g3 <- lmer(sat ~ 1 + female + avefemale + female:avefemale 
                + (1 + female|tableid), data=dat3, REML=FALSE)
anova(out3g2, out3g3)

summary(out3g3)

avefemale <- seq(0, 1, 0.01)
cust <- 66.5712 + ((-3.8363 - 3.646) * avefemale) + (4.6837*avefemale*avefemale)
plot(avefemale, cust, type="l", xlim=c(0, 1), ylim=c(60, 70), 
     xlab="Female Proportion in a Table", ylab="Customer Satisfaction")


plot(dat3$female, dat3$sat, type="n", xlab="Female", ylab="Customer Satisfaction")
abline(a = 65.84, b = -2.904, col="red")
abline(a = 64.75, b = -1.5, col="green")
abline(a = 63.65, b = -0.096, col="blue")
legend(0.01, 50, c("20%", "50%", "80%"), col=c("red", "green", "blue"), 
       lty=1, title="Percent Female in the Table")


sumout3 <- summary(out3g3)
coef(sumout3)


psi <- coef(sumout3)[,"Estimate"]
vpsi <- vcov(sumout3)
cvec1 <- c(1, 0, 0.2, 0) 
cvec2 <- c(1, 0, 0.5, 0) 
cvec3 <- c(1, 0, 0.8, 0) 
cvec4 <- c(0, 1, 0, 0.2) 
cvec5 <- c(0, 1, 0, 0.5) 
cvec6 <- c(0, 1, 0, 0.8) 
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6)
est <- t(cmat) %*% psi
se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat))
zval <- est/se
pval <- round(pnorm(-abs(zval)) * 2, 6)
cbind(est, se, zval, pval)

cvec1 <- c(1, 0, 0, 0) 
cvec2 <- c(1, 1, 0, 0) 
cvec3 <- c(0, 0, 1, 0) 
cvec4 <- c(0, 0, 1, 1) 
cmat <- cbind(cvec1, cvec2, cvec3, cvec4)
est <- t(cmat) %*% psi
se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat))
zval <- est/se
pval <- round(pnorm(-abs(zval)) * 2, 6)
cbind(est, se, zval, pval)



###################

dat3$avefemale <- ave(dat3$female, dat3$tableid)
dat3$difffemale <- dat3$female - dat3$avefemale

out3f1 <- lmer(sat ~ 1 + difffemale + avefemale + (1|tableid), 
               data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
summary(out3f1)

out3f2 <- lmer(sat ~ 1 + difffemale + avefemale 
               + (1 + difffemale|tableid), data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f1, out3f2)

out3f3 <- lmer(sat ~ 1 + difffemale + avefemale + difffemale:avefemale 
               + (1 + difffemale|tableid), data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f2, out3f3)

out3f4 <- lmer(sat ~ 1 + difffemale + avefemale + I(avefemale^2)
               + (1 + difffemale|tableid), data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f2, out3f4)

out3f5 <- lmer(sat ~ 1 + difffemale + avefemale + I(avefemale^2)
               + avefemale:difffemale + I(avefemale^2):difffemale
               + (1 + difffemale|tableid), data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f4, out3f5)

summary(out3f4)


avefemale <- seq(0, 1, 0.01)
cust <- 67.6715 - (13.1943 * avefemale) + (10.2163*avefemale*avefemale)
plot(avefemale, cust, type="l", xlim=c(0, 1), ylim=c(60, 70), 
     xlab="Female Proportion in a Table", ylab="Customer Satisfaction")