setwd("D:/Courses/Longitudinal Data Analysis/Data Sets/")

eye.data <- read.table("visual_acuity_CH30.txt",header=T,na.strings = "-99")
attach(eye.data)

plot(Trt[Subject==1],Reac.Time[Subject==1],
 type="l",ylim=c(min(Reac.Time),max(Reac.Time)))
for (i in 2:max(Subject)){
  lines(Trt[Subject==i],Reac.Time[Subject==i])
}

fit <- aov(Reac.Time ~ factor(Trt)+Error(factor(Subject)))
summary(fit)

fit <- aov(Reac.Time ~ Sex+factor(Trt)+Error(factor(Subject)))
summary(fit)

fit <- aov(Reac.Time ~ Sex*Eye*factor(Lens.Strength)+Error(factor(Subject)))
summary(fit)


# For a more successful coding, consider:

Sex1 <- 1*(Sex=="Male")-1*(Sex=="Female")
LS1 <- -3*(Lens.Strength==6)-1*(Lens.Strength==18)+1*(Lens.Strength==36)+3*(Lens.Strength==60)
LS2 <- 1*(Lens.Strength==6)-1*(Lens.Strength==18)-1*(Lens.Strength==36)+1*(Lens.Strength==60)
LS3 <- -1*(Lens.Strength==6)+3*(Lens.Strength==18)-3*(Lens.Strength==36)+1*(Lens.Strength==60)
Eye1 <- 1*(Eye=="Left")-1*(Eye=="Right")

fit <- aov(Reac.Time ~ Sex1*Eye1*(LS1+LS2+LS3)+Error(factor(Subject)))
summary(fit)


# What happens if you ignore the repeated measures structure?

fit1 <- aov(Reac.Time ~ Sex1*Eye1*(LS1+LS2+LS3))
summary(fit1)


# We can also use the "linear mixed model" command.

library('nlme')

fit <- lme(Reac.Time ~ factor(Trt), random = ~1|Subject)
summary(fit)

fit <- lme(Reac.Time ~ Sex*Eye*factor(Lens.Strength), random = ~1|Subject)
summary(fit)

# For a more successful coding, consider:

Sex1 <- 1*(Sex=="Male")-1*(Sex=="Female")
LS1 <- -3*(Lens.Strength==6)-1*(Lens.Strength==18)+1*(Lens.Strength==36)+3*(Lens.Strength==60)
LS2 <- 1*(Lens.Strength==6)-1*(Lens.Strength==18)-1*(Lens.Strength==36)+1*(Lens.Strength==60)
LS3 <- -1*(Lens.Strength==6)+3*(Lens.Strength==18)-3*(Lens.Strength==36)+1*(Lens.Strength==60)
Eye1 <- 1*(Eye=="Left")-1*(Eye=="Right")

fit <- lme(Reac.Time ~ Sex1*Eye1*(LS1+LS2+LS3), random = ~1|Subject)
summary(fit)

# To check assumptions, look at residuals.

print(fit$resid)

qqnorm(fit$resid[,1])
plot(Sex,fit$resid[,1])
plot(Eye,fit$resid[,1])
plot(factor(Lens.Strength),fit$resid[,1])


# We can also add a random slope on Eye:

fit <- lme(Reac.Time ~ Sex1*Eye1*(LS1+LS2+LS3), random = ~(1+Eye1)|Subject)
summary(fit)