dat1 <-read.table("lecture4ex1.csv", sep=",", header=TRUE)
dat1$lengthx <- dat1$length/10
dat1$avelengthx <- ave(dat1$lengthx, dat1$groupid)
dat1$difflengthx <- dat1$lengthx - dat1$avelengthx
dat1$avegold <- ave(dat1$goldcolor, dat1$groupid)
dat1$diffgold <- dat1$goldcolor - dat1$avegold

library(lme4)
library(r2mlm)

out1b0 <- lmer(consume ~ 1 + (1|groupid), data=dat1, REML=FALSE)
r2mlm(out1b0)

out1b1 <- lmer(consume ~ 1 + difflengthx + (1|groupid), data=dat1, REML=FALSE)
out1b2 <- lmer(consume ~ 1 + diffgold + (1|groupid), data=dat1, REML=FALSE)

anova(out1b0, out1b1)
r2mlm(out1b1)

anova(out1b0, out1b2)
r2mlm(out1b2)

out1b3 <- lmer(consume ~ 1 + difflengthx + diffgold + (1|groupid), data=dat1, REML=FALSE)
out1b4 <- lmer(consume ~ 1 + difflengthx + diffgold + difflengthx:diffgold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b3, out1b4)
r2mlm(out1b3)
r2mlm(out1b4)

# Use out1b1

out1b5 <- lmer(consume ~ 1 + difflengthx + volume + (1|groupid), data=dat1, REML=FALSE)
out1b6 <- lmer(consume ~ 1 + difflengthx + nfish + (1|groupid), data=dat1, REML=FALSE)
out1b7 <- lmer(consume ~ 1 + difflengthx + avelengthx + (1|groupid), data=dat1, REML=FALSE)
out1b8 <- lmer(consume ~ 1 + difflengthx + avegold + (1|groupid), data=dat1, REML=FALSE)

anova(out1b1, out1b5)
anova(out1b1, out1b6)
anova(out1b1, out1b7)
anova(out1b1, out1b8)

out1b9 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + (1|groupid), data=dat1, REML=FALSE)
out1b10 <- lmer(consume ~ 1 + difflengthx + avelengthx*volume + (1|groupid), data=dat1, REML=FALSE)
anova(out1b9, out1b10)

out1b11 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1|groupid), data=dat1, REML=FALSE)
out1b12 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx*nfish + (1|groupid), data=dat1, REML=FALSE)
out1b13 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume*nfish + (1|groupid), data=dat1, REML=FALSE)
anova(out1b11, out1b12)
anova(out1b11, out1b13)

out1b14 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx + avegold + (1|groupid), data=dat1, REML=FALSE)
out1b15 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx*avegold + (1|groupid), data=dat1, REML=FALSE)
out1b16 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume*avegold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b14, out1b15)
anova(out1b14, out1b16)

out1b17 <- lmer(consume ~ 1 + difflengthx + diffgold + volume + avelengthx + avegold + nfish + (1|groupid), data=dat1, REML=FALSE)
out1b18 <- lmer(consume ~ 1 + difflengthx + diffgold + volume + avelengthx + avegold*nfish + (1|groupid), data=dat1, REML=FALSE)
anova(out1b17, out1b18)

# Random slope then

out1b9 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + (1|groupid), data=dat1, REML=FALSE)
out1b19 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b9, out1b19)

out1b20 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + (1|groupid), data=dat1, REML=FALSE, 
                control = lmerControl(optimizer ="Nelder_Mead"))
out1b21 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + (1 + goldcolor|groupid), data=dat1, REML=FALSE, 
                control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1b20, out1b21)

# Use out1b19

out1b22 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:avelengthx + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
out1b23 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b19, out1b22)
anova(out1b19, out1b23)

out1b24 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
out1b25 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + difflengthx:nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b24, out1b25)

out1b26 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + avegold + (1 + difflengthx|groupid), data=dat1, REML=FALSE, 
                control = lmerControl(optimizer ="Nelder_Mead"))
