[R] Multi-response MCMCglmm (gaussian and zapoisson)

[Susanne Lachmuth]

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