[R-sig-ME] MCMCglmm query

[Jarrod Hadfield]

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.