# from Chapters 4 and 12-13 in Zuur et al, mixed models + some additions
rm(list=ls())
Owls=read.delim("Owls.txt")

# some simple plotting
boxplot(NegPerChick~Nest,data=Owls)
boxplot(NegPerChick~FoodTreatment,data=Owls)
boxplot(NegPerChick~SexParent,data=Owls)
plot(x=Owls$ArrivalTime,y=Owls$NegPerChick)

# a simple linear model to look at residuals
M.lm=lm(NegPerChick~SexParent*FoodTreatment+SexParent*ArrivalTime,data=Owls)
plot(M.lm,select=c(1))                                                           # only residuals~fitted

Owls$LogNeg<-log10(Owls$NegPerChick+1)                                           # transforming the response; note that adding 1 is arbitrary and should be assessed
M2.lm=lm(LogNeg~SexParent*FoodTreatment+SexParent*ArrivalTime,data=Owls)
plot(M2.lm,select=c(1))
E=rstandard(M2.lm)

op<-par(mar=c(3,3,2,2))
boxplot(E~Nest,data=Owls,axes=FALSE,ylim=c(-3,3))                                # distribution of residuals per nest, no large heterogeneity
abline(0,0)
axis(2)
text(1:27,-2.5, levels(Owls$Nest),cex=0.75,srt=65)
par(op)

# Step 2 of protocol in Zuur et al Mixed Models
library(nlme)                                                                    # one library used to fit mixed models
Form<-formula(LogNeg~SexParent*FoodTreatment+SexParent*ArrivalTime)              # create a forumla object
M.gls=gls(Form,data=Owls)                                                        # without the random effect nest
summary(M.gls)

M1.lme=lme(Form,random=~1|Nest,method="REML",data=Owls)                          # with the random effect nest, only for intercept, fitted using restricted/residual maximum lik
anova(M.gls,M1.lme)                                                              # comparing with/without random effect; such tests are not always correct, see eg Pinheiro & Bates

E2<-resid(M1.lme,type="normalized")                                              # plotting residuals~fitted values for the model
F2<-fitted(M1.lme)
op<-par(mfrow=c(2,2),mar=c(4,4,3,2))
MyYlab="Residuals"

plot(x=F2,y=E2,xlab="Fitted values",ylab=MyYlab)
boxplot(E2~SexParent,data=Owls,main="Sex of parent",ylab=MyYlab)                 # checking variances of residuals are appr the same for different groups
boxplot(E2~FoodTreatment,data=Owls,main="Food treatment",ylab=MyYlab)
plot(x=Owls$ArrivalTime,y=E2,main="Arrival time",ylab=MyYlab,xlab="Time (hours)")
par(op)

# step 7 of the protocol
summary(M1.lme)                                                                  # summary, interpret the parameters, intercept meaningless as such
M1b.lme=lme(LogNeg~SexParent*FoodTreatment+SexParent*I(ArrivalTime-24),random=~1|Nest,method="REML",data=Owls)   # using time-24 to get a meaningful intercept, ie value for midnight
summary(M1b.lme)                                                                 # new summary
anova(M1.lme)                                                                    # sequential P-values (removing terms one after each other); differs from summary() except from the last term

M1.Full=lme(Form,random=~1|Nest,method="ML",data=Owls)
M1.A=update(M1.Full,.~.-SexParent:FoodTreatment)
M1.B=update(M1.Full,.~.-SexParent:ArrivalTime)
anova(M1.Full,M1.A)                                                              # comparing models, AIC, BIC or P-values (LRT) can be used
anova(M1.Full,M1.B)

Form2<-formula(LogNeg~SexParent+FoodTreatment+SexParent*ArrivalTime)
M2.Full=lme(Form2, random= ~1| Nest, method = "ML", data = Owls)
M2.A=update(M2.Full, .~. -FoodTreatment)
M2.B=update(M2.Full, .~. -SexParent:ArrivalTime)
anova(M2.Full,M2.A)
anova(M2.Full,M2.B)

Form3 <- formula(LogNeg~SexParent+FoodTreatment+ArrivalTime)
M3.Full <- lme(Form3, random= ~1| Nest, method = "ML", data = Owls)
M3.A <- update(M3.Full, .~. -FoodTreatment)
M3.B <- update(M3.Full, .~. -SexParent)
M3.C <- update(M3.Full, .~. -ArrivalTime)
anova(M3.Full,M3.A)
anova(M3.Full,M3.B)
anova(M3.Full,M3.C)

Form4 <- formula(LogNeg ~ FoodTreatment + ArrivalTime)
M4.Full <- lme(Form4, random= ~1| Nest, method = "ML", data = Owls)
M4.A <- update(M4.Full, .~. -FoodTreatment)
M4.B <- update(M4.Full, .~. -ArrivalTime)
anova(M4.Full,M4.A)
anova(M4.Full,M4.B)

M5 <- lme(LogNeg ~ FoodTreatment + ArrivalTime, random= ~1| Nest, method = "REML", data = Owls)
summary(M5)

E2<-resid(M5,type="normalized")                                                                     # plotting residuals of "best" model
F2<-fitted(M5)
op<-par(mfrow=c(2,2),mar=c(4,4,3,2))
MyYlab="Residuals"

