[R-sig-ME] MCMCglmm query
Jarrod Hadfield
j.hadfield at ed.ac.uk
Wed Aug 24 19:40:44 CEST 2011
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