R: mixed zero-inflated poisson regression

 Dear Ivar,

 Your assumption on MCMCglmm is incorrect. You can specify different predictor for both the count and the binomial component. Have a look at the CourseNotes vignette for an example.

 Best regards,

 Thierry

 ----------------------------------------------------------------------------
 ir. Thierry Onkelinx
 Instituut voor natuur- en bosonderzoek
 team Biometrie&   Kwaliteitszorg
 Gaverstraat 4
 9500 Geraardsbergen
 Belgium

 Research Institute for Nature and Forest
 team Biometrics&   Quality Assurance
 Gaverstraat 4
 9500 Geraardsbergen
 Belgium

 tel. + 32 54/436 185
 Thierry.Onkelinx at inbo.be
 www.inbo.be

 To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
 ~ Sir Ronald Aylmer Fisher

 The plural of anecdote is not data.
 ~ Roger Brinner

 The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
 ~ John Tukey

 -----Oorspronkelijk bericht-----
 Van: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-
 bounces at r-project.org] Namens Ivar Herfindal
 Verzonden: maandag 14 november 2011 9:57
 Aan: r-sig-mixed-models at r-project.org
 Onderwerp: [R-sig-ME] mixed zero-inflated poisson regression

 Dear all

 It seems like I have reached the limit of my statistical skills, and now turn to the
 mixed-list for some assistance.

 I have a dataset with a response variable that appears to have a distribution that
 is zero-inflated poisson (ZIP) (or close to). The response variable is number of
 fledglings from nests with known mothers.
 The number of fledglings has a lot of zeros that comes from failed breeding
 attempts. The number of fledgings can therefore be considered as a two-step
 process:
 1. do the female breed, and
 2. if the female breed, what was the number of fledglings.

 I want to relate the number of fledglings to some predictor variables, and I am
 particularly interested in the combined effect of these variables on the
 probability to breed and number of fledglings given breeding. I have several
 observations per female (multiple years) and each nest location is also
 represented by several females (in different years). I therefore need/want to
 account for random factors. I guess I could run two separate models (one
 binomial and one poisson), but I need to make predictions on how my
 explanatory variables affect the number of fledglings irrespective of the
 mechanism behind.

 I have searched R-sites a lot, and it found some packages that can model ZIP.
 1. pscl (function zeroinfl)
 Pros: have the option to relate both the poisson and binomial process to the
 predictor variable, so I can evaluate how probability of breeding and number of
 fledglings depend on my predictor variables
 Cons: Can not include random factors.

 2. glmmADMB (function glmmadmb), MCMCglmm, and gamlss
 Pros: allows random factors
 Cons: 1. all packages seems to assume that all zeros have the same probability
 of belonging to the zero component, i.e. the logit-link structure does not depend
 on the predictor variable.
 2. For some of the packages (particularly MCMCglmm) I have some challenges
 about how to specify the model correctly.

 Does anyopne in the R(mixed) community have experience on this? I provide a
 simple example below to illustrate what I want.

 Ivar Herfindal

 ## A simple reproducible example
 library(pscl)
 library(glmmADMB)
 library(MCMCglmm)
 library(gamlss)
 set.seed(4242)

 # create some dummy data that have a lot of zeros in the dependent variable
 dataa<- cbind.data.frame(Y=rbinom(100, 1, 0.7)) dataa$Y<- dataa$Y *
 rpois(100, 3) # the response dataa$A<- rnorm(100) # explanatory variable
 dataa$gr<- factor(1:10) # grouping variable

 # fit models with different packages

 ## First non-mixed models
 ### pscl zeroinfl
 zimod<- zeroinfl(Y ~ A, data=dataa, dist="poisson") # Get estimate of A on
 poisson and binomial process separate, this is what I want (I believe so...)

 zimod2<- zeroinfl(Y ~ A|1, data=dataa, dist="poisson") # All zero counts have
 the same probability of belonging to the zero component

 ### glmmADMB
 admod<- glmmadmb(Y ~ A, data=dataa, family="poisson", zeroInflation=TRUE)
 # Same results as zimod2.
 # Q: Is it possible to obtain similar results as zimod, i.e.
 # separate the effect of A on the poisson and binomial process?

 ### MCMCglmm
 mcmod<- MCMCglmm(Y ~ A, data=dataa, family="zipoisson",
 rcov=~idh(trait):units)
 # more or less same as admod and zimod2, but "plot(mcmod$Sol)" reveals large
 variation.
 # Can change burn and nitt to get more stable(?)

 mcmod2<- MCMCglmm(Y ~ A, data=dataa, family="zipoisson",
 rcov=~idh(trait):units, burn=10000, nitt=20000) # Q: Is it possible to separate the
 effect of A on the poisson and binomial process?

 ### gamlss
 gamod<- gamlss(Y ~ A, data=dataa, family=ZIP) # same as zimod2 and admod

 ## Some attempts on mixed models with gr as random factor admodm<-
 glmmadmb(Y ~ A, data=dataa, family="poisson", zeroInflation=TRUE,
 random=~(1|gr)) mcmodm<- MCMCglmm(Y ~ A, data=dataa,
 family="zipoisson", rcov=~idh(trait):units, random=~idh(trait):gr) gamodm<-
 gamlss(Y ~ A + random(gr), data=dataa, family=ZIP) # For glmmadmb and gamlss
 I think the models are specified ok # For MCMCglmm I am not sure if it is
 specified correctly # But I still want to be able to separate the effect of A on the
 two processes.

 --
 --------------------------------------------------
 Ivar Herfindal (PhD, Researcher)
 Centre for Conservation Biology
 Department of Biology
 Norwegian University for Science and Technology
 N-7491 Trondheim, Norway

 email:     Ivar.Herfindal at bio.ntnu.no
 web:       http://www.ntnu.no/employees/ivar.herfindal
 phone:     +47 73596253
 Cellphone: +47 91625333
 fax:       +47 73596090