out1b27 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + avegold + difflengthx:avegold + (1 + difflengthx|groupid), data=dat1, REML=FALSE, 
                control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1b26, out1b27)

# Use out1b23

# Tear-down strategy

out1t0 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish
               + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))

dat1a <- dat1
dat1a$difflengthx <- dat1a$difflengthx * 10
dat1a$avelength <- dat1a$avelength * 10
out1t1a <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avegold + difflengthx:volume + difflengthx:nfish
               + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1a, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out1t0a <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish
               + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1a, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))

out1t2 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:volume + difflengthx:nfish
               + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out1t3 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:avegold + difflengthx:nfish
               + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out1t4 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume
               + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out1t5 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish
               + diffgold:avegold + diffgold:volume + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out1t6 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish
               + diffgold:avelengthx + diffgold:volume + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out1t7 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish
               + diffgold:avelengthx + diffgold:avegold + diffgold:nfish
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out1t8 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish
               + diffgold:avelengthx + diffgold:avegold + diffgold:volume
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))

anova(out1t1a, out1t0a)
anova(out1t2, out1t0)
anova(out1t3, out1t0)
anova(out1t4, out1t0)
anova(out1t5, out1t0)
anova(out1t6, out1t0)
anova(out1t7, out1t0)
anova(out1t8, out1t0)

out1t9 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE)
out1t10 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + (1 + diffgold|groupid), data=dat1, REML=FALSE)
out1t11 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish 
               + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t10, out1t9)
anova(out1t11, out1t9)

# Use out1t11

out1t12 <- lmer(consume ~ 1 + difflengthx + avelengthx + avegold + volume + nfish 
                + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
out1t13 <- lmer(consume ~ 1 + difflengthx + diffgold + avegold + volume + nfish 
                + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
out1t14 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + nfish 
                + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
out1t15 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + nfish 
                + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
out1t16 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume 
                + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t12, out1t11)
anova(out1t13, out1t11)
anova(out1t14, out1t11)
anova(out1t15, out1t11)
anova(out1t16, out1t11)



dat2 <-read.table("lecture7ex1.csv", sep=",", header=TRUE)
# dat2$eeiq <- dat2$eeiq/100
# dat2$eeworkexp <- dat2$eeworkexp/10
# dat2$erexp <- dat2$erexp/10

library(lme4)
out2m0 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) 
               + ersex + erexp + eesex:ersex 
               + (1 + eesex|erid), data=dat2, REML=FALSE)
summary(out2m0)

out2m1 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) 
               + ersex + erexp + eesex:ersex 
               + (1 + eesex|erid), data=dat2, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out2m2 <- lmer(score ~ 1 + eesex + eesalary + I(eeiq - 100) 
               + ersex + erexp + eesex:ersex 
               + (1 + eesex|erid), data=dat2, REML=FALSE)
out2m3 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp 
               + ersex + erexp + eesex:ersex 
               + (1 + eesex|erid), data=dat2, REML=FALSE)
out2m4 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) 
               + ersex + eesex:ersex 
               + (1 + eesex|erid), data=dat2, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m1, out2m0)
anova(out2m2, out2m0)
anova(out2m3, out2m0)
anova(out2m4, out2m0)


out2m5 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) 
               + ersex + eesex:ersex 
               + (1 + eesex|erid), data=dat2, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out2m6 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) 
               + ersex + eesex:ersex 
               + (1 + eesex + eeworkexp|erid), data=dat2, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out2m7 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq/10 - 10) 
               + ersex + eesex:ersex 
               + (1 + eesex + I(eeiq/10 - 10)|erid), data=dat2, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))

anova(out2m5, out2m6)
anova(out2m5, out2m7)

out2m8 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) + eesalary
               + ersex + eesex:ersex 
               + (1 + eesex|erid), data=dat2, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
out2m9 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) + eesalary
               + ersex + eesex:ersex 
               + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m8, out2m9)

# out2b9

out2m10 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) + eesalary
               + ersex + eesex:ersex + eesalary:ersex
               + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m9, out2m10)

