library(nlme)
library(lattice)

data(Orthodont)
help(Orthodont)
head(Orthodont,n=20)
attach(Orthodont)

xyplot(distance~age|Subject,
	panel=function(x,y){
		panel.lmline(x,y,type="l")
		panel.xyplot(x,y,col=1,pch=16)})

age.means <- tapply(distance,age,mean)
plot(c(8,14),c(20,28),xlab="Age",ylab="Mean distance",pch=16,type="n")
points(c(8,10,12,14),age.means,type="b",col=1,pch=16)

I <- Sex=="Male"
age.means.m <- tapply(distance[I],age[I],mean)
points(c(8,10,12,14),age.means.m,pch=16,col=4)
I <- Sex=="Female"
age.means.m <- tapply(distance[I],age[I],mean)
points(c(8,10,12,14),age.means.m,pch=16,col=3)
legend("topleft",legend=c("Males","Females","Combined"),pch=16,col=c(4,3,1))


lmList.obj <- lmList(distance~age|Subject,data=Orthodont) # fits a linear regression model with different slopes and intercepts
										# for each child
coef(lmList.obj)
plot(intervals(lmList.obj))
plot(coef(lmList.obj)[,2],(coef(lmList.obj)[,1]),xlab= "Estimated slope",ylab="Estimated intercept") 

lmList.obj <- lmList(distance~I(age-11)|Subject,data=Orthodont)
coef(lmList.obj)
plot(intervals(lmList.obj))

lm.obj<-lm(distance~age+Sex,data=Orthodont)
summary(lm.obj)
anova(lm.obj)

lme.obj <- lme(distance~I(age-11)+Sex,data=Orthodont,random=~1|Subject)
lme.obj$coefficients  # gives estimated fixed and predicted random effects
lme.obj$contrasts # shows aliasing and/or contrast

# Next line : shows the fixed effect model only predictions (column 3) 
#	and the mixed effects model predictions (column 4)
cbind(Subject,age,lme.obj$fitted)  

summary(lme.obj)
intervals(lme.obj)

re <- random.effects(lme.obj)
hist(unlist(re))
qqnorm(unlist(re))

library(lattice)
xyplot(lme.obj$fitted[,1]~age,groups=Subject,col=1,ylab="Fixed effects prediction",pch=16)
xyplot(lme.obj$fitted[,2]~age,groups=Subject,col=1:50,ylab="Full LMM prediction")
xyplot(lme.obj$fitted[,2]~age|Subject,
	panel=function(x,y){
		panel.lmline(x,y,type="l")
		panel.xyplot(x,y,col=1,pch=16)})

lme.obj <- lme(distance~age+Sex,data=Orthodont,random=~1+age|Subject)
summary(lme.obj)
anova(lme.obj)

xyplot(lme.obj$fitted[,1]~age,groups=Subject,col=1,ylab="Fixed effects prediction",pch=16)
xyplot(lme.obj$fitted[,2]~age,groups=Subject,col=1:50,ylab="Full LMM prediction")
xyplot(lme.obj$fitted[,2]~age|Subject,ylab="Full LMM prediction")
xyplot(lme.obj$fitted[,2]~age|Subject,
	panel=function(x,y){
		panel.lmline(x,y,type="l")
		panel.xyplot(x,y,col=1,pch=16)})

lme.obj.1 <- lme(distance~I(age-11)+Sex,data=Orthodont,random=~1|Subject,method="ML")
lme.obj.2 <- lme(distance~I(age-11)*Sex,data=Orthodont,random=~1|Subject,method="ML")
summary(lme.obj.2)
anova(lme.obj.1,lme.obj.2)

############################# variance structures

lme.obj <- lme(distance~I(age-11)+Sex,data=Orthodont,random=pdDiag(~1+age))
summary(lme.obj)

lme.obj <- lme(distance~I(age-11)+Sex,data=Orthodont,random=pdSymm(~1+age))
summary(lme.obj)