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

[Jarrod Hadfield]

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.

Cheers,

Jarrod



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.