plot(x=F2,y=E2,xlab="Fitted values",ylab=MyYlab)
boxplot(E2~SexParent,data=Owls,main="Sex of parent",ylab=MyYlab)
boxplot(E2~FoodTreatment,data=Owls,main="Food treatment",ylab=MyYlab)
plot(x=Owls$ArrivalTime,y=E,main="Arrival time",ylab=MyYlab,xlab="Time (hours)")
par(op)

library(lattice)

xyplot(E2~ArrivalTime|SexParent*FoodTreatment,data=Owls,
  ylab="Residuals",xlab="Arrival time (hours)",
      panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})


library(mgcv)                                                                          # using mixed additive models to assess linearity of relationships
M6 <- gamm(LogNeg ~ FoodTreatment + +s(ArrivalTime), random=list(Nest=~1),data=Owls)
        
summary(M6)
plot(M6$gam)
anova(M6$gam)
summary(M6$gam)


M7=gamm(NegPerChick~FoodTreatment+
       s(ArrivalTime,by=as.numeric(FoodTreatment=="Deprived"))+
       s(ArrivalTime,by=as.numeric(FoodTreatment=="Satiated")),
        random=list(Nest=~1),data=Owls)

M8=gamm(NegPerChick~FoodTreatment+
       s(ArrivalTime,by=as.numeric(SexParent=="Female"))+
       s(ArrivalTime,by=as.numeric(SexParent=="Male")),
        random=list(Nest=~1),data=Owls)

AIC(M6$lme,M7$lme,M8$lme)

# using model selection based on all models and AICc
# treats random effects as one parameter, prediction at population level (ie average over nests, not for a given nest)
require(MuMIn)
M1.lme=lme(LogNeg~SexParent*FoodTreatment+SexParent*I(ArrivalTime-24),random=~1|Nest,method="ML",data=Owls)
dredge(M1.lme)

############################################
# chapter 12 of Zuur et al
# Using the number of calls as the response; use an offset (ie not estimating the corresponding parameter) to get a ratio
# a Poisson or quasipoisson (Poisson with an extra dispersion parameter, fitted using quasilikelihood) used by default a link function = log
# hence log(Ncalls) ~ log(Broodsize) + other variables is the same as log(Ncalls/Broodsize) ~ other variables 
# First models fitted without random effects (ie generalized linear models)
Owls$NCalls<-Owls$SiblingNegotiation
Owls$LBroodSize<-log(Owls$BroodSize)

Form<-formula(NCalls~offset(LBroodSize)+ SexParent*FoodTreatment+ SexParent*ArrivalTime)
O1<-glm(Form,family=quasipoisson,data=Owls)
drop1(O1,test="F")
O2<-update(O1,.~. -SexParent*ArrivalTime)
drop1(O2,test="F")

Form<-formula(NCalls~offset(LBroodSize)+ FoodTreatment+ ArrivalTime)
O3<-glm(Form,family=quasipoisson,data=Owls)
drop1(O3,test="F")

# chapter 13 of Zuur, using random effects, ie generalized linear mixed models
library(lme4)                                                         # nlme library could have been used together with MASS (function glmmPQL) but lme4 is used here
Owls$NCalls<-Owls$SiblingNegotiation
Owls$LBroodSize<-log(Owls$BroodSize)
Owls$fNest<-factor(Owls$Nest)

# doing a stepwise model selection
O1.lmer<-lmer(NCalls~offset(LBroodSize)+
 SexParent*FoodTreatment+SexParent*ArrivalTime+(1|fNest),data=Owls,family=poisson)
summary(O1.lmer)
display(O1.lmer)

O2.lmer<-lmer(NCalls~offset(LBroodSize)+ SexParent*FoodTreatment+
              SexParent+ArrivalTime+(1|fNest),data=Owls, family=poisson)
anova(O1.lmer,O2.lmer,test="F")

O3.lmer<-lmer(NCalls~offset(LBroodSize)+ FoodTreatment+SexParent+
              ArrivalTime+(1|fNest),data=Owls,family=poisson)
summary(O3.lmer)
anova(O2.lmer,O3.lmer,test="F")

O4.lmer<-lmer(NCalls~offset(LBroodSize)+
              FoodTreatment+
              ArrivalTime+(ArrivalTime|fNest),data=Owls,
              family=poisson)
summary(O4.lmer)
anova(O3.lmer,O4.lmer,test="F")

# the library MuMIn can be used here too
dredge(O1.lmer)

# checking residuals
E4=residuals(O4.lmer,type="normalized")
F4=fitted(O4.lmer)
plot(F4,E4)
xyplot(E4~F4|FoodTreatment*SexParent,data=Owls)

# checking linearity
library(mgcv)
O4.gamm<-gamm(NCalls~offset(LBroodSize)+
              FoodTreatment+s(ArrivalTime),
              random=list(fNest=~1),data=Owls,
              family=poisson)

summary(O4.gamm$gam)
anova(O4.gamm$gam)
plot(O4.gamm$gam)

summary(O4.gamm$lme)

E4<-resid(O4.gamm$lme,type="normalized")

O4.gamm<-gamm(NCalls~offset(LBroodSize)+
              FoodTreatment+s(ArrivalTime),
              random=list(fNest=~1),data=Owls,
              family=quasi(link=log,variance="mu"))
              
              
summary(O4.gamm$gam)

intervals(O4.gamm$lme,which="var-cov")

eta <- fitted(O4.gamm$lme)+Owls$LBroodSize
fv  <- exp(eta)
res.raw <- Owls$NCalls - fv
res.P <- (Owls$NCalls - fv) / sqrt(fv )
sd(res.P)