Susanne Lachmuth
2011-Feb-17 16:08 UTC
[R] Multi-response MCMCglmm (gaussian and zapoisson)
Dear MCMCglmm users,
I am currently struggling with the specification of a proper prior and model
formula for a multi-response MCMCglmm with two of the three response variables
being Gaussian and the third being za-poisson. The model includes several fixed
effects and three nested random effects.
In general, I would prefer to fit a model with a fixed effect of trait and
suppressed intercept for getting trait specific intercepts. Further, I wish to
measure covariances between the response variables and consequently to specify
completely parameterized covariance matrices.
My current model looks like this:
prior<-list(R=list(V=diag(4),nu=0.004),G=list(G1=list(V=diag(4)/4,n=4),G2=list(V=diag(4)/4,n=4),G3=list(V=diag(4)/4,n=4)))
m1<-MCMCglmm(fixed=cbind(resp1,resp2,resp3) ~
trait:(expl1 + expl2+ expl3+ expl4+ expl5...)+trait-1,
random=~us(trait):rand1+us(trait): rand1: rand2+us(trait): rand1: rand2: rand3,
rcov=~us(trait):units,family=c("gaussian","gaussian","zapoisson"),nitt=20000,burnin=5000,thin=10,prior=prior,data=dat)
However, for zero-altered models it is recommended in the MCMCglmm course notes
to constrain over-dispersion to be the same for both processes (the
zero-alteration and poisson process) by using a trait by unit interaction in the
R-structure. Additionally, the intercept should not be suppressed for getting
the differences between the zero-altered regression coefficients and the Poisson
regression coefficients. This allows identification of zero-inflation or
zero-deflation in response to the explanatory variables.
I therefore fit an additional model for the zapoisson response only, looking
like this:
prior2<-list(R=list(V=diag(1)/4,nu=0.002),G=list(G1=list(V=diag(1)/4,nu=0.002),G2=list(V=diag(1)/4,nu=0.002),G3=list(V=diag(1)/4,nu=0.002)))
m2<-MCMCglmm(fixed=resp3 ~
trait:(expl1 + expl2+ expl3+ expl4+ expl5?.)+trait,
random=~trait:rand1+trait: rand1: rand2+ trait: rand1: rand2: rand3,
rcov=~trait:units,
family="zapoisson",nitt=20000,burnin=1000,thin=10,data=dat,
prior=prior2)
The results show that there is ?significant? zero-inflation and deflation in
response to some variables.
My main questions are:
Are the two priors specified correctly?
Does it make sense to include the zapoisson response (resp3) in model m1 and is
the model formula (in particular the R-structure) appropriate?
An alternative might be to analyze resp1 and resp2 (both size variables of
plants) with model m1 and fit an extra model (like m2) for resp 3 (Number of
flowers) with the size parameters (i.e. responses of m1) as covariates. We are
interested in a trade-off between growth and reproduction.
Any help would be greatly appreciated.
Many thanks,
Susanne Lachmuth
--------------------
Susanne Lachmuth
Department of Plant Ecolgy
Martin Luther University Halle-Wittenberg
Germany
