library(nlme)

OrthoMF <- Orthodont
OrthoMF$Sex1[Orthodont$Sex =="Female"] <- 1
OrthoMF$Sex1[Orthodont$Sex == "Male"] <- 0

#or use ifelse, Thanks John!

#What's the difference between the fits with Sex and Sex1?

fit1 <- lme(distance~Sex1*I(age-11), data=OrthoMF, random=~1 | Subject, method="REML")
summary(fit1)

fit3 <- lme(distance~Sex*I(age-11), data=OrthoMF, random=~1 | Subject, method="REML")
summary(fit3)

#Fit with both F and M

fit1b <- lme(distance ~ Sex1 * I(age - 11), data = OrthoMF, 
    random = ~1 | Subject, weights = varIdent(form = ~1 | Sex1))

plot(augPred(fit1b))

#Return to diagnostics
qqnorm(resid(fit1b))
qqline(resid(fit1b))
plot(fit1b, resid(.,type="p")~fitted(.) | Sex1, abline=0, lty=2)


#Thanks, Peter!
#What do you notice about younger girls?

plot(NA, xlim=c(8,14), ylim=c(15,35), xlab="Age", ylab="Fitted Values")

for(i in 1:length(levels(Orthodont$Subject))){
lines(Orthodont$age[(4*i-3):(4*i)], fitted(fit1b)[(4*i-3):(4*i)], col=2+2*(Orthodont$Sex[4*i-3]=="Male"))
text(Orthodont$age[4*i-3*(Orthodont$Sex[4*i-3]=="Male")], fitted(fit1b)[4*i-3*(Orthodont$Sex[4*i-3]=="Male")], 
Orthodont$Subject[4*i-3], cex=0.6)
}
legend(9, 35, c("Males", "Females"), col=c(4,2), lty=1, bty="n", horiz=T)

#Now let's pretend that we know even more about these Subjects ... 

hold1 <- c(rep("A",32),rep("B",24),rep("C",8))
hold2 <- c(rep("A",24),rep("B",16),rep("C",4))

grp <- c(hold1, hold2)

#Make a grp a factor for ANOVA
OrthoMF$grp <- as.factor(grp)

fit1c <- lme(distance~Sex1*I(age-11) + grp, data=OrthoMF, random=~1 | Subject, method="REML")
summary(fit1c)

#What about this?

hold3 <- c(rep(1,32),rep(2,24),rep(3,8))
hold4 <- c(rep(1,24),rep(2,16),rep(3,4))

grp2 <- c(hold3, hold4)

#grp2 is treated as a continuous variable 1,2,3,4!
OrthoMF$grp2 <- grp2

fit1d <- lme(distance~Sex1*I(age-11) + grp2, data=OrthoMF, random=~1 | Subject, method="REML")
summary(fit1d)