# This file posted as  milkprotein.R
#
# This code is applied to the milk protein 
# data from DHLZ, page 8.

set1 <- read.table("c:/documents and settings/kkoehler/IASTATE/my documents/courses/st565/data/dhlz.example1_4.data", 
       col.names=c("diet","cow","week","protein"))
set1

 
# Create factors

set1$dietf <- as.factor(set1$diet)
set1$weekf <- as.factor(set1$week)
set1$cowf <- as.factor(set1$cow)

#  Sort the data set by subject

i <- order(set1$cow,set1$week)
set1 <- set1[i,]


# Delete the list called i         
rm(i)  


# Compute sample means

means <- tapply(set1$protein,
    list(set1$week,set1$diet),mean)
means


#  Make a profile plot of the means
#  Unix users should insert the motif( )
#  command

x.axis <- unique(set1$week)

par(fin=c(6.0,6.0),pch=18,mkh=.1,mex=1.5,
            cex=1.2,lwd=3)
matplot(c(1,19), c(3.0,4.0), type="n", 
     xlab="Time(weeks)", ylab="Protein (percent)",
     main= "Observed Protein Means")   
matlines(x.axis,means,type='l',lty=c(1,3,7))
matpoints(x.axis,means, pch=c(16,17,15))    
legend(1,2.45,legend=c("Barley diet",
    'Barley+lupins','Lupin diet'),lty=c(1,3,7),bty='n')


#  Use the gls( ) function to fit a
#  model where the errors have a 
#  compound symmetry covariance structure
#  within cows. 

options(contrasts=c("contr.treatment","contr.poly"))

set1.glscs <- gls(protein ~ dietf+ week+
     dietf*week+week^2+dietf*week^2 +
     week^3 + dietf*week^3,
  data=set1,
  correlation = corCompSymm(form=~1|cow),
  method=c("REML"))
summary(set1.glscs)
anova(set1.glscs)



# Try an auto regressive covariance 
# structures across weeks within cows

set1.glsar <- gls(protein ~ dietf+ week+
     dietf*week+week^2+dietf*week^2 +
     week^3 + dietf*week^3,
  data=set1,
  correlation = corAR1(form=~1|cow),
  method=c("REML"))
summary(set1.glsar)
anova(set1.glsar)


#  Try a general correlation structure

set1.glss <- gls(protein ~ dietf + week+
     dietf*week +week^2+dietf*week^2 + 
      week^3 + dietf*week^3,
  data=set1,
  correlation = corSymm(form=~1|cow),
  method=c("REML"))
summary(set1.glss)
anova(set1.glss)


#  Try an AR(1) correlation structure
#  with hterogeneous variances

set1.glss <- gls(protein ~ dietf + week+
     dietf*week +week^2+dietf*week^2 + 
      week^3 + dietf*week^3,
  data=set1,
  correlation = corAR1(form=~1|cow),
  weights = varIdent(form = ~ 1 | week),
  method=c("REML"))
summary(set1.glsarh)
anova(set1.glss)


#  Compare the fit of various covariance
#  structures.

anova(set1.glss, set1.glscs)
anova(set1.glss, set1.glsar)
anova(set1.glss, set1.glsar, set1.glsarh)


# To compare the continuous week model to the 
# model where we fit a different mean at each 
# time point, we must compare likelihood values
# instead of REML likelihood values.


set1.glsarmle <- gls(protein ~ dietf+
     weekf+dietf*weekf,
  data=set1,
  correlation = corAR1(form=~1|cow),
  method=c("ML"))

set1.glscarmle <- gls(protein ~ dietf+ 
     week+ dietf*week + week^2 + 
     dietf*week^2 + week^3 + dietf*week^3,
  data=set1,
  correlation = corAR1(form=~1|cow),
  method=c("ML"))


anova(set1.glsarmle,set1.glscarmle)