mixed zero-inflated poisson regression

Hi Dr. Jarrod.

Didn't you want to say:

"prior<-list(R=list(V=diag(2), nu=0.002, fix=2)..." ?

Thank you.
Denise.

2011/11/14 Jarrod Hadfield <j.hadfield at ed.ac.uk>
Hi,

ZIP models generally mix poorly and often need to be run for long  
periods of time. However, you seem to be estimating a residual  
variance for the zero-inflation bit and this is not identifiable. In  
this case you will never get convergence or good mixing. I recommend  
you fix the residual variance of the zero-inflation process at one  
using a prior like:

prior<-list(R=list(V=diag(2), nu=0.002, fix=1), G=list(G1=list(V=1,  
nu=1, alpha.mu=0, alpha.V=1000)))

I have used a parameter expansion for the random effects. However,  
your specification is a bit odd because you assume that the random  
effects for the poisson and zero-inflation have equal variance and a  
correlation of 1. ~us(trait):gr  or  ~idh(trait):gr are probably  
more reasonable and you *could* use the priors:

prior<-list(R=list(V=diag(2), nu=0.002, fix=1),  
G=list(G1=list(V=diag(2), nu=2, alpha.mu=c(0,0),  
alpha.V=diag(2)*1000)))

for ~us(trait):gr, and

prior<-list(R=list(V=diag(2), nu=0.002, fix=1),  
G=list(G1=list(V=diag(2), nu=1, alpha.mu=c(0,0),  
alpha.V=diag(2)*1000)))

~idh(trait):gr.

~gr might make more sense if  you use a zero-altered poisson  
("zapoisson") rather than a zero-inflated poisson.

Cheers,

Jarrod




On 14 Nov 2011, at 16:17, Ivar Herfindal wrote:

Dear Thierry

Thank you for your quick response. I have previously tried to follow  
the CourseNote you suggested without success. I gave it another try  
today, and it seems I have manged to generate such a zero-inflated  
model as I want during the day. The syntax looks something like this:

mcmodm1 <- MCMCglmm(Y ~ trait -1 + at.level(trait, 1:2):A, data=dataa,
  family="zipoisson", rcov=~idh(trait):units, nitt=30000, random=~gr)

However, the main challenge is that the model seems quite unstable  
withmy "real" data. Even with 100000 iterations, two identical  
specified models gives quite different output, and plotting the  
mcmodm1$Sol suggest not very stable estimates. This occurs also on  
my simulated data.

Do anyone have any suggestion on how to make the model more stable  
(e.g. by the rcov or random arguments)? Or do I just have to run it  
as long as it takes to get "fairly precise" estimates?

Ivar

On 14.11.2011 10:40, ONKELINX, Thierry wrote:
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