#=====================================================
# Systematic model selection (cf. slide 2 to week 46)
#=====================================================
# 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)

# Steps 1 and 2 (cf. slide 2)

# Choice of covariance structure using REML with a "wide model" 
# for the fixed effects (more models for the covariance 
# structure could have been considered)
#(NB: Default "method" in lme is "REML")
fit.fev.3=lme(fev~base*factor(drug)*hour,random=~1|patient,data=fev)
fit.fev.4=lme(fev~base*factor(drug)*hour,random=~1+hour|patient,data=fev)
fit.fev.5=lme(fev~base*factor(drug)*hour,random=~1|patient, correlation=corAR1(),data=fev)
fit.fev.6=lme(fev~base*factor(drug)*hour,random=~1+hour|patient,correlation=corAR1(),data=fev)

anova(fit.fev.3,fit.fev.4,fit.fev.6)
anova(fit.fev.3,fit.fev.5,fit.fev.6)

# The covariance structure given as in model "fit.fev.6" gives 
# the best fit to the data (of the models considered here!), 
# so we will continue to use this

# Step 3
# Choice of fixed effects structure using ML-estimation 
# with the covariance structure determined above
#NB: Default method in lme() is REML

#All two-way interactions and the three-way interaction
# (NB default ''method'' value is "REML")
fit.fev.7=lme(fev~base*factor(drug)*hour, random=~1+hour|patient,
              correlation=corAR1(),data=fev,method="ML")

#Should we drop the three-way interaction?
drop1(fit.fev.7,scope=~base*factor(drug)*hour,data=fev,test="Chisq")
fit.fev.8=lme(fev~base+factor(drug)+hour+base:factor(drug)+base:hour+factor(drug):hour, 
              random=~1+hour|patient, correlation=corAR1(),data=fev,method="ML")

#Should we drop any two-way interactions?
drop1(fit.fev.8,scope=~base+factor(drug)+hour+base:factor(drug)+base:hour+factor(drug):hour,data=fev,test="Chisq")

fit.fev.9=lme(fev~base+factor(drug)+hour+factor(drug):hour+base:hour, 
              random=~1+hour|patient, correlation=corAR1(),data=fev,method="ML")
#Should we drop more two-way interactions?
drop1(fit.fev.9,scope=~base+factor(drug)+hour+factor(drug):hour+base:hour,data=fev,test="Chisq")
fit.fev.10=lme(fev~base+factor(drug)+hour+factor(drug):hour, 
               random=~1+hour|patient, correlation=corAR1(),data=fev,method="ML")
drop1(fit.fev.10,scope=~base+factor(drug)+hour+factor(drug):hour,data=fev,test="Chisq")

# The interaction factor(drug):hour seems important.

# Step 4
# "Final model" fitted by REML:
fit.fev.final=lme(fev~base+factor(drug)+hour+factor(drug):hour, 
                  random=~1+hour|patient, correlation=corAR1(),data=fev)
summary(fit.fev.final)

#=================================
# Prediction of random effects
# Example: FEV data
#=================================

# Above we found the "final model" fit.fev.final
# Predict and plot the random effects
ranef.fev=ranef(fit.fev.final)
ranef.fev
plot(ranef.fev)
plot(ranef.fev[,1],ranef.fev[,2],xlab="u1",ylab="u2")

#===========================
# Marginal linear model
#=========================


# Then consider marginal linear models

# Reorder the rows of the dataframe so that the observations in a cluster 
# are contiguous rows (a requirement for the gee-command below)
ind=order(fev$patient)
fev=fev[ind,] 

# Use the gee-package
library(gee)

# Fit marginal model assuming "working independence":
fit.indep=gee(fev~base+factor(drug)+hour+factor(drug):hour, id=patient,corstr="independence", data=fev)
summary(fit.indep)

# Fit marginal model assuming "exchangeable correlation":
fit.exch=gee(fev~base+factor(drug)+hour+factor(drug):hour, id=patient,corstr="exchangeable", data=fev)
summary(fit.exch)

# Fit marginal model assuming "AR-1 correlation":
fit.ar1=gee(fev~base+factor(drug)+hour+factor(drug):hour, id=patient,corstr="AR-M", Mv=1, data=fev)
summary(fit.ar1)