dat1 <-read.table("lecture8ex1.csv", sep=",", header=TRUE)
psych::describe(dat1)
library(lme4)
out1a <- lmer(jsat ~ 1 + (1|pid), data=dat1, REML=FALSE)
summary(out1a)
out1b <- lmer(jsat ~ 1 + time + (1|pid), data=dat1, REML=FALSE)
summary(out1b)
anova(out1a, out1b)

out1c <- lmer(jsat ~ 1 + time + (1 + time|pid), data=dat1, REML=FALSE)
summary(out1c)

dat1$timel <- dat1$time - 3
out1c1 <- lmer(jsat ~ 1 + timel + (1 + timel|pid), data=dat1, REML=FALSE)
summary(out1c1)

anova(out1b, out1c)


library(ggplot2)
dat1_2 <- dat1[dat1$pid%%10 == 0,]
ggplot(dat1_2, aes(x=time, y=jsat, group=pid)) + geom_smooth(method=lm, se=FALSE)

dat1_2 <- dat1[dat1$pid%%3 == 0,]
dat1_2$age <- dat1_2$age0 + dat1_2$time
ggplot(dat1_2, aes(x=age, y=jsat, group=pid)) + geom_smooth(method=lm, se=FALSE)

dat1_2 <- dat1[dat1$pid%%5 == 0,]
dat1_2$tenure <- dat1_2$tenure0 + dat1_2$time
ggplot(dat1_2, aes(x=tenure, y=jsat, group=pid)) + geom_smooth(method=lm, se=FALSE)


# cor(coef(out2)$pid)

out1d <- lmer(jsat ~ 1 + time + married + (1 + time|pid), data=dat1, REML=FALSE)
summary(out1d)

dat1$avemarried <- ave(dat1$married, dat1$pid)
dat1$diffmarried <- dat1$married - dat1$avemarried
out1e <- lmer(jsat ~ 1 + time + diffmarried + avemarried + (1 + time|pid), data=dat1, REML=FALSE)
summary(out1e)

out1g <- lmer(jsat ~ 1 + time + married + avemarried + (1 + time|pid), data=dat1, REML=FALSE)
summary(out1e)

out1h <- lmer(jsat ~ 1 + time + avemarried + I(age0-25) + I(tenure0-5) 
              + fftj + officejob + (1 + time|pid), data=dat1, REML=FALSE)
summary(out1h)


out1i <- lmer(jsat ~ 1 + time + time*avemarried + I(age0-25) + I(tenure0-5) 
              + fftj + officejob + (1 + time|pid), data=dat1, REML=FALSE)
out1j <- lmer(jsat ~ 1 + time + avemarried + time*I(age0-25) + I(tenure0-5) 
              + fftj + officejob + (1 + time|pid), data=dat1, REML=FALSE)
out1k <- lmer(jsat ~ 1 + time + avemarried + I(age0-25) + time*I(tenure0-5) 
              + fftj + officejob + (1 + time|pid), data=dat1, REML=FALSE)
out1l <- lmer(jsat ~ 1 + time + avemarried + I(age0-25) + I(tenure0-5) 
              + time*fftj + officejob + (1 + time|pid), data=dat1, REML=FALSE)
out1m <- lmer(jsat ~ 1 + time + avemarried + I(age0-25) + I(tenure0-5) 
              + fftj + time*officejob + (1 + time|pid), data=dat1, REML=FALSE)
anova(out1h, out1i)
anova(out1h, out1j)
anova(out1h, out1k)
anova(out1h, out1l)
anova(out1h, out1m)

dat1a <- dat1[!duplicated(dat1$pid),]
apply(dat2a, 2, mean)

out1n <- lmer(jsat ~ 1 + time + time*I(avemarried-0.3725) 
              + time*I(age0-28.375) + time*I(tenure0-2.054) 
              + time*I(fftj-0.4625) + time*I(officejob-0.33) 
              + (1 + time|pid), data=dat1, REML=FALSE)
summary(out1n)

# Probing Interaction between time and age0

plot(dat1$time, dat1$jsat, type="n", xlab="Year of Measurement", ylab="Job Satisfaction")
abline(a = 47.31, b = -1.22, col="red")
abline(a = 50.21, b = -0.52, col="green")
abline(a = 53.11, b = 0.18, col="blue")
legend(0, 90, c(25, 35, 45), col=c("red", "green", "blue"), 
       lty=1, title="Age at Year 0")

plot(dat1$age0, dat1$jsat, type="n", xlab="Age at Year 0", ylab="Job Satisfaction")
abline(a = 49.27 - (28.375*0.58), b = 0.58, col="red")
abline(a = 48.52 - (28.375*0.72), b = 0.72, col="green")
abline(a = 47.77 - (28.375*0.86), b = 0.86, col="blue")
abline(a = 47.02 - (28.375*1.00), b = 1.00, col="black")
legend(18, 90, c("0", "1", "2", "3"), col=c("red", "green", "blue", "black"), 
       lty=1, title="Year of Measurement")


sumout1n <- summary(out1n)
coef(sumout1n)

