[R-sig-ME] MCMCglmm query

Jarrod Hadfield j.hadfield at ed.ac.uk
Wed Aug 24 19:40:44 CEST 2011


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, 

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.


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