setwd("D:/Courses/Longitudinal Data Analysis/Data Sets/")

diet.data <- read.table("Vit_E_CH27.txt",header=T,na.strings = "-99")
attach(diet.data)

# Suppose we want to look at the full picture of how weight increases.

plot(Week[Animal==1],Weight[Animal==1],type="l",xlab="Week",ylab="Weight",
  col=Group[(Animal==1)&(Week==1)],ylim=c(min(Weight),max(Weight)))
for (i in 2:15){
  lines(Week[Animal==i],Weight[Animal==i],type="l",col=Group[(Animal==i)&(Week==1)])
}

# A natural summary is the average weight, by week, in each group.

plot(c(1,3,4,5,6,7),tapply(Weight[Group==1],Week[Group==1],mean),type="l",
  col=1,ylim=c(min(Weight),max(Weight)))
lines(c(1,3,4,5,6,7),tapply(Weight[Group==2],Week[Group==2],mean),col=2)
lines(c(1,3,4,5,6,7),tapply(Weight[Group==3],Week[Group==3],mean),col=3)

# For a more robust summary, use the median weight, by week, in each group.

plot(c(1,3,4,5,6,7),tapply(Weight[Group==1],Week[Group==1],median),type="l",
  col=1,ylim=c(min(Weight),max(Weight)))
lines(c(1,3,4,5,6,7),tapply(Weight[Group==2],Week[Group==2],median),col=2)
lines(c(1,3,4,5,6,7),tapply(Weight[Group==3],Week[Group==3],median),col=3)

# We might also want to see how the variability changes across time.

plot(c(1,3,4,5,6,7),tapply(Weight[Group==1],Week[Group==1],sd),type="l",
  col=1)
lines(c(1,3,4,5,6,7),tapply(Weight[Group==2],Week[Group==2],sd),col=2)
lines(c(1,3,4,5,6,7),tapply(Weight[Group==3],Week[Group==3],sd),col=3)


# What is the extent of correlation in weight by weeks?

cor(cbind(Weight[Week==1],Weight[Week==3],Weight[Week==4],
          Weight[Week==5],Weight[Week==6],Weight[Week==7]))

# Review of ANOVA via a comparison of the weights at weeks 1 and 7 and
# the change in weight.

plot(Group[Week==1],Weight[Week==1],xlab="Group",ylab="Weight")
tapply(Weight[Week==1],Group[Week==1],mean)
tapply(Weight[Week==1],Group[Week==1],sd)
fit <- aov(Weight[Week==1] ~ factor(Group[Week==1]))
summary(fit)
plot(fit)

plot(Group[Week==7],Weight[Week==7],xlab="Group",ylab="Weight")
tapply(Weight[Week==7],Group[Week==7],mean)
tapply(Weight[Week==7],Group[Week==7],sd)
fit <- aov(Weight[Week==7] ~ factor(Group[Week==7]))
summary(fit)
plot(fit)

plot(Group[Week==7],Weight[Week==7]-Weight[Week==1],xlab="Group",ylab="Weight")
tapply(Weight[Week==7]-Weight[Week==1],Group[Week==7],mean)
tapply(Weight[Week==7]-Weight[Week==1],Group[Week==7],sd)
sel <- ((Week==1)|(Week==7))
plot(Week[(Animal==1)&(sel==T)],Weight[(Animal==1)&(sel==T)],type="l",
  xlab="Week",ylab="Weight",
  col=Group[(Animal==1)&(Week==1)],ylim=c(min(Weight),max(Weight)))
for (i in 2:15){
  lines(Week[(Animal==i)&(sel==T)],Weight[(Animal==i)&(sel==T)],type="l",
    col=Group[(Animal==i)&(Week==1)])
}
fit <- aov(Weight[Week==7]-Weight[Week==1] ~ factor(Group[Week==7]))
summary(fit)
plot(fit)

# Review of regression via the relationship of weight at week 7 to
# weight at week 1.

plot(Weight[Week==1],Weight[Week==7],xlab="Initial Weight",ylab="Final Weight")

fit <- lm(Weight[Week==7] ~ Weight[Week==1])
summary(fit)
plot(fit)

# Analysis of covariance for weight at week 7.

fit <- lm(Weight[Week==7] ~ Weight[Week==1]+factor(Group[Week==1]))
summary(fit)
plot(fit)

aov(fit)