########## R script: Orthodont ##########

# For doing a longitudinal data analysis of the
# orthodontic data.

# Last changed: 02 APR 2007

# Define function wait() for sending prompts
# to the screen. 

wait <- function()
{
        cat("Hit return to continue\n")
        ans <- readline()
        invisible()
}

# Load in the `nonlinear linear mixed effects' models
# library of R functions.

library(nlme)

# Load in the `Orthodont' data set (this comes with `nlme')
# The response variable is:
#
# distance = distances from the pituitary to the
#             pterygomaxillary fissure (mm)

data(Orthodont)

# Plot the data

library(lattice)
pobj <-
xyplot(distance ~ age |Sex, data =Orthodont , groups =Subject,
        strip = function(...) strip.default(..., style = 1),
        subscripts = T, panel = function(x, y, subscripts,groups) 
        {
           panel.grid()
           panel.xyplot(x, y,pch=16,col="black")
           panel.superpose(x, y, subscripts, groups, type = "l")
        }
        ,xlab = "age (years)", ylab ="distance",
         bg="white")
         
print(pobj)
wait()

# All subjects seem to have roughly the same slope, so
# fit a random intercept model with and indicator for
# `Sex'.

fit <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)

cat("\n\n\n\n\n\n")
print(intervals(fit))
cat("\n\n\n\n\n\n")
wait()

# Examine the residuals

print(plot(fit))
 
########## End of Orthodont ##########