"plasma" <- 
matrix(c(4.2999999999999998, 3.7000000000000002, 4., 3.6000000000000001, 4.0999999999999996, 3.7999999999999998, 
	3.7999999999999998, 4.4000000000000004, 5., 3.7000000000000002, 3.7000000000000002, 4.4000000000000004,
	4.7000000000000002, 4.2999999999999998, 5., 4.5999999999999996, 4.2999999999999998, 3.1000000000000001,
	4.7999999999999998, 3.7000000000000002, 5.4000000000000004, 3., 4.9000000000000004, 4.7999999999999998,
	4.4000000000000004, 4.9000000000000004, 5.0999999999999996, 4.7999999999999998, 4.2000000000000002, 
	6.5999999999999996, 3.6000000000000001, 4.5, 4.5999999999999996, 3.2999999999999998, 2.6000000000000001,
	4.0999999999999996, 3., 3.7999999999999998, 2.2000000000000002, 3., 3.8999999999999999, 4., 
	3.1000000000000001, 2.6000000000000001, 3.7000000000000002, 3.1000000000000001, 3.2999999999999998, 
	4.9000000000000004, 4.4000000000000004, 3.8999999999999999, 3.1000000000000001, 5., 3.1000000000000001,
	4.7000000000000002, 2.5, 5., 4.2999999999999998, 4.2000000000000002, 4.2999999999999998, 4.0999999999999996,
	4.5999999999999996, 3.5, 6.0999999999999996, 3.3999999999999999, 4., 4.4000000000000004, 3., 
	2.6000000000000001, 3.1000000000000001, 2.2000000000000002, 2.1000000000000001, 2., 2.3999999999999999,
	2.7999999999999998, 3.3999999999999999, 2.8999999999999999, 2.6000000000000001, 3.1000000000000001, 
	3.2000000000000002, 3., 4.0999999999999996, 3.8999999999999999, 3.1000000000000001, 3.2999999999999998,
	2.8999999999999999, 3.2999999999999998, 3.8999999999999999, 2.2999999999999998, 4.0999999999999996, 
	4.7000000000000002, 4.2000000000000002, 4., 4.5999999999999996, 4.5999999999999996, 3.7999999999999998,
	5.2000000000000002, 3.1000000000000001, 3.7000000000000002, 3.7999999999999998, 2.6000000000000001, 
	1.8999999999999999, 2.2999999999999998, 2.7999999999999998, 3., 2.6000000000000001, 2.5, 2.1000000000000001,
	3.3999999999999999, 2.2000000000000002, 2.2999999999999998, 3.2000000000000002, 3.2999999999999998, 
	2.6000000000000001, 3.7000000000000002, 3.8999999999999999, 3.1000000000000001, 2.6000000000000001, 
	2.7999999999999998, 2.7999999999999998, 4.0999999999999996, 2.2000000000000002, 3.7000000000000002, 
	4.5999999999999996, 3.3999999999999999, 4., 4.0999999999999996, 4.4000000000000004, 3.6000000000000001,
	4.0999999999999996, 2.7999999999999998, 3.2999999999999998, 3.7999999999999998, 2.2000000000000002, 
	2.8999999999999999, 2.8999999999999999, 2.8999999999999999, 3.6000000000000001, 3.7999999999999998, 
	3.1000000000000001, 3.6000000000000001, 3.2999999999999998, 1.5, 2.8999999999999999, 3.7000000000000002,
	3.2000000000000002, 2.2000000000000002, 3.7000000000000002, 3.7000000000000002, 3.1000000000000001, 
	2.6000000000000001, 2.2000000000000002, 2.8999999999999999, 2.7999999999999998, 2.1000000000000001, 
	3.7000000000000002, 4.7000000000000002, 3.5, 3.2999999999999998, 3.3999999999999999, 4.0999999999999996,
	3.2999999999999998, 4.2999999999999998, 2.1000000000000001, 2.3999999999999999, 3.7999999999999998, 2.5,
	3.2000000000000002, 3.1000000000000001, 3.8999999999999999, 3.3999999999999999, 3.6000000000000001, 
	3.3999999999999999, 3.7999999999999998, 3.6000000000000001, 2.2999999999999998, 2.2000000000000002, 
	4.2999999999999998, 4.2000000000000002, 2.5, 4.0999999999999996, 4.2000000000000002, 3.1000000000000001,
	1.8999999999999999, 3.1000000000000001, 3.6000000000000001, 3.7000000000000002, 2.6000000000000001, 
	4.0999999999999996, 3.7000000000000002, 3.3999999999999999, 4.0999999999999996, 4.2000000000000002, 4.,
	3.1000000000000001, 3.7999999999999998, 2.3999999999999999, 2.2999999999999998, 3.6000000000000001, 
	3.3999999999999999, 3.1000000000000001, 3.8999999999999999, 3.7999999999999998, 3.6000000000000001, 3.,
	3.5, 4., 4., 2.7000000000000002, 3.1000000000000001, 3.8999999999999999, 3.7000000000000002, 
	2.3999999999999999, 4.7000000000000002, 4.7999999999999998, 3.6000000000000001, 2.2999999999999998, 3.5,
	4.2999999999999998, 3.5, 3.2000000000000002, 4.7000000000000002, 3.6000000000000001, 3.7999999999999998,
	4.2000000000000002, 4.4000000000000004, 3.7999999999999998, 3.5, 4.2000000000000002, 2.5, 
	3.1000000000000001, 3.7999999999999998, 4.4000000000000004, 3.8999999999999999, 4., 4., 3.7000000000000002,
	3.5, 3.7000000000000002, 3.8999999999999999, 4.2999999999999998, 2.7999999999999998, 3.8999999999999999,
	4.7999999999999998, 4.2999999999999998, 3.3999999999999999, 4.9000000000000004, 5., 4., 2.7000000000000002,
	3.6000000000000001, 4.4000000000000004, 3.7000000000000002, 3.5, 4.9000000000000004, 3.8999999999999999,
	4., 4.2999999999999998, 4.9000000000000004, 3.7999999999999998, 3.8999999999999999, 4.7999999999999998,
	3.5, 3.2999999999999998, 3.7999999999999998)
, nrow = 33, ncol = 8
,  dimnames = list(character(0)
, c("T0.0", "T0.5", "T1.0", "T1.5", "T2.0", "T2.5", "T3.0", "T4.0")
)
)

