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