[R-sig-ME] priors for a multi-response model (MCMCglmm)

Ned Dochtermann ned.dochtermann at gmail.com
Wed Oct 13 20:04:33 CEST 2010

Thanks a lot for the comments.
Regarding 1 and 2 below, sorry about that--those were actually typos from
trying to simplify the code and make it more generic. Both aspects of the
code were actually specified as you suggest; sorry for the sloppiness.

3 and 4 really look like the key issues for this analysis (besides the
number of parameters being estimated which has been a concern throughout).
Unfortunately those points suggest that the best alternative is to estimate
the covariance matrix using a Poisson distribution, despite the known
zero-inflation. Under the family statement in the help for MCMCglmm a
ztpoisson distribution is mentioned however no zero-truncated distribution
is mentioned in the course notes. Is this something that was previously
available but has been removed?


Ned Dochtermann
Department of Biology
University of Nevada, Reno

ned.dochtermann at gmail.com

-----Original Message-----
From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk] 
Sent: Tuesday, October 12, 2010 2:43 AM
To: Ned Dochtermann
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] priors for a multi-response model (MCMCglmm)

Hi Ned,

I think you haven't specified the model that you want to fit.

1) you have "random=~us(trait):units, rcov=~us(trait):units" which is  
fitting the same terms for the random effects and the residuals. I  
think you want something like  "random=~us(trait):individual,  
rcov=~us(trait):units"? At the moment any separation of these effects  
is coming entirely from the prior.

2) depending on what your responses are "~trait+period+day" for the  
fixed terms is probably not appropriate. You probably want to interact  
period and day with trait, so that the effect of these two predictors  
can have independent effects on each response

3) zero-inflated poisson are treated as multi-response. The first  
latent variable is for the poisson counts, and the second is for the  
zero-inflation. The residual variance for the zero-inflated part  
cannot be estimated from the data (for the same reasons that the  
residual variance in a binary model is not identified) and neither can  
the residual covariance between zero-inflation and the poisson counts  
(you only see one of the processes in any one observation).  I  
generally place strong priors on them by fixing the residual variance  
to one and the residual covariance to zero. I achieve this by fitting  
an idh residual structure idh(trait):units which fixes the residual  
covariance to zero, and fix the 2nd diagonal element of the residual  
matrix (the zero-inflation variance) to one in the prior:  
"prior=list(R=diag(2), nu=...,  fix=2)"

4)  If you have 3 ZIP responses you have 6 "traits" to worry about.   
It is not possible to place the appropriate constraints on this  
matrix: the sub-diagonals cannot be estimated (corresponding to the  
covariance between the poisson counts and the zero-inflation within  
traits) and neither can the even elements of the diagonal (the 3 zero- 
inflation terms). You may be able to specify a prior which has a  
desirable marginal effect on the things you want to calculate, but I  
think this would be hard and a lot of work.

5) These are very parameter rich models, and I would avoid them unless  
you have a lot of individuals measured many times.



On 11 Oct 2010, at 21:06, Ned Dochtermann wrote:

> Hi all,
> While I'm still struggling with properly specifying priors in  
> general, I've
> run into a specific problem I can't quite muddle through. I'm trying  
> to
> estimate the covariances among several behaviors with repeated  
> measures per
> individual. I initially did so using the following structure:
> multi.prior<- 
> list(G=list(G1=(list(V=diag(3),nu=3))),R=list(V=diag(3),nu=3))
> #I know this assumes unit variance
> multi.trait.p1<-MCMCglmm(cbind(trait1,trait2,trait3)~trait,
> random=~us(trait):units, rcov=~us(trait):units,
> family=c("poisson","poisson","poisson), data=Compiled,  
> prior=multi.prior,
> verbose=FALSE)
> or if including fixed factors:
> multi.trait.p2<-MCMCglmm(cbind(trait1,trait2,trait3)~trait+period+day,
> random=~us(trait):units, rcov=~us(trait):units,
> family=c("poisson","poisson","poisson), data=Compiled,  
> prior=multi.prior,
> verbose=FALSE)
> (I actually run both a lot longer than the defaults but I've left  
> that out
> here as it isn't relevant)
> Both of these models seem to work well, they give reasonable answers  
> and
> satisfy a variety of diagnostics.
> However, in looking back over the data I realized the data had  
> pretty severe
> zero-inflation. Thus, I've tried to rerun the analyses using zero- 
> inflated
> models. Based on the MCMCglmm course notes I thought that the first  
> step for
> the priors would be to expand both G and R to diag(4):
> multi.prior.zip<- 
> list(G=list(G1=(list(V=diag(4),nu=4))),R=list(V=diag(4),nu=
> 4)) #the last 'nu' is wrong based on how ZIP model priors are  
> specified in
> the course notes
> multi.trait.zip<-MCMCglmm(cbind(trait1,trait2,trait3)~trait,
> random=~us(trait):units, rcov=~us(trait):units,
> family=c("zipoisson","zipoisson","zipoisson),
> data=Compiled,prior=multi.prior.zip,verbose=FALSE)
> This doesn't work as "V is the wrong dimension for some priorG/priorR
> elements". I also suspect it is more generally wrong due to the  
> random and
> rcov statements and issues with estimating aspects of the zero- 
> inflation and
> poisson covariances; however I'm specifically interested in  
> estimating the
> covariance matrix so I don't want to use an idh specification here.
> I'd like to get the covariance matrix from a ZIP model but I'm not  
> sure what
> all the errors in the above coding are nor the solutions. Basically  
> I know
> both the specification of the prior and the specification of the  
> model are
> wrong. Any help would be greatly appreciated.
> Thanks,
> Ned
> --
> Ned Dochtermann
> Department of Biology
> University of Nevada, Reno
> ned.dochtermann at gmail.com
> http://wolfweb.unr.edu/homepage/mpeacock/Dochter/
> --
> _______________________________________________
> 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