#LMM 01 ONE-WAY RANDOM ANOVA
library(nlme)      # {nlme} for lme() & other functions
library(ape)       # {ape} for varcomp()
library(help=nlme) # prototype for finding package index

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#PINHEIRO & BATES MIXED-EFFECTS MODELS 

#RAIL EXAMPLE LINEAR APPROACH
data(Rail)
plot(Rail)
Rail                   #GROUPED DATA OBJECT
rail=data.frame(Rail)
rail                   #CORRESPONDING DATA.FRAME

FM=lm(travel~Rail, data=rail)
summary(FM)
anova(FM)

FM2=lm(travel~Rail, data=Rail)
summary(FM2) 
anova(FM2)             #lm() WORKS WITH EITHER DATA TYPE

Y=rail$travel
X=factor(rail$Rail,ordered=F) #CONVERTING TO UNORDERED FACTORS 
X                             #WITH TREATMENTS CONTRASTS
FM3=lm(Y~X)
summary(FM3)      #RAIL 2 IS "INTERCEPT" IN THE REPORT
model.matrix(FM3) #NOTE DIFFERENCE IN COLUMN ARRANGEMENT FOR FACTORS
anova(FM3)

FM4=lm(Y~X-1)     #CELL MEANS MODEL CONTRASTS
model.matrix(FM4)
anova(FM4)        #NOTE DIFFERENT DF FOR MODEL HERE FM4 VS FM3

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#lmList() LINEAR MODELS FOR GROUPS
setwd("c:/DATA/Models")
K=read.table("Railframe.txt")
K

FM5=lm(travel~factor(Rail)-1,data=K)  #FM5 = CELL MEANS MODEL ANOVA
summary(FM5)
anova(FM5)  #NOTE 6 DF FOR Rail here as opposed to 5 df in FM5 

lmList(travel~1|Rail,data=K) #SEPARATE REGRESSION COEFFICIENTS FOR EACH GROUP
             #GIVES CELL MEANS MODEL RESULT

railfactor=factor(K$Rail)
contrasts(railfactor)=contr.sum #CONVERTS FACTOR CONTRASTS
FM6=lm(travel~railfactor,data=K)
summary(FM6)
anova(FM6)  #df 5 for railfactor IN contr.sum MODEL ANOVA IN R

lmList(Rail) #SEPARATE REGRESSION COEFFICIENTS FOR EACH GROUP
             #GIVES CELL MEANS MODEL RESULT

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#lme() MIXED-EFFECTS MODELS APPROACH

FMe=lme(travel~1, data=Rail, random=~1|Rail)
summary(FMe)
anova(FMe)

coef(FM3)    # LEAST SQUARES COEFFICIENTS treatments MODEL IN R
lmList(Rail) # LEAST SQUARES COEFFICIENTS cell means MODEL IN R
coef(FMe)    # REML COEFFICIENTS OF FIXED EFFECTS  
ranef(FMe)   # REML COEFFICIENTS OF RANDOM EFFECTS (i.e, BLUP's)

              #VARIANCE COMPONENTS:
VarCorr(FMe)  #SQUARE OF StdDev IN Summary.lme
varcomp(FMe)  #Alternate Function in {ape}
intervals(FMe) #APPROXIMATE CONFIDENCE INTERVALS

predict(FMe)

R=Rail$Rail
T=Rail$travel
PD=predict(FM)
PE=predict(FMe)

plot(T,R,col="black")              #DATA POINTS
points(PD,R,col="red", pch=20)     #PREDICTIONS FROM lm()
points(PE,R,col="blue", pch="+")   #PREDICTIONS FROM lme()