out2m11 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) + eesalary
                + ersex + erexp + eesex:ersex
                + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, 
                control = lmerControl(optimizer ="Nelder_Mead"))
out2m12 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) + eesalary
                + ersex + erexp + eesex:ersex + eesex:erexp
                + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, 
                control = lmerControl(optimizer ="Nelder_Mead"))
out2m13 <- lmer(score ~ 1 + eesex + eeworkexp + I(eeiq - 100) + eesalary
                + ersex + erexp + eesex:ersex + eesalary:erexp
                + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, 
                control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m11, out2m12)
anova(out2m11, out2m13)

# Keey out2b9
summary(out2m9)
 

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

dat1 <-read.table("lecture4ex1.csv", sep=",", header=TRUE)
dat1$lengthx <- dat1$length/10
dat1$avelengthx <- ave(dat1$lengthx, dat1$groupid)
dat1$difflengthx <- dat1$lengthx - dat1$avelengthx
dat1$avegold <- ave(dat1$goldcolor, dat1$groupid)
dat1$diffgold <- dat1$goldcolor - dat1$avegold

library(lme4)
out1b0 <- lmer(consume ~ 1 + (1|groupid), data=dat1, REML=FALSE)
summary(out1b0)
out1b9 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume 
               + (1|groupid), data=dat1, REML=FALSE)
summary(out1b9)

out1b19 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume 
             + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
summary(out1b19)
out1b91 <- lmer(consume ~ 1 + difflengthx + avelengthx 
                + (1|groupid), data=dat1, REML=FALSE)
summary(out1b91)
out1b92 <- lmer(consume ~ 1 + volume 
                + (1|groupid), data=dat1, REML=FALSE)
summary(out1b92)
out1b23 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume 
                + difflengthx:volume + (1 + difflengthx|groupid), 
                data=dat1, REML=FALSE)
summary(out1b23)


out1i0 <- lmer(consume ~ 1 + (1|groupid), data=dat1, REML=FALSE)
summary(out1i0)
out1i1 <- lmer(consume ~ 1  + difflengthx + (1|groupid), data=dat1, REML=FALSE)
summary(out1i1)
out1i2 <- lmer(consume ~ 1  + avelengthx + (1|groupid), data=dat1, REML=FALSE)
summary(out1i2)
out1i12 <- lmer(consume ~ 1  + difflengthx + avelengthx + (1|groupid), data=dat1, REML=FALSE)
summary(out1i12)
out1n1 <- lmer(consume ~ 1  + lengthx + (1|groupid), data=dat1, REML=FALSE)
summary(out1n1)
out1n12 <- lmer(consume ~ 1  + lengthx + avelengthx + (1|groupid), data=dat1, REML=FALSE)
summary(out1n12)

alln <- table(dat1$groupid)
bign <- sum(alln)
k <- length(alln)
ntilde <- (bign^2 - sum(alln^2))/(bign * (k-1))
ntilde


dat1$zconsume <- scale(dat1$consume)
dat1$zlength <- scale(dat1$length)
dat1$avezlength <- ave(dat1$zlength, dat1$groupid)
dat1$diffzlength <- dat1$zlength - dat1$avezlength

volume <- dat1[!duplicated(dat1$groupid), "volume"]
mvolume <- mean(volume)
sdvolume <- sd(volume)
dat1$zvolume <- (dat1$volume - mvolume) / sd(volume)

out1z23 <- lmer(zconsume ~ 1 + diffzlength + avezlength + zvolume 
                + diffzlength:zvolume + (1 + diffzlength|groupid), 
                data=dat1, REML=FALSE)
summary(out1z23)

out1zb23 <- lmer(consume ~ 1 + diffzlength + avezlength + zvolume 
                + diffzlength:zvolume + (1 + diffzlength|groupid), 
                data=dat1, REML=FALSE)
summary(out1zb23)
out1bz23 <- lmer(zconsume ~ 1 + difflengthx + avelengthx + volume 
                + difflengthx:volume + (1 + difflengthx|groupid), 
                data=dat1, REML=FALSE)
