## GLMM Examples                   
require(lme4)

## fix glmer to allow for # obs >= # random effects
## execute this:
fixInNamespace(glmer_finalize, 'lme4')
## and change the following line:
## if (!(length(levels(dm$flist[[1]])) < ncol(dm$Zt)))
## to
##  if (!(length(levels(dm$flist[[1]])) <= ncol(dm$Zt)))
## save it and close the editor window/frame to return to R.

## this hack does only work for some versions of lme4!!

#####################################################################
## example from glm talk continued
data(fabric, package='npmlreg')
## counts of the number of faults in rolls of fabric of different lengths
## x = log(leng)
str(fabric)

summary(mod1 <- glm(y ~ x, fabric, family = poisson))     
summary(mod1 <- glm(y ~ x, fabric, family = quasipoisson)) 

(mp.fit <- glmer(y~x + (1|x), fabric, family=poisson)) ## this works alright


####################################################################
## Owls Data Set, see
## http://glmm.wikidot.com/examples
## second example

## get data
download.file("http://www.highstat.com/Book2/ZuurDataMixedModelling.zip",
              "tmp.zip")
unzip("tmp.zip",files="Owls.txt")
Owls <- read.table("Owls.txt",header=TRUE)

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

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

head(Owls)

## fit first model candidate
g1 <- glmer(SiblingNegotiation~FoodTreatment*SexParent+offset(log(BroodSize))+
            (1|Nest),family=poisson,data=Owls)
print(summary(g1))

## test overdispersion
rdev <- sum(residuals(g1)^2)
mdf <- length(fixef(g1))
rdf <- nrow(Owls)-mdf  ##  residual df [NOT accounting for random effects]
rdev/rdf
#Lots of overdispersion because it is much bigger than 1

#Significance test of overdispersion
(prob.disp <- pchisq(rdev,rdf,lower.tail=FALSE,log.p=TRUE))

#P-value is 10E-377.

## second model candidate: use quasipoisson
g2 <- glmer(SiblingNegotiation~FoodTreatment*SexParent+offset(log(BroodSize))+
            (1|Nest),family=quasipoisson,data=Owls)
print(summary(g2))

## third model candidate: add per observation random effect
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))

## plot residuals: notice zero inlated data
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")

## calculate and plot predictions
library(ggplot2)

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()

## 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)))
pvar1 <- diag(mm %*% tcrossprod(vcov(g1),mm))
pvar2 <- diag(mm %*% tcrossprod(vcov(g2),mm))
tvar1 <- pvar1+VarCorr(g1)$Nest
tvar2 <- pvar2+VarCorr(g2)$Nest
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)))
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))

## functions to help comparing the various model fits
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
}

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

## and glmmPQL
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")

## aggregate and display results
mtab1 <- modelTab(g1,g2,g3,horiz=TRUE)
library("nlme")
mtab2 <- modelTab(g4,g5,g6,horiz=TRUE)
cbind(mtab1,mtab2)
detach("package:nlme")