[R-sig-ME] Multi-response MCMCglmm (Gaussian and zapoisson)

Jarrod Hadfield j.hadfield at ed.ac.uk
Mon Feb 21 22:23:23 CET 2011

Hi Susanne,

A multi-response model with two Gaussian and one zero-altered Poisson  
(ZAP) responses cannot be currently fitted in MCMCglmm due to  
restrictions on the R-structure. The residual variance of the ZA part  
of the ZAP cannot be estimated (it is a binary variable) and so needs  
to be set to some value, or as you indicate set to the same residual  
variance as the P part. The covariance between the ZA and the P part  
is also non-identifiable from the data (no datum can be both zero and  
non-zero) and it is customary to set it to zero. A fully parameterised  
residual (co)variance matrix for the model should ideally look like:

   V(1)  C(1,2) C(1,3) C(1,4)
  C(1,2)  V(2)  C(2,3) C(2,4)
  C(1,3) C(2,3)  V(3)    0
  C(1,4) C(2,4)   0    V(3)

Currently this is not possible, and the "us" variance function that  
you use estimates V(4) and C(3,4) for which all the information is  
coming from the prior. You can probably go one step further and  
restrict V(4)=1 using parameter expanded priors, but you are still  
left with estimating C(3,4). There may well be transformations that  
give interpretable marginal distributions that are independent of  
C(2,3) and V(4) but I do not know them.

The (co)variance matrices for the random effects can be fully estimated.

More generally, these models are VERY complex and you would need a lot  
of data, with lots of replication at the appropriate levels.  
Personally, I would try and come up with a simple, biologically  
plausible model first and see how you go from there. The priors are  
likely to be informative and in the case of the residual (co)variance  
matrix are improper. Remember that the marginal prior on an individual  
variance in a covariance matrix has degree of belief parameter  
nu-dim(V)+1 which for your random effect covariance matrices is 1  
(informative) and -2.996 for your residual covariance matrix (improper).

m2 looks more reasonable, but again has many fixed and random effects  
so I would worry whether the data-set is large enough and well  
structured enough to estimate all parameters with reasonable  
precision. You may want to consider parameter expanded priors in this  
case. Note also, that in the first set of CourseNotes with a  
zero-altered Poisson section (v2.05-v2.08, online from Aug-Nov 2010) I  
made a mistake and said that if a model of the form

trait:(expl1 + expl2+ expl3+ expl4+ expl5….)+trait

is fitted then the significant "za" terms indicate zero-deflation/inflation.

This was wrong and is corrected in subsequent updates (v2.12 is  
current). The "za" terms in teh parameterisation:

trait*(expl1 + expl2+ expl3+ expl4+ expl5….)

indicate zero-deflation/inflation when significant. The model is the  
same but the interpretation is different. I'm sorry for this.



Quoting Susanne Lachmuth <susanne.lachmuth at botanik.uni-halle.de>:

> 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 (reproduction) with the size parameters (i.e. responses  
> of m1) as covariates. We are interested in a trade-off of growth and  
> reproduction.
> Any help would be greatly appreciated.
> Many thanks,
> Susanne Lachmuth
> -------------------------------
> Susanne Lachmuth
> Department of Plant Ecolgy
> University of Halle-Wittenberg
> Germany
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

More information about the R-sig-mixed-models mailing list