#=============================
# Prediction of random effects
# Example: Fecal fat
#=============================

# Read data on fecal fat: 
data.name="http://www.uio.no/studier/emner/matnat/math/STK3100/data/fecfat.txt"
fecfat=read.table(data.name,header=T)
head(fecfat)

# Plot data:

treat.mat=matrix(1:4,nrow=4,ncol=6)
fecfat.mat=matrix(fecfat$fecfat,nrow=4)
matplot(treat.mat,fecfat.mat,type='l',lty=1,col=1:6,
        xlab="Treatment",ylab="Fecal Fat",xaxt='n')
axis(1,at=1:4,labels=c("None","Tablet","Capsule","Coated"))
legend(2,71,lty=1,col=1:6,legend=paste(rep("Patient",6),1:6),cex=0.5)
# Use the nlme-library 
library(nlme)

# Fit model with random effect for patient
fit.fec=lme(fecfat~factor(pilltype), random=~1|subject,data=fecfat)
summary(fit.fec)

# Predict and plot random effects
ranef.fec=ranef(fit.fec)
ranef.fec
plot(ranef.fec)