#First save tlc.txt to a folder on your PC
#eg, mine is is saved in C:/Users/dzhao1/OneDrive - University of Oklahoma/LMM
#Set up working directory
setwd("C:/Users/dzhao1/OneDrive - University of Oklahoma/LMM")
TLC.wide=read.table("tlc.txt", col.names=c("id","trt","y0","y1","y4","y6"))
#reshape the data from wide format to long format
TLC.long <- reshape(TLC.wide, idvar="id", varying=c("y0","y1","y4","y6"), 
                        v.names="y", timevar="time",time=c(0,1,4,6), direction="long")
#week is the categorical version of time
TLC.long$week=as.factor(TLC.long$time)
#make 'A' the reference level of trt
TLC.long$trt <- factor(TLC.long$trt, levels = c("P","A"))
attach(TLC.long)

# Spaghetti Plot of individiual trajectories
mean.na.rm=function(x) mean(x,na.rm=T)
interaction.plot(time, id, y, type="b", fun=mean.na.rm,
                 xlab="Time (in weeks)", ylab="Blood Lead Levels", 
                 main="Spaghetti Plot", col=c(1:100), legend=F)
#mean profile plot
interaction.plot(time, trt, y, type="b", fun=mean.na.rm,
                 xlab="Time (in weeks)", ylab="Blood Lead Levels", 
                 main="Mean Response Plot")
#######################################################
#install package
install.packages("lme4") 
#load package
library("lme4")
#Fit a LMM with fixed week, trt, and random patient,REML estimation
TLC1 = lmer(y ~ week + trt + (1 | id), TLC.long, REML=T)
#Summary of the results (fixed effects, AIC, loglik, etc)
summary(TLC1)
#Get 95% CI for fixed effects and variances
confint(TLC1)
#estimated random effects
ranef(TLC1)
#individual effect
coef(TLC1)
#individual prediction (fixed effects + random effects)
predict(TLC1)
#Mean prediction (fixed effects, without random effects)
predict(TLC1, re.form=NA)
#Residual plot
plot(TLC1)
qqnorm(resid(TLC1))
qqline(resid(TLC1))