########## R script: pigweights ##########

# For doing a longitudinal data analysis of the
# Victorian pig weights data.

# Last changed: 02 APR 2007

# Define function wait() for sending prompts
# to the screen. (Very useful for scripts such
# as this).

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

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

library(nlme)

# Read in data into data frame named 'pigweights'

pigweights <- read.table("pigweights.txt",header=T)

# Convert `pigweights' to a `grouped data object'

pigweights <- groupedData(weight~weeks|idnum,data=pigweights)

# Plot the data using a `trellis' plot.

print(plot(pigweights))
wait()

# Fit random intercept model

fit1 <- lme(weight~weeks,random=~1,data=pigweights)

# Print estimates and 95% confidence intervals for the
# main parameters.

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

# Fit random intercept & slope model

fit2 <- lme(weight~weeks,data=pigweights)

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

# Compute likelihood ratio statistic for 
# comparing two models

LR.statist <- 2*(fit2$logLik - fit1$logLik)

cat("\n\n\n\n\n\n")
cat("Likelihood ratio statistic for comparing\n")
cat("random intercept (null) and random intercept & slope \n")
cat("(alternative) models is ")
print(LR.statist)
cat("\n\n\n\n\n\n")

print(plot(fit1,main="Random intercept model"))
wait()

print(plot(fit2,main="Random intercept & slope model"))

########## End of pigweights ##########