#LMM 060 PROTOCOL FOR LINEAR MIXED MODELS
library(nlme)
setwd("c:/DATA/Models")
Owls=read.table("owls.txt",header=T)
Owls
#LOG TRANSFORMATION OF RESPONSE VARIABLE TO MATCH ZUUR
Owls$LogNeg=log10(Owls$NegPerChick+1) 

# Step 1A - FITTING THE FULL LINEAR MODEL FIXED PART
FMform=formula(LogNeg~SexParent*FoodTreatment+SexParent*ArrivalTime)
FM=lm(FMform,data=Owls)
summary(FM)

#STANDARD VALIDATION PLOTS
op=par(mfrow = c(2, 2))
plot(FM)
par(op)

#HOMEMADE PLOTS
op=par(mfrow = c(3,2)) 
plot(Owls$SexParent, resid(FM),pch=19,col='green',xlab="SexParent",ylab="Residuals")
plot(Owls$FoodTreatment, resid(FM),col='red',xlab="FoodTreatment",ylab="Residuals",pch=20)
plot(Owls$ArrivalTime, resid(FM),col='blue',xlab="ArrivalTime",ylab="Residuals",pch=20)
hist(resid(FM),col="brown",main='',xlab="Residuals")
boxplot(resid(FM)~Nest,data=Owls,col='green',ylab='Residuals')
boxplot(rstandard(FM)~Nest,data=Owls,col='red',ylab='Standardized Residuals')
par(op)

#STEP 1B - SEPARATE LINEAR MODELS FOR EACH NEST
#NOTE: SexParent doesn't have two levels for each nest 
#      so must be dropped from formula here
FMLform=formula(LogNeg~FoodTreatment+ArrivalTime|Nest)
OG=groupedData(FMLform,data=Owls)
OG
FMG=lmList(OG)
plot(FMG)
plot(getGroups(FMG),resid(FMG),col="red",ylab='Residuals')
plot(intervals(FMG))

#STEP 2 - FITTING THE FULL MODEL FIXED PART WITH gls()
FMgls=gls(FMform,data=Owls)
summary(FMgls)

#STEPS 3-5 - CHOOSE A VARIANCE STRUCTURE using REML
FMe=lme(FMform,random=~1|Nest,data=Owls)
summary(FMe)
anova(FMgls,FMe)

#STEP 6 - VALIDATE RANDOM PART OF FULL MODEL
#STANDARD VALIDATION PLOT
plot(FMe)
#HOMEMADE PLOTS
op=par(mfrow = c(3,2)) 
plot(Owls$SexParent, resid(FMe),pch=19,col='green',xlab="SexParent",ylab="Residuals")
plot(Owls$FoodTreatment, resid(FMe),col='red',xlab="FoodTreatment",ylab="Residuals",pch=20)
plot(Owls$ArrivalTime, resid(FMe),col='blue',xlab="ArrivalTime",ylab="Residuals",pch=20)
hist(resid(FMe),col="brown",main='',xlab="Residuals")
boxplot(resid(FMe)~Nest,data=Owls,col='green',ylab='Residuals')
boxplot(residuals(FMe,type='pearson')~Nest,data=Owls,col='red',ylab='Standardized Residuals')
par(op)

#STEPS 7-8 - FIND OPTIMAL REDUCED MODEL FIXED STRUCTURE using ML
FMe.ml=update(FMe,method="ML") 
summary(FMe.ml)
anova(FMe.ml,type="marginal")
RMe1.ml=update(FMe.ml,.~.-SexParent:FoodTreatment,data=Owls)
summary(RMe1.ml)
anova(RMe1.ml,type="marginal")
anova(FMe.ml,RMe1.ml)
RMe2.ml=update(RMe1.ml,.~.-SexParent:ArrivalTime,data=Owls)
summary(RMe2.ml)
anova(RMe2.ml,type="marginal")
anova(RMe1.ml,RMe2.ml)
RMe3.ml=update(RMe2.ml,.~.-SexParent,data=Owls)
summary(RMe3.ml)
anova(RMe3.ml,type="marginal")
anova(RMe2.ml,RMe3.ml)

#STEP 9 REFIT OPTIMAL MODEL with REML AND VALIDATE
RMe3=update(RMe3.ml,method="REML")
summary(RMe3)
anova(RMe3)
#STANDARD VALIDATION PLOT
plot(RMe3)
#HOMEMADE PLOTS
op=par(mfrow = c(3,2)) 
plot(Owls$SexParent, resid(RMe3),pch=19,col='green',xlab="SexParent",ylab="Residuals")
plot(Owls$FoodTreatment, resid(RMe3),col='red',xlab="FoodTreatment",ylab="Residuals",pch=20)
plot(Owls$ArrivalTime, resid(RMe3),col='blue',xlab="ArrivalTime",ylab="Residuals",pch=20)
hist(resid(RMe3),col="brown",main='',xlab="Residuals")
boxplot(resid(RMe3)~Nest,data=Owls,col='green',ylab='Residuals')
boxplot(residuals(RMe3,type='pearson')~Nest,data=Owls,col='red',ylab='Standardized Residuals')
par(op)



#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#ADDITIONAL STEPS 3-5 - CHOOSE A VARIANCE STRUCTURE using REML
FMe=lme(FMform,random=~1|Nest,data=Owls)
summary(FMe)
anova(FMe)