psi <- coef(sumout1n)[,1]
vpsi <- vcov(sumout1n)
cvec1 <- c(1, 0, 0, -3.375, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 0, 0, 1.625, 0, 0, 0, 0, 0, 0, 0, 0)  
cvec3 <- c(1, 0, 0, 6.625, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec4 <- c(0, 1, 0, 0, 0, 0, 0, 0, -3.375, 0, 0, 0) 
cvec5 <- c(0, 1, 0, 0, 0, 0, 0, 0, 1.625, 0, 0, 0) 
cvec6 <- c(0, 1, 0, 0, 0, 0, 0, 0, 6.625, 0, 0, 0)  
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, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)  
cvec3 <- c(1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec4 <- c(1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec5 <- c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec6 <- c(0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0)  
cvec7 <- c(0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0)  
cvec8 <- c(0, 0, 0, 1, 0, 0, 0, 0, 3, 0, 0, 0)  
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6, cvec7, cvec8)
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)

# Probing interaction between time and tenure0

plot(dat1$time, dat1$jsat, type="n", xlab="Year of Measurement", ylab="Job Satisfaction")
abline(a = 52.06, b = -2.89, col="red")
abline(a = 49.34, b = -0.81, col="green")
abline(a = 46.62, b = 1.27, col="blue")
legend(0, 90, c(0, 2, 4), col=c("red", "green", "blue"), 
       lty=1, title="Tenure at Year 0")

plot(dat1$tenure0, dat1$jsat, type="n", xlab="Tenure at Year 0", ylab="Job Satisfaction")
abline(a = 49.27 - (2.054*-1.36), b = -1.36, col="red")
abline(a = 48.52 - (2.054*-0.32), b = -0.32, col="green")
abline(a = 47.77 - (2.054*0.72), b = 0.72, col="blue")
abline(a = 47.02 - (2.054*1.76), b = 1.76, col="black")
legend(0, 90, c("0", "1", "2", "3"), col=c("red", "green", "blue", "black"), 
       lty=1, title="Year of Measurement")

cvec1 <- c(1, 0, 0, 0, -2.054, 0, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 0, 0, 0, -0.054, 0, 0, 0, 0, 0, 0, 0)  
cvec3 <- c(1, 0, 0, 0, 1.946, 0, 0, 0, 0, 0, 0, 0) 
cvec4 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, -2.054, 0, 0) 
cvec5 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, -0.054, 0, 0) 
cvec6 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 1.946, 0, 0)  
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, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)  
cvec3 <- c(1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec4 <- c(1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec5 <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0) 
cvec6 <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0)  
cvec7 <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0)  
cvec8 <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 3, 0, 0)  
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6, cvec7, cvec8)
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)

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

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

library(lme4)
out2a <- lmer(mem ~ 1 + (1|pid), data=dat2, REML=FALSE)
summary(out2a)

dat2$timec <- dat2$time - 1

dat2$calyear <- dat2$time - 3
dat2[dat2$cohort == 0, "calyear"] <- dat2[dat2$cohort == 0, "calyear"] + 1

out2b <- lmer(mem ~ 1 + timec + (1|pid), data=dat2, REML=FALSE)
summary(out2b)

out2b1 <- lmer(mem ~ 1 + calyear + (1|pid), data=dat2, REML=FALSE)
summary(out2b1)

out2c <- lmer(mem ~ 1 + time + I(time^2) + (1|pid), data=dat2, REML=FALSE)
summary(out2c)

mytime <- 0:3
mypred <- 64.7526 - 2.7848*mytime + 0.1632*mytime^2
plot(mytime, mypred, type="l")

out2d <- lmer(mem ~ 1 + timec + (1 + timec|pid), data=dat2, REML=FALSE)
anova(out2b, out2d)
summary(out2d)

dat2_2 <- dat2[dat2$pid%%5 == 0,]
library(ggplot2)
dat2_2$cohort <- factor(dat2_2$cohort, labels=c("2551", "2550"))
ggplot(dat2_2, aes(x=time, y=mem, group=pid, colour=cohort)) + geom_smooth(method=lm, se=FALSE, size=0.5)

dat2_2 <- dat2[dat2$pid%%5 == 0,]
library(ggplot2)
dat2_2$cohort <- factor(dat2_2$cohort, labels=c("2551", "2550"))
dat2_2$calyear <- dat2_2$calyear + 2552
ggplot(dat2_2, aes(x=calyear, y=mem, group=pid, colour=cohort)) + geom_smooth(method=lm, se=FALSE, size=0.5)



dat2$avestress <- ave(dat2$stress, dat2$pid)
dat2$diffstress <- dat2$stress - dat2$avestress
out2e <- lmer(mem ~ 1 + timec + diffstress + avestress + (1 + timec|pid), data=dat2, REML=FALSE)
summary(out2e)

# out2f <- lmer(mem ~ 1 + timec + diffstress + avestress + (1 + timec + diffstress|pid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))

dat2a <- dat2[!duplicated(dat2$pid),]
apply(dat2a, 2, mean)

