dat1 <-read.table("lecture4ex1.csv", sep=",", header=TRUE)
library(lme4)
out1r <- lmer(consume ~ 1 + I(length - 15) + (1|groupid), data=dat1, REML=FALSE)
out1f <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1|groupid), data=dat1, REML=FALSE)
anova(out1r, out1f)


library(lme4)
out1m0 <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1|groupid), data=dat1, 
               REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))

out1m0a <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1 + I(length - 15)|groupid), 
                data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1m0b <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1 + goldcolor|groupid), 
                data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))

anova(out1m0, out1m0a)
anova(out1m0, out1m0b)
summary(out1m0a)

ranef1m0a <- ranef(out1m0a)
plot(ranef1m0a)



dat2 <-read.table("lecture4ex2.csv", sep=",", header=TRUE)
out2m0 <- lmer(score ~ 1 + eesex + I(scale(iq)) + (1|erid), data=dat2, 
               REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out2m0a <- lmer(score ~ 1 + eesex + I(scale(iq)) + (1 + eesex|erid), data=dat2, 
                REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out2m0b <- lmer(score ~ 1 + eesex + I(scale(iq)) + (1 + I(scale(iq))|erid), data=dat2, 
                REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m0, out2m0a)
anova(out2m0, out2m0b)

out2m0a <- lmer(score ~ 1 + eesex + I(iq - 100) + (1 + eesex|erid), data=dat2, REML=FALSE)
summary(out2m0a)

ranef2m0a <- ranef(out2m0a)
plot(ranef2m0a)

out2m0c <- lmer(score ~ 1 + eesex + I(iq - 100) + ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
summary(out2m0c)



dat3 <-read.table("lecture4ex3.csv", sep=",", header=TRUE)
out3m0 <- lmer(sat ~ 1 + I(age - 40) + female + (1|tableid), data=dat3, 
               REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out3m0a <- lmer(sat ~ 1 + I(age - 40) + female + (1 + I(age - 40)|tableid), data=dat3, 
                REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out3m0b <- lmer(sat ~ 1 + I(age - 40) + female + (1 + female|tableid), data=dat3, 
                REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3m0, out3m0a)
anova(out3m0, out3m0b)



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

out2m1 <- lmer(score ~ 1 + eesex + I(iq - 100) + (1 + eesex|erid), data=dat2, REML=FALSE)
out2m1a <- lmer(score ~ 1 + eesex + I(iq - 100) + ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m1, out2m1a)


out2m1b <- lmer(score ~ 1 + eesex + I(iq - 100) + ersex + eesex:ersex 
                + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m1a, out2m1b)
summary(out2m1b)

library(ggplot2)
eesex <- c("Male", "Male", "Female", "Female")
ersex <- c("Male", "Female", "Male", "Female")
ratings <- c(75.54, 75.55, 75.13, 76.21)
mydat <- data.frame(eesex, ersex, sat)

ggplot(mydat, aes(factor(eesex), ratings, fill = ersex)) + 
  geom_bar(stat="identity", position = "dodge") + 
  labs(x = "Interviewee", y = "Interview Ratings", fill = "Interviewer")


sumout2 <- summary(out2m1b)
coef(sumout2)

psi <- coef(sumout2)[,"Estimate"]
cvec1 <- c(0, 1, 0, 0, 1) 
t(cvec1) %*% psi


vcov(sumout2)

vpsi <- vcov(sumout2)
sqrt(t(cvec1) %*% vpsi %*% cvec1)

psi <- coef(sumout2)[,"Estimate"]
vpsi <- vcov(sumout2)
cvec1 <- c(1, 0, 0, 0, 0) 
cvec2 <- c(1, 0, 0, 1, 0) 
cvec3 <- c(0, 1, 0, 0, 0) 
cvec4 <- c(0, 1, 0, 0, 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)

cvec1 <- c(1, 0, 0, 0, 0) 
cvec2 <- c(1, 1, 0, 0, 0) 
cvec3 <- c(0, 0, 0, 1, 0) 
cvec4 <- c(0, 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)

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

out1m1 <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1 + I(length - 15)|groupid), 
                data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))

out1m1intcept <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) 
                      + (1 + I(length - 15)|groupid), data=dat1, REML=FALSE, 
                      control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m1, out1m1intcept)
summary(out1m1intcept)


out1m1slope <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) 
                    + I(length - 15):I(volume - 1) + I(length - 15):I(nfish - 8) 
                    + (1 + I(length - 15)|groupid), 
                    data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))

anova(out1m1intcept, out1m1slope)
summary(out1m1slope)

