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

[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:

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


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