out2g <- lmer(mem ~ 1 + timec + diffstress + I(avestress - 50.473) 
              + I(cohort - 0.501) + I(female - 0.495) 
              + I(act - 50.219) + I(ext - 49.534) 
              + (1 + timec|pid), data=dat2, REML=FALSE)
summary(out2g)

out2h <- lmer(mem ~ 1 + time + diffstress + time*I(avestress - 50.473) 
              + time*I(cohort - 0.501) + time*I(female - 0.495) 
              + time*I(act - 50.219) + time*I(ext - 49.534) 
              + (1 + time|pid), data=dat2, REML=FALSE)
summary(out2h)

# Probing Interaction between time and cohort

plot(dat2$time, dat2$mem, type="n", xlab="Year of University", ylab="Feeling of Membership")
abline(a = 66.54 - (-2.95), b = -2.95, col="red")
abline(a = 60.93 - (-0.80), b = -0.80, col="blue")
legend(1, 110, c(2550, 2551), col=c("red", "blue"), 
       lty=1, title="Class of")

sumout2h <- summary(out2h)
coef(sumout2h)

psi <- coef(sumout2h)[,1]
vpsi <- vcov(sumout2h)
cvec1 <- c(1, 0, 0, 0, 0.499, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 0, 0, 0, -0.501, 0, 0, 0, 0, 0, 0, 0, 0)   
cvec3 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0.499, 0, 0, 0) 
cvec4 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, -0.501, 0, 0, 0) 
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, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)  
cvec3 <- c(1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec4 <- c(1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec5 <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec6 <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0)  
cvec7 <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0)  
cvec8 <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 3, 0, 0, 0)  
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6, cvec7, cvec8)
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)

# Probing Interaction between time and activity

plot(dat2$time, dat2$mem, type="n", xlab="Year of University", ylab="Feeling of Membership")
abline(a = 64.35 - (-5.56), b = -5.56, col="red")
abline(a = 63.75 - (-1.96), b = -1.96, col="green")
abline(a = 63.15 - (1.64), b = 1.64, col="blue")
legend(1, 110, c(40, 50, 60), col=c("red", "green", "blue"), 
       lty=1, title="High School Activity")

plot(dat2$act, dat2$mem, type="n", xlab="High School Activity", ylab="Feeling of Membership")
abline(a = 63.74 - (50.219*-0.06), b = -0.06, col="red")
abline(a = 61.86 - (50.219*0.30), b = 0.30, col="green")
abline(a = 59.98 - (50.219*0.66), b = 0.66, col="blue")
abline(a = 58.10 - (50.219*1.02), b = 1.02, col="black")
legend(60, 50, c("1", "2", "3", "4"), col=c("red", "green", "blue", "black"), 
       lty=1, title="Year of University")


cvec1 <- c(1, 0, 0, 0, 0, 0, -10.219, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 0, 0, 0, 0, 0, -0.219, 0, 0, 0, 0, 0, 0)   
cvec3 <- c(1, 0, 0, 0, 0, 0, 9.781, 0, 0, 0, 0, 0, 0)   
cvec4 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -10.219, 0) 
cvec5 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.219, 0) 
cvec6 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9.781, 0) 
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, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)  
cvec3 <- c(1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec4 <- c(1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec5 <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0) 
cvec6 <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0)  
cvec7 <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0)  
cvec8 <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 3, 0)  
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6, cvec7, cvec8)
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)

# Probing Interaction between time and extraversion

plot(dat2$time, dat2$mem, type="n", xlab="Year of University", ylab="Feeling of Membership")
abline(a = 59.93 - (0.12), b = 0.12, col="red")
abline(a = 63.93 - (-1.98), b = -1.98, col="green")
abline(a = 67.93 - (-4.08), b = -4.08, col="blue")
legend(3, 50, c(40, 50, 60), col=c("red", "green", "blue"), 
       lty=1, title="Extraversion")

plot(dat2$act, dat2$mem, type="n", xlab="Extraversion", ylab="Feeling of Membership")
abline(a = 63.74 - (49.534*0.40), b = 0.40, col="red")
abline(a = 61.86 - (49.534*0.19), b = 0.19, col="green")
abline(a = 59.98 - (49.534*-0.02), b = -0.02, col="blue")
abline(a = 58.10 - (49.534*-0.23), b = -0.23, col="black")
legend(50, 50, c("1", "2", "3", "4"), col=c("red", "green", "blue", "black"), 
       lty=1, title="Year of University")

cvec1 <- c(1, 0, 0, 0, 0, 0, 0, -9.534, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 0, 0, 0, 0, 0, 0, 0.466, 0, 0, 0, 0, 0)   
cvec3 <- c(1, 0, 0, 0, 0, 0, 0, 10.466, 0, 0, 0, 0, 0)   
cvec4 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -9.534) 
cvec5 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.466) 
cvec6 <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10.466) 
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, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec2 <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)  
cvec3 <- c(1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec4 <- c(1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
cvec5 <- c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0) 
cvec6 <- c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1)  
cvec7 <- c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2)  
cvec8 <- c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 3)  
cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6, cvec7, cvec8)
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)