library(nlme)

# 
# Plotting the raw phosphate values with line plots for each subject.
# Separate plots are given for the 'control' group and the 'obese' group.
# 
par(mfrow=c(1,2)) 
matplot(matrix(c(0.0,0.5,1.0,1.5,2.0,3.0,4.0,5.0),ncol=1), 
t(plasma[1:13,]),type="l",col=1,lty=1, 
xlab="Time (hours after glucose challenge)",ylab="Plasma inorganic phosphate",ylim=c(1,7)) 
title("Control") 
matplot(matrix(c(0.0,0.5,1.0,1.5,2.0,3.0,4.0,5.0),ncol=1), 
t(plasma[14:33,]),type="l",col=1,lty=1, 
xlab="Time (hours after glucose challenge)",ylab="Plasma inorganic phosphate",ylim=c(1,7)) 
title("Obese") 
# 
# We see a quadratic-type trend across time in each group.
# The curve shapes are somewhat different between the groups:  
# Is there a time-by-group interaction?
#

#
# Is the association between response values at two timepoints
# the same no matter how spaced the times are apart?
#
pairs(plasma[1:13,])  # control subjects
pairs(plasma[14:33,]) # obese subjects
# 
# It doesn't look that way.
# We should not assume "compound symmetry".
# 

#
# Putting the data matrix in "long" form 
#
group<-rep(c(0,1),c(104,160)) 
# 
time<-c(0.0,0.5,1.0,1.5,2.0,3.0,4.0,5.0) 
time<-rep(time,33) 
# 
subject<-rep(1:33,rep(8,33)) 
plasma.dat<-cbind(subject,time,group,as.vector(t(plasma))) 
dimnames(plasma.dat)<-list(NULL,c("Subject","Time","Group","Plasma")) 
# 
plasma.df<-data.frame(plasma.dat) 
plasma.df$Group<-factor(plasma.df$Group,levels=c(0,1),labels=c("Control","Obese")) 
# 
attach(plasma.df) 
#
# 'plasma.df' is the data frame in "long" form.


