[R-sig-ME] MCMCglmm query
James Higham
jhigham at dpz.eu
Tue Aug 30 15:23:39 CEST 2011
Dear Jarrod,
Thanks very much for your fast and detailed reply.
Apologies for the delay in response; I wanted to make sure I could get
these changes working before responding.
The answer to the sign question makes perfect sense - thanks a lot.
Thanks also for the tips on offsetting in MCMCglmm and the other
suggestions and corrections. I've played around with it and including
the logged offset as a covariate and setting the co-efficient to 1 with
very low variance seems to work well for both the regular poisson and
zero-inflated poisson models. Likewise, the suggested syntax for
applying the offset only to the non-zi part of the model seems to be
working well.
I believe that the variables are zero-inflated based on comparison of
predicted zero numbers under the poisson and those observed in the
dataset, but I'll take another look using the MCMCglmm checks.
Thanks again so much for your help - it's very kind of you to take the
time to ensure that so many people are able to use your package
successfully for their projects.
All the best
James
On Wed, Aug 24, 2011 at 7:40 PM, Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote:
> Hi,
>
> The coefficients are opposite in sign because the zero-inflation
> predictors are predicting excess zeros not excess non-zero's. Given
> the small value of traitzi_ConsortCount
> (-3.8491) are you sure the data are zero-inflated and not just
> over-dispersed? See the CourseNotes 5.3-5.5.
>
> Also, offset(Obstime) will not result in an offset in MCMCglmm, you
> need to fit Obstime as a covariate and then fix the associated
> regression coefficient to one in the prior (see
> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q4/005004.html).
> You probably want to log it (Obstime), and you may only want to fit it
> for the non zi term using something like
> at.level(trait,1):log(Obstime).
>
> The residual variance for the zi term cannot be estimated so it is
> more usual to set it tome value (e.g. 1) in the prior:
>
> R=list(V=diag(2), nu=a, fix=2)
>
> where a is some number. Generally a=2 is not so non-informative - I
> use parameter expanded priors if I want relative non-informativeness.
>
> I would also worry about temporal autocorrelation within subjects,
> something which MCMCglmm is not too flexible in handling.
>
> Jarrod
>
>
>
>
>
> On 24 Aug 2011, at 17:34, James Higham wrote:
>
>
> > Dear MCMCglmm users,
> >
> > I have been playing with MCMCglmm for days and have got one or two
> > things that I am confused about. I am new both to MCMCglmm, but also
> > to
> > MCMC methods more generally. I do not have access to a statistician
> > for
> > further advice and would be extremely grateful to anyone who might
> > be
> > kind enough to spare some time to offer any.
> >
> > Briefly, my data are in the form of behavioural data counts for
> > different days during primate female reproductive cycles. There are
> > multiple days of observation for each cycle (so multiple
> > observations of
> > the same female within a cycle) with a total of 31 different female
> > cycles. These females come from multiple social groups. I want to
> > undertake univariate mixed models to investigate how the count of a
> > given behaviour (e.g. MatingCount) offset for the number of
> > observation
> > scans undertaken on that day (offset response variable) may be
> > related
> > to factors such as cyclephase (fertile, pre-fertile, post-fertile,
> > measured hormonally, fixed) and cycletype (conceptive vs
> > nonconceptive,
> > fixed), while controlling for having multiple days from the same
> > female
> > cycles (FemaleID, random) and multiple females from the same social
> > groups (Group, random). Some response variables are zero-inflated
> > (e.g.
> > females don’t mate at all on most days of the cycle) and as it’s not
> > clear whether the sources of all zeros are the same for all
> > individuals
> > (e.g. low ranking females may want to mate on more days of the cycle
> > than they do but are excluded from doing so by high ranking females)
> > zero-inflated poisson models seem appropriate.
> >
> >
> > I have the following model structure following previous posts
> >
> >
> > (e.g.
> >
> > https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q2/002283.html)
> >
> >
> >
> >
> > zip1<-MCMCglmm(MatingCount~offset(Obstime)+trait+trait:cyclephase+trait:cycletype-1,random=~idh(trait):FemaleID+idh(trait):Group,prior=prior,thin=10,nitt=100000,rcov=~idh(trait):units,pl=TRUE,family="zipoisson",data=data)
> >
> > I have followed previous threads (e.g.
> >
> > https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/005546.html)
> > in
> > setting the priors as
> >
> >
> > prior <- list(R=list(V=diag(2),n=2), G=list(G1=list(V=diag(2), n=2),
> > G2=list(V=diag(2), n=2)))
> >
> > This is producing models that appear to have good mixing (from
> > plots,
> > Geldman diagnostic showing all point est. 1.01-1.02).
> >
> >
> > Here are example results for “ConsortCount”
> > post.mean
> > l-95% CI u-95% CI eff.samp pMCMC
> > traitConsortCount 5.4271
> > 3.5951 7.2812 9700.00 0.00144 **
> > traitzi_ConsortCount -3.8491
> > -6.1093 -1.4728 125.49 0.00722 **
> > traitConsortCount:cyclephase -0.3890 -0.5794
> > -0.2273 6091.00 < 1e-04 ***
> > traitzi_ConsortCount:cyclephase 2.0585 1.4343
> > 2.7447 30.13 < 1e-04 ***
> > traitConsortCount:cycletype -0.0166 -0.3515
> > 0.3127 8189.13 0.91938
> > traitzi_ConsortCount:cycletype 0.2823 -0.5088
> > 1.1460 8189.13 0.91938
> >
> > My problem with this output is that I do not understand why the
> > poisson
> > and zi results for cyclephase have the opposite sign. The variable
> > is
> > coded such that 0 = fertile phase, 1 = pre-fertile phase (we'd
> > predict
> > counts to be lower than in the fertile phase), 2 = post-fertile
> > phase
> > (prediction is that counts are lower still). I would expect consorts
> > both to be less likely as the cyclephase value increases, but also
> > that
> > it would occur less often when it does occur as cyclephase value
> > increases (i.e. both zi and poisson parts would have negative
> > terms).
> > Following this, I have split each of the variables into two,
> > removing
> > all the zeros from one dataset and modelling this part as a poisson
> > mixed model, and in the other part converting all values >0 into 1
> > and
> > modelling this part as a binomial mixed model. These models were
> > undertaken using the lmer function in lme4 and give negative terms
> > (as
> > predicted) for both binomial and poisson parts of each variable. I
> > have
> > also confirmed that this appears to be the case by plotting mean
> > values
> > of e.g. ConsortCount by cyclephase in both the new binomial and
> > poisson
> > datasets separately – both show that mean ConsortCount is highest in
> > phase 0, then phase 1, then phase 2. I could present the lmer
> > results,
> > but I would prefer to use MCMCglmm if possible(models seem to be
> > struggling to converge in lmer if I include the offset variable, the
> > poisson part of the data would be overdispersed which I believe
> > MCMCglmm
> > accounts for, and it seems more elegant to analyse the data using
> > one
> > model in any case).
> >
> >
> >
> > If anyone can see an obvious explanation for this that I’m missing
> > I’d
> > be very grateful. Also, if anyone has got this far and has noticed
> > anything awry in my model specifications I’d be extremely grateful
> > if
> > they could let me know!
> >
> > Apologies if anything above doesn't make sense and/or I've
> > misunderstood
> > something.
> >
> >
> > Best wishes
> >
> > James Higham
> > Post-doc
> > German Primate Center
> >
> >
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
> --
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
>
More information about the R-sig-mixed-models
mailing list