out1m1trimmed <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) 
                      + I(volume - 1):I(length - 15) + (1 + I(length - 15)|groupid), 
                      data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))

anova(out1m1intcept, out1m1trimmed)

summary(out1m1trimmed)

quantile(dat1$length)
quantile(dat1$volume)

mean(dat1$length) - sd(dat1$length)
mean(dat1$length)
mean(dat1$length) + sd(dat1$length)

mean(dat1$volume) - sd(dat1$volume)
mean(dat1$volume)
mean(dat1$volume) + sd(dat1$volume)

plot(dat1$volume, dat1$consume, type="n", xlab="Volume", ylab="Daily Consumption")
abline(a = 38.21 - 13.72*1, b = 13.78, col="red")
abline(a = 42.81 - 14.73*1, b = 14.73, col="green")
abline(a = 47.41 - 15.68*1, b = 15.68, col="blue")
legend(0.1, 60, c("9 cm", "14 cm", "19 cm"), col=c("red", "green", "blue"), 
       lty=1, title="Fish Length")


plot(dat1$length, dat1$consume, type="n", xlab="Fish Length", ylab="Daily Consumption")
abline(a = 36.27 - 0.825*15, b = 0.825, col="red")
abline(a = 40.75 - 0.882*15, b = 0.882, col="green")
abline(a = 45.22 - 0.939*15, b = 0.939, col="blue")
legend(4, 60, c("0.5 cubic m", "0.8 cubic m", "1.1 cubic m"), 
       col=c("red", "green", "blue"), lty=1, title="Tank Volume")

sumout1 <- summary(out1m1trimmed)
coef(sumout1)

psi <- coef(sumout1)[,"Estimate"]
cvec1 <- c(1, -6, 0, 0, 0, 0) 
t(cvec1) %*% psi


vcov(sumout1)

vpsi <- vcov(sumout1)
sqrt(t(cvec1) %*% vpsi %*% cvec1)


cvec1 <- c(1, -6, 0, 0, 0, 0) 
cvec2 <- c(1, -1, 0, 0, 0, 0) 
cvec3 <- c(1, 4, 0, 0, 0, 0) 
cvec4 <- c(0, 0, 0, 1, 0, -6) 
cvec5 <- c(0, 0, 0, 1, 0, -1) 
cvec6 <- c(0, 0, 0, 1, 0, 4) 
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.5, 0, 0) 
cvec2 <- c(1, 0, 0, -0.2, 0, 0) 
cvec3 <- c(1, 0, 0, 0.1, 0, 0) 
cvec4 <- c(0, 1, 0, 0, 0, -0.5) 
cvec5 <- c(0, 1, 0, 0, 0, -0.2) 
cvec6 <- c(0, 1, 0, 0, 0, 0.1) 
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 <-read.table("lecture4ex3.csv", sep=",", header=TRUE)
out3m1 <- lmer(sat ~ 1 + I(age - 40)*female + I(numperson-4) + outdoor 
               + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)

quantile(dat3$age, c(0.25, 0.50, 0.75))

plot(dat3$age, dat3$sat, type="n", xlab="Age", ylab="Consumer Satisfaction")
abline(a = 65.02 - 0.24*40, b = 0.24, col="red")
abline(a = 63.04 - 0.28*40, b = 0.28, col="blue")
legend(0.1, 85, c("Male", "Female"), col=c("red", "blue"), 
       lty=1, title="Sex")

plot(dat3$age, dat3$sat, type="n", xlab="Age", ylab="Consumer Satisfaction")
abline(a = 65.02 - 0.24*40, b = 0.24, col="red")
abline(a = 63.04 - 0.28*40, b = 0.28, col="blue")
legend(0.1, 85, c("Male", "Female"), col=c("red", "blue"), 
       lty=1, title="Sex")

sumout3 <- summary(out3m1)
coef(sumout3)

psi <- coef(sumout3)[,"Estimate"]
vpsi <- vcov(sumout3)
cvec1 <- c(1, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 0, 1, 0, 0, 0) 
cvec3 <- c(0, 1, 0, 0, 0, 0) 
cvec4 <- c(0, 1, 0, 0, 0, 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)

cvec1 <- c(1, -15, 0, 0, 0, 0) 
cvec2 <- c(1, -5, 0, 0, 0, 0) 
cvec3 <- c(1, 5, 0, 0, 0, 0) 
cvec4 <- c(0, 0, 1, 0, 0, -15) 
cvec5 <- c(0, 0, 1, 0, 0, -5) 
cvec6 <- c(0, 0, 1, 0, 0, 5) 
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)