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