#ADDING CORRELATION STRUCTURES FOR Nest USING weights OPTION IN gls()
FMgls1=gls(FMform,data=Owls,weights=varIdent(form=~1|Nest))
summary(FMgls1)
anova(FMgls1)
anova(FMe,FMgls1)  #FMe strongly preferred over FMgls1
anova(FMgls,FMgls1)  #FMgls1 with varIdent variance structure slightly preferred by test & AIC but not BIC

#ADDING varIdent CORRELATION STRUCTURE FOR Nest USING weights OPTION IN gls() & lme()
FMe2=lme(FMform,random=~1|Nest,data=Owls,weights=varIdent(form=~1|Nest),control=lmeControl(tolerance=1e-3))
#method fails to converge here

#ADDING COMPOUND SYMMETRIC (NESTED) CORRELATION STRUCTURE USING weights OPTION IN gls()
FMgls2=gls(FMform,data=Owls,correlation=corCompSymm(form=~1|Nest))
summary(FMgls2)  #IDENTICAL RESULTS TO FMe
anova(FMgls2)

#ADDING corAR1 AUTOCORRELATION FOR FoodTreatment WITHIN Nest TO glm(FM)
FMgls3=gls(FMform,data=Owls,correlation=corAR1(form=~1|FoodTreatment))
summary(FMgls3)
anova(FMgls3)
anova(FMe,FMgls3)  #FMgls3 STRONGLY PREFERRED

#ADDING corAR1 AUTOCORRELATION FOR Nest/FoodTreatment TO lme(FM)
FMe1=lme(FMform,random=~1|Nest,data=Owls,correlation=corAR1(form=~1|Nest/FoodTreatment))
summary(FMe1)
anova(FMe1)
anova(FMe1,FMe) #FMe1 STRONGLY SUPPORTED 
anova(FMe1,FMgls3) #NESTED FMe1 WITH AUTOCORRELATION STRONGLY PREFERRED OVER gls() RESULT

#ADDING corAR1 AUTOCORRELATION FOR Nest/FoodTreatment TO lme(FM)
FMe2=lme(FMform,random=~1|Nest/FoodTreatment,data=Owls,correlation=corAR1(form=~1|Nest/FoodTreatment))
summary(FMe2)
anova(FMe2)
anova(FMe2,FMe1)  #MODELS ARE NEARLY EQUIVALENT, FMe1 with lower AIC & BIC

#ADDING Nest/FoodTreatment 2-level lme(FM) ONLY
FMe3=lme(FMform,random=~1|Nest/FoodTreatment,data=Owls)
summary(FMe3)
anova(FMe3)
anova(FMe3,FMe) #COMPARING ONE VS TWO LEVEL NESTING 
anova(FMe3,FMe1)  
anova(FMe3,FMe2)

AIC(FM,FMgls3,FMe1,FMe2,FMe3)

#STEP 6 - VALIDATE RANDOM PART OF FULL MODEL
#STANDARD VALIDATION PLOT
plot(FMe1)
#HOMEMADE PLOTS
op=par(mfrow = c(3,2)) 
plot(Owls$SexParent, resid(FMe1),pch=19,col='green',xlab="SexParent",ylab="Residuals")
plot(Owls$FoodTreatment, resid(FMe1),col='red',xlab="FoodTreatment",ylab="Residuals",pch=20)
plot(Owls$ArrivalTime, resid(FMe1),col='blue',xlab="ArrivalTime",ylab="Residuals",pch=20)
hist(resid(FMe1),col="brown",main='',xlab="Residuals")
boxplot(resid(FMe1)~Nest,data=Owls,col='green',ylab='Residuals')
boxplot(residuals(FMe1,type='pearson')~Nest,data=Owls,col='red',ylab='Standardized Residuals')
par(op)

#STEPS 7-8 - FIND OPTIMAL REDUCED MODEL FIXED STRUCTURE using ML
FMe1.ml=update(FMe1,method="ML") 
summary(FMe1.ml)
anova(FMe1.ml,type="marginal")
RMe1.ml=update(FMe1.ml,.~.-SexParent:FoodTreatment,data=Owls)
summary(RMe1.ml)
anova(RMe1.ml,type="marginal")
anova(FMe1.ml,RMe1.ml)
RMe2.ml=update(RMe1.ml,.~.-SexParent:ArrivalTime,data=Owls)
summary(RMe2.ml)
anova(RMe2.ml,type="marginal")
anova(RMe1.ml,RMe2.ml)
RMe3.ml=update(RMe2.ml,.~.-SexParent,data=Owls)
summary(RMe3.ml)
anova(RMe3.ml,type="marginal")
anova(RMe2.ml,RMe3.ml)

AIC(FMe1.ml,RMe1.ml,RMe2.ml,RMe3.ml)

#STEP 9 REFIT OPTIMAL MODEL with REML AND VALIDATE
RMe3=update(RMe3.ml,method="REML")
summary(RMe3)
anova(RMe3)
#STANDARD VALIDATION PLOT
plot(RMe3)
#HOMEMADE PLOTS
op=par(mfrow = c(3,2)) 
plot(Owls$SexParent, resid(RMe3),pch=19,col='green',xlab="SexParent",ylab="Residuals")
plot(Owls$FoodTreatment, resid(RMe3),col='red',xlab="FoodTreatment",ylab="Residuals",pch=20)
plot(Owls$ArrivalTime, resid(RMe3),col='blue',xlab="ArrivalTime",ylab="Residuals",pch=20)
hist(resid(RMe3),col="brown",main='',xlab="Residuals")
boxplot(resid(RMe3)~Nest,data=Owls,col='green',ylab='Residuals')
boxplot(residuals(RMe3,type='pearson')~Nest,data=Owls,col='red',ylab='Standardized Residuals')
par(op)