library(coda)
library(MCMCglmm)
library(reshape)
load("../DATA/Owls.rda")



Owls <- rename(Owls,c(SiblingNegotiation="NCalls"))
## setup to facilitate model variations below:
offvec <- c(1,1,2,rep(1,5)) ## indicates whether fixed effect is the offset
                            ## (log brood size) or not

## Specifying the fixed effect is a little bit tricky.
## For zero-inflated models, \code{MCMCglmm} internally constructs an 
## augmented data frame: if the original data frame is 

## \begin{verbatim}
## resp      f1  x1
## 1         A   1.1
## 0         A   2.3
## 3         B   1.7
## \end{verbatim}

## where \code{response} is a Poisson response, \code{f1} is a factor,
## and \code{x1} is a continuous predictor, then the augmented data frame
## will look like this:
## \begin{verbatim}
## resp      trait   f1  x1
## 1         resp    A   1.1
## 0         resp    A   2.3
## 3         resp    B   1.7
## 1         zi_resp A   1.1
## 0         zi_resp A   2.3
## 1         zi_resp B   1.7
## \end{verbatim}

## \code{MCMCglmm} provides a special helper function, \code{at.level}, which
## enables specification of covariates that affect only the count part of
## the model, or only the binary part of the model: the first would be
## specified as an interaction of \code{at.level(trait,1)} with a covariate
## (i.e., the covariate affects only responses for level~1 of trait, which
## are the count responses).  

## Here is the model specified by Zuur et al.,
## which includes an offset effect of brood size and most but not all of
## the interactions between \code{FoodTreatment}, \code{SexParent},
## and \code{ArrivalTime} on the count aspect of the model, but only a single
## fixed (intercept) effect on the binary (zero-inflation). 

fixef2 <- NCalls~trait-1+
  at.level(trait,1):logBroodSize+
  at.level(trait,1):((FoodTreatment+ArrivalTime)*SexParent)

prior_overdisp  <- list(R=list(V=diag(c(1,1)),nu=0.002,fix=2),
                        G=list(list(V=diag(c(1,1e-6)),nu=0.002,fix=2)))

prior_overdisp_broodoff <- within(prior_overdisp,
                                  { B <- list(mu=c(0,1)[offvec],
                                              V=diag(c(1e8,1e-6)[offvec]))})

st1 <- system.time(ofit_MCMCglmm <- MCMCglmm(fixef2,
                                     rcov=~idh(trait):units,
                                     random=~idh(trait):Nest,
                                     prior=prior_overdisp_broodoff,
                                     data=Owls,
                                     family="zipoisson",
                                     verbose=FALSE))

if (FALSE) {
  library(glmmADMB)
  st1 <- system.time(gfit1 <-
                     glmm.admb(fixed=NCalls~(FoodTreatment+ArrivalTime)*SexParent,
                               offset="logBroodSize",
                               random=~1,
                               group="Nest",
                               data=Owls,
                               zeroInflation=TRUE,
                               family="poisson"))
}
                  
## need to write results section
ofit_MCMCglmm_mcmc <- as.mcmc(ofit_MCMCglmm$Sol)
results <- list(obj=NULL,
                coef=colMeans(ofit_MCMCglmm_mcmc)[offvec==1],
                confint.mcmc=(HPDinterval(ofit_MCMCglmm_mcmc))[offvec==1,],
                confint.quantile=t(apply(ofit_MCMCglmm_mcmc,2,
                  quantile,c(0.025,0.975)))[offvec==1,],
                time=c(fit=st1,mcmc=st1),
                terminfo=NULL,
                convinfo=list(effsize=effectiveSize(ofit_MCMCglmm_mcmc),
                  geweke=geweke.diag(ofit_MCMCglmm_mcmc)))

save("results",file="fit.RData")