setwd("C:/Users/Gary/Dropbox/Post doc/R workshop//Workshop 2/example")
setwd("C:/Users/t8s/Documents/Dropbox/Post doc/R workshop/Workshop 2/example")

#install.packages("nlme")

library(lme4)
#library(nlme)
library(lattice)

alcohol <- read.csv("alcohol_pp.csv", header = TRUE, sep =",")

xyplot(alcuse ~ age | id, data = alcohol[(1:27),], 
       panel = function(x,y){
         panel.xyplot(x,y)
         panel.lmline(x,y)
       })


#gm.reg1 <- lmer(alcuse ~ 1 + (1 | id), data = alcohol, REML = FALSE)
#summary(gm.reg1)

gm.reg1 <- lmer(alcuse ~ age_14 + (1| id), data = alcohol, REML = FALSE)
summary(gm.reg1)

gm.reg2 <- lmer(alcuse ~ age_14 + (1 + age_14| id), data = alcohol, REML = FALSE)
summary(gm.reg2)

anova(gm.reg1, gm.reg2)
gm.reg2.cov <- VarCorr(gm.reg2)
gm.reg2.cov[1]

alcohol$predscore.gm.reg2 <- fitted(gm.reg2)
xyplot(predscore.gm.reg2 ~ age_14, data = alcohol, groups = id, type = c("p","l"))

xyplot(predscore.gm.reg2 ~ age_14, 
       data = alcohol[alcohol$id %in% c(1:15),], groups = id, type = c("p","l"))

#the seq function generates a sequence of numbers
#from 0 to 2, with an increment of 0.1
x <- seq(0,2,0.1)
y <- 0.6244+2*(-0.0684)*x+0.1512*x^2
plot(y ~ x, xlab = "age14", ylab = "Between-individual variance", type = "l")
  
re.gm.reg2 <- ranef(gm.reg2, condVar = TRUE)
plot(re.gm.reg2[[1]], xlab = "Intercept", ylab = "Slope of Age")
abline(h = 0, col = "red")
abline(v = 0, col = "red")

gm.reg3 <- lmer(alcuse ~ age_14 + coa + age_14:coa + peer + peer:age_14 + (1 + age_14 | id), data = alcohol, REML = FALSE)
summary(gm.reg3)

gm.reg4 <- lmer(alcuse ~ age_14 + coa + peer + peer:age_14 + (1 + age_14 | id), data = alcohol, REML = FALSE)
summary(gm.reg4)

#Predicting alcohol use at different level of peer use and COA
#age is centered to 14 in the model
predictor_age14 <- c(0:2)

#calculate the predicted values
low_peer_not_coa <- -0.112+0.381*predictor_age14
low_peer_coa <- 0.458+0.381*predictor_age14
high_peer_not_coa <- 0.897+0.162*predictor_age14
high_peer_coa <- 1.468 + 0.162*predictor_age14

#original coding for age - 14,15 and 16
predictor_age <- c(14:16)

#plot the predicted alcohol consumption for different level of peer use, COA and age.
plot(low_peer_not_coa ~ predictor_age, type = "b", 
     ylab = "Predicted alcohol use", xlab = "Age", pch = 2, ylim = c(-0.8,2))
lines(low_peer_coa ~ predictor_age, type = "b", pch = 0)
lines(high_peer_not_coa ~ predictor_age, type = "b", pch = 17)
lines(high_peer_coa ~ predictor_age, type = "b", pch = 15)
legend("bottomright", pch= c(2,0,17,15), c("Low peer, Not COA",
                                           "Low peer, COA", "High peer, Not COA","High peer, COA"))

#the level 1 residuals can be extracted using the residuals function
resid <- residuals(gm.reg4)
qqmath(resid)

resid.std <- resid/sd(resid)
plot(alcohol$id, resid.std)
abline(h = 0, col = "red")

#use the ranef to extract the individual lv 2 residuals (random effect)
ran <- ranef(gm.reg4)
ran <- ran[[1]]

#the random effects are defined for EACH individuals
#we need to extract only 1 record for each individuals
#and plot residuals against participant ID
unique.alcohol <- alcohol[!duplicated(alcohol$id),]