# 
# A mixed model with a quadratic time effect:
#
Time.sq <- Time^2
plasma.lme1<-lme(Plasma~Time+Time.sq+Group,random=~Time|Subject, 
data=plasma.df,method="ML") 
summary(plasma.lme1) 
# 
# This model assumes the effect of Time on phosphate
# is the same for each group (no Time-by-Group interaction).


# Fitting a an ordinary regression model that 
# ignores the correlation between the repeated measurements:
#
#multiple regression model-independence model
summary(lm(Plasma~Time+Time.sq+Group,data=plasma.df)) 
# 
# Notice the deceptively small standard error for
# the Group effect, which leads to erroneous significance


#
# fitted values for the subjects (separated by group)
# based on the above mixed model:
#
predictions<-matrix(predict(plasma.lme1),ncol=8,byrow=T) 
par(mfrow=c(1,2)) 
matplot(matrix(c(0.0,0.5,1,1.5,2,3,4,5),ncol=1), 
t(predictions[1:13,]),type="l",lty=1,col=1, 
xlab="Time (hours after glucose challenge)",ylab="Plasma inorganic phosphate",ylim=c(1,7)) 
title("Control") 
matplot(matrix(c(0.0,0.5,1,1.5,2,3,4,5),ncol=1), 
t(predictions[14:33,]),type="l",lty=1,col=1, 
xlab="Time (hours after glucose challenge)",ylab="Plasma inorganic phosphate",ylim=c(1,7)) 
title("Obese") 
# 


# 
# A mixed model with a quadratic time effect
# and a Time-by-Group interaction:
#
plasma.lme2<-lme(Plasma~Time*Group+Time.sq,random=~Time|Subject, 
data=plasma.df,method="ML") 
summary(plasma.lme2) 
#
# This model assumes the effect of Time on phosphate
# may be different for each group (Time-by-Group interaction).

# Which mixed model is preferred?
#
anova(plasma.lme1,plasma.lme2) 
# 

#
# fitted values for the subjects (separated by group)
# based on the mixed model with interaction:
#
predictions<-matrix(predict(plasma.lme2),ncol=8,byrow=T) 
par(mfrow=c(1,2)) 
matplot(matrix(c(0.0,0.5,1,1.5,2,3,4,5),ncol=1), 
t(predictions[1:13,]),type="l",lty=1,col=1, 
xlab="Time (hours after glucose challenge)",ylab="Plasma inorganic phosphate",ylim=c(1,7)) 
title("Control") 
matplot(matrix(c(0.0,0.5,1,1.5,2,3,4,5),ncol=1), 
t(predictions[14:33,]),type="l",lty=1,col=1, 
xlab="Time (hours after glucose challenge)",ylab="Plasma inorganic phosphate",ylim=c(1,7)) 
title("Obese") 

#
#
# Checking the normality assumptions for the random effects 
# and for the random error term:
# 
res.int<-random.effects(plasma.lme2)[,1] 
res.int 
res.slope<-random.effects(plasma.lme2)[,2] 
par(mfrow=c(1,3)) 
qqnorm(res.int,ylab="Estimated random intercepts",main="Random intercepts") 
qqnorm(res.slope,ylab="Estimated random slopes",main="Random slopes") 
resids<-resid(plasma.lme2) 
qqnorm(resids,ylab="Estimated residuals",main="Residuals") 
#