Skip to content

priors for a multi-response model (MCMCglmm)

4 messages · Ned Dochtermann, Jarrod Hadfield

Original
Ned Dochtermann
# 
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/
--
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:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      Ned Dochtermann
    

    

      
      #
    

  


  

      
Jarrod,
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?



Thanks, 
Ned

--
Ned Dochtermann
Department of Biology
University of Nevada, Reno

ned.dochtermann at gmail.com
website
--



-----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.

Cheers,

Jarrod

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

        

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jarrod Hadfield
    

    

      
      #
    

  


  

      
Hi Ned,

The zero-truncated poisson ("ztpoisson") is implemented and as far as  
I know works fine. In terms of model syntax it will be the same as a  
normal poisson model. The zero's would have to be discarded though.

Cheers,

Jarrod

        
On 13 Oct 2010, at 19:04, Ned Dochtermann wrote: