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