summary(out1bz23)





############ Using MSEM

library(lavaan)
dat <-read.table("lecture7ex3.csv", sep=",", header=TRUE)

model1 <- '
level: 1
sat ~ groupflow
level: 2
sat ~ groupflow
'

fit1 <- sem(model = model1, data = dat, cluster = "groupid")
summary(fit1)

dat$aveflow <- ave(dat$groupflow, dat$groupid)
dat$diffflow <- dat$groupflow - dat$aveflow
library(lme4)
out2 <- lmer(sat ~ 1 + diffflow + aveflow + (1|groupid), data=dat, REML=FALSE)
summary(out2)

# Formative Model
model2 <- '
level: 1
sat ~ diffflow
level: 2
sat ~ aveflow
'
fit2 <- sem(model = model2, data = dat, cluster = "groupid")
summary(fit2)

model3 <- '
level: 1
sat ~ emostability
level: 2
sat ~ emostability
'
fit3 <- sem(model = model3, data = dat, cluster = "groupid")
summary(fit3)

dat$aveemo <- ave(dat$emostability, dat$groupid)
dat$diffemo <- dat$emostability - dat$aveemo
out4 <- lmer(sat ~ 1 + diffemo + aveemo + (1|groupid), data=dat, REML=FALSE)
summary(out4)

model4 <- '
level: 1
sat ~ diffemo
level: 2
sat ~ aveemo
'
fit4 <- sem(model = model4, data = dat, cluster = "groupid")
summary(fit4)


# Compare fit1 with mlm

model5 <- '
level: 1
groupflow ~ aw*perceivedability
sat ~ bw*groupflow + perceivedability
level: 2
groupflow ~ ab*perceivedability + experience
sat ~ bb*groupflow + perceivedability + experience
abw := aw*bw
abb := ab*bb
'
fit5 <- sem(model = model5, data = dat, cluster = "groupid")
summary(fit5)

model6 <- '
level: 1
sat ~~ emostability + groupflow + perceivedability
emostability ~~ groupflow + perceivedability
groupflow ~~ perceivedability
level: 2
sat ~~ emostability + groupflow + perceivedability + experience
emostability ~~ groupflow + perceivedability + experience
groupflow ~~ perceivedability + experience
perceivedability ~~ experience
'
fit6 <- sem(model = model6, data = dat, cluster = "groupid")
summary(fit6, standardize=TRUE)

standardizedsolution(fit6)

### Find ICC from ANOVA

iccanova <- function(y, x, data) {
  data[,x] <- factor(data[,x])
  f <- paste(y, "~", x)
  out <- do.call("aov", list(as.formula(f), data=as.name("data")))
  result <- summary(out)[[1]]
  msb <- result[1, "Mean Sq"]
  msw <- result[2, "Mean Sq"]
  dfb <- result[1, "Df"]
  icc1 <- (msb - msw)/(msb + dfb*msw)
  icc2 <- (msb - msw)/msb
  c(icc1=icc1, icc2=icc2)
}

library(lavaan)

model1 <- '
    level: 1
        consume ~~ length
    level: 2
        consume ~~ length + nfish + volume
        length ~~ nfish + volume
        nfish ~~ volume
'
fit1 <- sem(model = model1, data = dat1, cluster = "groupid")
summary(fit1, standardize=TRUE)


icciq2 <- iccanova("iq", "erid", data=dat2)

# IQ have very low ICC so fix L2 variance of IQ as 0
model2 <- '
    level: 1
        score ~~ iq
    level: 2
        score ~~ score
        iq ~~ 0*iq
'
fit2 <- sem(model = model2, data = dat2, cluster = "erid")
summary(fit2, standardize=TRUE)





iccage3 <- iccanova("age", "tableid", data=dat3)

model3 <- '
    level: 1
        sat ~~ age
    level: 2
        sat ~~ numperson
        age ~~ 0*age
'
fit3 <- sem(model = model3, data = dat3, cluster = "tableid")
summary(fit3, standardize=TRUE)