ran1.std <- ran$`(Intercept)`/sd(ran$`(Intercept)`)
plot(unique.alcohol$id, ran1.std)
abline(h = 0, col = "red")

qqmath(ran1.std)

#slope residual
ran2.std <- ran$age_14/sd(ran$age_14)
plot(unique.alcohol$id, ran2.std)
abline(h = 0, col = "red")

qqmath(ran2.std)

par(mfrow = c(1,1))

#level 1 residuals vs age
plot(alcohol$age_14, resid)
abline(h = 0)

par(mfrow = c(2,2))
#random intercepts vs COA
plot(unique.alcohol$coa, ran1.std)
abline(h = 0)

#random slopes vs COA
plot(unique.alcohol$coa, ran2.std)
abline(h = 0)

#random intercepts vs peer
plot(unique.alcohol$peer, ran1.std)
abline(h = 0)

#random slopes vs peer
plot(unique.alcohol$peer, ran2.std)
abline(h = 0)

###################################################################
library(lme4)

#Treating time more flexibly

reading <- read.csv("reading_pp.csv", header = TRUE, sep =",")

xyplot(piat ~ age| id, data = reading[reading$id %in% c(1:9),],
       panel = function(x,y) {
         panel.xyplot(x,y, pch = 16)
         panel.lmline(x,y, lty = 4)
       })

xyplot(piat ~ agegrp | id, data = reading[reading$id %in% c(1:9),],
       panel = function(x,y) {
         panel.xyplot(x,y, pch = 16)
         panel.lmline(x,y, lty = 4)
       })

reading$agegrp.c <- reading$agegrp - 6.5
reading$age.c <- reading$age - 6.5

gm.reading1 <- lmer(piat ~ age.c + (1 + age.c | id), data = reading, REML = FALSE)
summary(gm.reading1)

gm.reading2 <- lmer(piat ~ agegrp.c + (1 + agegrp.c | id), data = reading, REML = FALSE)
summary(gm.reading2)

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

unemploy <- read.csv("unemployment_pp.csv", head = TRUE, sep =",")

gm.unemp1 <- lmer(cesd ~ months + (1 + months | id), data = unemploy, REML = FALSE)
summary(gm.unemp1)

gm.unemp2 <- lmer(cesd ~ months + unemp + (1 + months | id), data = unemploy, REML = FALSE)
summary(gm.unemp2)

gm.unemp3 <- lmer(cesd ~ months + unemp + months:unemp + (1 + months | id), data = unemploy, REML = FALSE)
summary(gm.unemp3)

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

phf_data <- read.csv("phf.csv", head = TRUE, sep = ",")

xyplot(phf ~ age50 | id, data = phf_data[phf_data$id %in% c(1:10),], 
       panel = function(x,y){
         panel.xyplot(x,y)
         panel.lmline(x,y)
       })

gm.phf1 <- lmer(phf ~ age50 + (1 | id), data = phf_data, REML = FALSE)
summary(gm.phf1)

gm.phf2 <- lmer(phf ~ age50 + (1 + age50 | id), data = phf_data, REML = FALSE)
summary(gm.phf2)

k <- VarCorr(gm.phf2)
k[[1]]

x <- seq(-20,30,1)
y <- 28.819+2*0.42*x + 0.0922*x^2
plot(y ~ x, type = "l")

gm.phf3 <- lmer(phf ~ age50 + female + grade + grade:age50 + female:age50 + (1 + age50 | id), data = phf_data, REML = FALSE)
summary(gm.phf3)

gm.phf4 <- lmer(phf ~ age50 + female + grade + grade:age50 + (1 + age50 | id), data = phf_data, REML = FALSE)
summary(gm.phf4)

predictor_age50 <- seq(-11,26,5)

M_high <- 53.27-0.212*predictor_age50
M_intermediate <- 52.593-0.248*predictor_age50
M_low <- 50.666-0.325*predictor_age50
F_high <- 51.533-0.212*predictor_age50
F_intermediate <- 50.856-0.248*predictor_age50
F_low <- 48.929 - 0.325*predictor_age50

predictor_age <- seq(39,76,5)

#plot the predicted alcohol consumption for different level of peer use, COA and age.
plot(M_high ~ predictor_age, type = "b", 
     ylab = "Predicted physical functioning", xlab = "Age", pch = 0, ylim = c(40,60))
