#Exercise 22 # In this exercise we will look at how a linear mixed model can be used to analyze the dataset "Orthodont"

library(nlme); data(Orthodont)
#Y=distance, X1=age, X2=Sex. The dataset is grouped after the variable Subject. 
#> head(Orthodont)
# distance age Subject  Sex
#     26.0   8     M01 Male
#     25.0  10     M01 Male
#     29.0  12     M01 Male
#     31.0  14     M01 Male
#     21.5   8     M02 Male
#     22.5  10     M02 Male


#A: Transform Subject to the values 1-27
Orthodont$ID = as.factor(as.numeric(Orthodont$Subject))

#B: Exploratory analysis
plot(Orthodont)     #Plot Age againts Distance for each Subject
s = split(Orthodont$distance, Orthodont$age)
boxplot(s[1:4])

#C: Plot distance against age. For every group, add a regression line
plot(Orthodont$age, Orthodont$distance,col=Orthodont$ID)
for(i in 1:27){
    od = Orthodont[Orthodont$ID==i,]
    abline(lm(distance~age,data=od),col=i)
}

#D, perform a Linear Mixed Effects model:
fit1 = lme(distance ~ age,data=Orthodont,random=~1|ID)

#plot the data and the results from the LME (small clarification of the code)
plot(Orthodont$age,Orthodont$distance,col=Orthodont$ID)
alpha_hat = fit1$coef$fixed[1]
beta_hat  = fit1$coef$fixed[2]
abline(a=alpha_hat, b=beta_hat,lwd=4)
for (i in 1:27){
    b_hat_i = fit1$coef$random$ID[i,]
    abline(a=alpha_hat+b_hat_i, b=beta_hat,col=i)
}


#E
summary(fit1)
#Random effects:
# Formula: ~1 | ID
#        (Intercept) Residual
#StdDev:    2.114724 1.431592

#	d^2		= 2.114772^2 = 4.47
#	Sigma^2 = 1.431592^2 = 2.05


#F: Add the fixed covariate "sex".
fit2 = lme(distance ~ age+Sex,data=Orthodont,random=~1|ID)

#Is it a significant covariate based on the Wald test?
summary(fit2)
#-2.321023/0.7614168 = -3.048295, p < 0.05. SexFemale is significant and important in the model.