[R-sig-ME] MCMCglmm query

[James Higham]

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