###################################################
### chunk number 1: 
###################################################
options(continue=" ")


###################################################
### chunk number 2:  eval=FALSE
###################################################
## download.file("http://www.highstat.com/Book2/ZuurDataMixedModelling.zip",
##               file="tmp.zip")
## unzip("tmp.zip",files="Owls.txt")


###################################################
### chunk number 3: 
###################################################
Owls <- read.table("Owls.txt",header=TRUE)


###################################################
### chunk number 4: 
###################################################
library(lattice)
print(bwplot(reorder(Nest,NegPerChick)~NegPerChick|FoodTreatment:SexParent,
       data=Owls))


###################################################
### chunk number 5: 
###################################################
print(dotplot(reorder(Nest,NegPerChick)~NegPerChick|FoodTreatment:SexParent,
       data=Owls))


###################################################
### chunk number 6: 
###################################################
library(lme4)


###################################################
### chunk number 7: fit1
###################################################
g1 <- glmer(SiblingNegotiation~FoodTreatment*SexParent+offset(log(BroodSize))+
            (1|Nest),family=poisson,data=Owls)
print(summary(g1))


###################################################
### chunk number 8: 
###################################################
rdev <- sum(residuals(g1)^2)
mdf <- length(fixef(g1))
rdf <- nrow(Owls)-mdf  ##  residual df [NOT accounting for random effects]
rdev/rdf


###################################################
### chunk number 9: 
###################################################
(prob.disp <- pchisq(rdev,rdf,lower.tail=FALSE,log.p=TRUE))


###################################################
### chunk number 10: fit2
###################################################
Owls$obs <- 1:nrow(Owls) ## add observation number to data
g2 <- glmer(SiblingNegotiation~FoodTreatment*SexParent+offset(log(BroodSize))+
            (1|Nest)+(1|obs),family=poisson,data=Owls)
print(summary(g2))


###################################################
### chunk number 11: 
###################################################
plot(fitted(g2),residuals(g2))
rvec <- seq(0,30,length=101)
lines(rvec,predict(loess(residuals(g2)~fitted(g2)),newdata=rvec),
      col=2,lwd=2)
abline(h=0,col="gray")


###################################################
### chunk number 12: 
###################################################
library(ggplot2)


###################################################
### chunk number 13: 
###################################################
G0 <- ggplot(Owls,aes(x=reorder(Nest,NegPerChick),
                      y=NegPerChick))+
  xlab("Nest")+ylab("Negotiations per chick")+coord_flip()+
  facet_grid(FoodTreatment~SexParent)
## boxplot display
G1 <- G0+ geom_boxplot()
## dotplot display (I prefer this one)
G2 <- G0+stat_sum(aes(size=factor(..n..)),alpha=0.5)+
  theme_bw()


###################################################
### chunk number 14: 
###################################################
## set up prediction frame
pframe0 <- with(Owls,expand.grid(SexParent=levels(SexParent),
                                FoodTreatment=levels(FoodTreatment)))
## construct model matrix
mm <- model.matrix(~FoodTreatment*SexParent,data=pframe0)
## predictions from each model; first construct linear
##  predictor, then transform to raw scale
pframe1 <- data.frame(pframe0,eta=mm%*%fixef(g1))
pframe1 <- with(pframe1,data.frame(pframe1,NegPerChick=exp(eta)))
pframe2 <- data.frame(pframe0,eta=mm%*%fixef(g2))
pframe2 <- with(pframe2,data.frame(pframe2,NegPerChick=exp(eta)))


###################################################
### chunk number 15: 
###################################################
pvar1 <- diag(mm %*% tcrossprod(vcov(g1),mm))
pvar2 <- diag(mm %*% tcrossprod(vcov(g2),mm))


###################################################
### chunk number 16: 
###################################################
tvar1 <- pvar1+VarCorr(g1)$Nest
tvar2 <- pvar2+VarCorr(g2)$Nest


###################################################
### chunk number 17: 
###################################################
pframe1 <- data.frame(pframe1,pse=sqrt(pvar1),tse=sqrt(tvar1))
pframe1 <- with(pframe1,
                data.frame(pframe1,
                           plo=exp(eta-1.96*pse),
                           phi=exp(eta+1.96*pse),
                           tlo=exp(eta-1.96*tse),
                           thi=exp(eta+1.96*tse)))
pframe2 <- data.frame(pframe2,pse=sqrt(pvar2),tse=sqrt(tvar2))
pframe2 <- with(pframe2,
                data.frame(pframe2,
                           plo=exp(eta-1.96*pse),
                           phi=exp(eta+1.96*pse),
                           tlo=exp(eta-1.96*tse),
                           thi=exp(eta+1.96*tse)))


###################################################
### chunk number 18: 
###################################################
print(G2 + 
      geom_hline(data=pframe1,aes(yintercept=NegPerChick),col="red")+
      geom_hline(data=pframe2,aes(yintercept=NegPerChick),col="blue")+
      geom_rect(aes(xmin=0,xmax=28,ymin=tlo,ymax=thi,x=NULL),
                data=pframe1,fill="red",alpha=0.3)+
      geom_rect(aes(xmin=0,xmax=28,ymin=tlo,ymax=thi,x=NULL),
                data=pframe2,fill="blue",alpha=0.3))


###################################################
### chunk number 19: 
###################################################
Fixef <- function(x) {
  if (inherits(x,"mer") || inherits(x,"lme")) { fixef(x) } else
  if (inherits(x,"glm")) { coef(x)  } else
  stop("unknown model type")
}
         
modelTab <- function(...,horiz=FALSE,round=3,cFun=Fixef) {
  L1 <- list(...)
  n <- length(L1)
  coefs <- lapply(L1,cFun)
  allcoefs <- Reduce(union,lapply(coefs,names))
  cmat <- matrix(NA,nrow=n,ncol=length(allcoefs),
                 dimnames=list(NULL,allcoefs))
  for (i in 1:n) {
    cmat[i,names(coefs[[i]])] <- coefs[[i]]
  }
  if (horiz) cmat <- t(cmat)
  cmat <- round(cmat,round)
  cmat
}


###################################################
### chunk number 20: 
###################################################
library(geepack)
g3 <- geeglm(SiblingNegotiation~FoodTreatment*SexParent+offset(log(BroodSize)),
       corstr="exchangeable",id=Nest,
       family=poisson,data=Owls)


###################################################
### chunk number 21: 
###################################################
g4 <- MASS::glmmPQL(SiblingNegotiation~FoodTreatment*SexParent+
                    offset(log(BroodSize)),
                    random=~1|Nest,
                    family=poisson,data=Owls)
g5 <- update(g4,family=quasipoisson)
g6 <- update(g4,random=~1|Nest/obs)
detach("package:nlme")


###################################################
### chunk number 22: 
###################################################
mtab1 <- modelTab(g1,g2,g3,horiz=TRUE)
library("nlme")
mtab2 <- modelTab(g4,g5,g6,horiz=TRUE)
cbind(mtab1,mtab2)
detach("package:nlme")


###################################################
### chunk number 23:  eval=FALSE
###################################################
## library(MCMCglmm)
## MCMCglmm()