dat1 <-read.table("lecture4ex1.csv", sep=",", header=TRUE)
library(psych)
describe(dat1)
aggregate(consume ~ groupid, data=dat1, FUN=mean)
aggregate(consume ~ groupid, data=dat1, FUN=range)
stripchart(consume ~ groupid, vertical = TRUE, data = dat1[dat1$groupid < 10,])

library(beeswarm)
dat1_1 <- dat1[dat1$groupid < 10,]
beeswarm(consume ~ groupid, data=dat1_1, col=rainbow(9))

library(lme4)
out1m0 <- lmer(consume ~ 1 + (1|groupid), data=dat1, REML=FALSE)
summary(out1m0)


dat2 <-read.table("lecture4ex2.csv", sep=",", header=TRUE)
describe(dat2)

out2m0 <- lmer(score ~ 1 + (1|erid), data=dat2, REML=FALSE)
summary(out2m0)


dat3 <-read.table("lecture4ex3.csv", sep=",", header=TRUE)
describe(dat3)

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

out1m1 <- lmer(consume ~ nfish + volume + (1|groupid), data=dat1, REML=FALSE)
summary(out1m1)

out1m1_1 <- lmer(consume ~ I(nfish - mean(nfish)) + I(volume - mean(volume)) + (1|groupid), data=dat1, REML=FALSE)
summary(out1m1_1)

summary1m1_1 <- summary(out1m1_1)
coef1m1_1 <- coef(summary1m1_1)
coef1m1_1
tvalue1m1_1 <- coef1m1_1[,"t value"]
tvalue1m1_1
pnorm(abs(tvalue1m1_1), lower.tail=FALSE) * 2


out2m1 <- lmer(score ~ 1 + ersex + (1|erid), data=dat2, REML=FALSE)
summary(out2m1)


out3m1 <- lmer(sat ~ 1 + I(numperson - 4) + outdoor + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)


out1m2 <- lmer(consume ~ I(length - 15) + goldcolor + (1|groupid), data=dat1, REML=FALSE)
summary(out1m2)

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

out3m2 <- lmer(sat ~ 1 + I(age - 40) + female + (1|tableid), data=dat3, REML=FALSE)
summary(out3m2)

out1m3 <- lmer(consume ~ I(nfish - mean(nfish)) + I(volume - mean(volume)) + I(length - 15) + goldcolor + (1|groupid), data=dat1, REML=FALSE)
summary(out1m3)

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

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

library(ggplot2)

dat3_1 <- dat3[dat3$tableid%%5 == 0,]
out <- ggplot(dat3_1, aes(x=age, y=sat, group=tableid)) + geom_smooth(method=lm, se=FALSE)
btwslope <- 0.267-0.318
out + geom_abline(intercept=64.2875+(btwslope*40), slope=btwslope, color="red")

dat3$avefemale <- ave(dat3$female, dat3$tableid)
out3m4 <- lmer(sat ~ female + I(avefemale - 0.5) + (1|tableid), data=dat3, REML=FALSE)
summary(out3m4)

dat3_1 <- dat3[dat3$tableid%%5 == 0,]
out <- ggplot(dat3_1, aes(x=female, y=sat, group=tableid)) + geom_smooth(method=lm, se=FALSE)
btwslope <- -1.5293-1.1253
out + geom_abline(intercept=65.0817+(btwslope*0.5), slope=btwslope, color="red")


out3m5 <- lmer(sat ~ female + I(age - 40) + I(avefemale - 0.5) + I(aveage - 40) + I(numperson - 4) + outdoor + (1|tableid), data=dat3, REML=FALSE)
summary(out3m5)