#=================================
# Random intercept and slope model
# Example: FEV data
#=================================

# Read FEV data from the web
data.name="http://www.stat.ufl.edu/~aa/glm/data/FEV2.dat"
fev=read.table(data.name,header=T)
head(fev)

# Plot data (using quite "primitive" commands):
plot(1:8,fev[fev$patient==1,"fev"],type='l',ylim=c(0,5),
     col='red',ylab="FEV",xlab="Hours")
for(i in 2:24) lines(1:8,fev[fev$patient==i,"fev"],col='red')
for(i in 25:48) lines(1:8,fev[fev$patient==i,"fev"],col='blue')
for(i in 49:72) lines(1:8,fev[fev$patient==i,"fev"])

# Define placebo as reference group for treatment:
fev$drug=relevel(as.factor(fev$drug),ref="p")

# Use the nlme library 
library(nlme)

# Fit for illustration model with random intercept 
fit.fev.1=lme(fev~base+factor(drug)+hour,random=~1|patient, data=fev) 
summary(fit.fev.1)

# Fit for illustration model with random intercept and slope 
fit.fev.2=lme(fev~base+factor(drug)+hour, random=~1+hour|patient, data=fev) 
summary(fit.fev.2)