lines(M_intermediate ~ predictor_age, type = "b", pch = 1)
lines(M_low ~ predictor_age, type = "b", pch = 2)
lines(F_high ~ predictor_age, type = "b", pch = 15)
lines(F_intermediate ~ predictor_age, type = "b", pch = 16)
lines(F_low ~ predictor_age, type = "b", pch = 17)
legend("topright", pch= c(0,1,2,15,16,17), c("Male, High grade",
                                           "Male, Intermediate grade", "Male, Low grade",
                                           "Female, High grade", "Female, Intermediate grade",
                                           "Female, Low grade"))

phf_data$predscore.gm.phf4 <- fitted(gm.phf4)

xyplot(predscore.gm.phf4 ~ age, 
       data = phf_data[phf_data$id %in% c(1:50),], groups = id, type = c("p","l"))

#the level 1 residuals can be extracted using the residuals function
resid <- residuals(gm.phf4)
qqmath(resid)

resid.std <- resid/sd(resid)
plot(phf_data$id, resid.std)
abline(h = 0, col = "red")

#use the ranef to extract the individual lv 2 residuals (random effect)
ran <- ranef(gm.phf4)
ran <- ran[[1]]

unique.phf <- phf_data[!duplicated(phf_data$id),]

ran1.std <- ran$`(Intercept)`/sd(ran$`(Intercept)`)
plot(unique.phf$id, ran1.std)
abline(h = 0, col = "red")

qqmath(ran1.std)

#slope residual
ran2.std <- ran$age50/sd(ran$age50)
plot(unique.phf$id, ran2.std)
abline(h = 0, col = "red")

qqmath(ran2.std)

par(mfrow = c(1,1))

#level 1 residuals vs age
plot(phf_data$age50, resid)
abline(h = 0)

par(mfrow = c(2,2))
#random intercepts vs female
plot(unique.phf$female, ran1.std)
abline(h = 0)

#random slopes vs female
plot(unique.phf$female, ran2.std)
abline(h = 0)

#random intercepts vs grade
plot(as.numeric(unique.phf$grade), ran1.std)
abline(h = 0)

#random slopes vs peer
plot(as.numeric(unique.phf$grade), ran2.std)

abline(h = 0)

##########################################################################################
phf_data$phf2 <- phf_data$phf^2

gm.phf5 <- lmer(phf2 ~ age50 + female + grade:age50 + (1 + age50 | id), data = phf_data, REML = FALSE)
summary(gm.phf5)

#the level 1 residuals can be extracted using the residuals function
resid <- residuals(gm.phf5)
qqmath(resid)

resid.std <- resid/sd(resid)
plot(phf_data$id, resid.std)
abline(h = 0, col = "red")

#use the ranef to extract the individual lv 2 residuals (random effect)
ran <- ranef(gm.phf5)
ran <- ran[[1]]

unique.phf <- phf_data[!duplicated(phf_data$id),]

ran1.std <- ran$`(Intercept)`/sd(ran$`(Intercept)`)
plot(unique.phf$id, ran1.std)
abline(h = 0, col = "red")

qqmath(ran1.std)

#slope residual
ran2.std <- ran$age50/sd(ran$age50)
plot(unique.phf$id, ran2.std)
abline(h = 0, col = "red")

qqmath(ran2.std)


#########################################################################################
# Non-linear trend
#########################################################################################

phf_data$age50_2 <- phf_data$age50^2
gm.phf6 <- lmer(phf ~ age50 + age50_2 + (1 + age50 + age50_2 | id), data = phf_data, REML = FALSE)
summary(gm.phf6)

gm.phf7 <- lmer(phf ~ age50 + age50_2 + (1 + age50 | id), data = phf_data, REML = FALSE)
summary(gm.phf7)

phf_data$predscore.gm.phf7 <- fitted(gm.phf7)
xyplot(predscore.gm.phf7 ~ age50, data = phf_data, groups = id, type = c("p","l"))

xyplot(predscore.gm.phf7 ~ age50, data = phf_data[phf_data$id %in% c(1:40),], groups = id, type = c("p","l"))
