Hello,
Thank you for another thorough response.
I've decided to stick with the ordinal model, although I can't work out
where to put the pl=TRUE to work out the range of latent variables.
I just used Vanimal/(Vanimal+Vresidual) to get heritability, rather than
Vanimal/(Vanimal+Vresidual+1). I have now tried the latter though and
the answer is very similar.
Unfortunately running the model for longer and giving it a longer burnin
has not reduced the autocorrelation. When you say 'parameter expansion',
does that mean of G1 in the prior?
I've carried on using this prior:
prior2<-list(G=list(G1=list(V=matrix(p.var/2),n=1)),
R=list(V=matrix(p.var/2),n=1))
in which I think the variance is set to 1.
Thank you again: all help is very much appreciated.
Helen
On 12/07/2012 14:13, Jarrod Hadfield wrote:
> Hi,
>
> Quoting Helen Ward on Wed, 11 Jul 2012 11:58:50
> +0100:
>
>> Dear Jarrod and David,
>>
>> Thank you both very much for your responses: I'm sorry it's taken me
>> so long to say that, but I wanted to have a proper go at everything
>> before I replied, which took longer than I anticipated!
>>
>> Thank you especially for your zipoisson code, Jarrod. When I
>> initially attempted to use your prior and model MCMCglmm asked me for
>> a stronger prior, so I changed it to
>>
>>
>> Prior1=list(R=list(V=diag(2),nu=1, fix=2), G=list(G1=G1)) : is that a
>> reasonable thing to have done? The model
>>
>>
> I think it is better to diagnose the problem rather than upping nu. In
> this case I think the residual variance of the Poisson process must be
> going to zero under the flat prior? Perhaps the counts are
> under-dispersed with respect to the Poisson, further motivating an
> ordinal model?
>
>
>
> model1.1<-MCMCglmm(Toes~trait,random=~us(trait):animal,family="zipoisson",rcov=~idh(trait):units,pedigree=ped,data=data,prior=Prior1,nitt=500000,thin=500,burnin=200000,verbose=FALSE)
>
>>
>> ran with this prior, however I admit I am now flummoxed as to how to
>> proceed and go about estimating heritability.
>>
>>
>> For this reason, (and because you both recommended it!), I went on to
>> try an ordinal model using the following code:
>>
>> p.var<-var(data$Toes,na.rm=TRUE)
>>
>> prior2<-list(G=list(G1=list(V=matrix(p.var/2),n=1)),
>> R=list(V=matrix(p.var/2),n=1))
>>
>>
> IMPORTANT: The residual variance is not identifiable in an ordinal
> model so all the information is coming from the prior. This is OK, but
> I would fix the variance at 1 so as to avoid interpretation
> difficulties (and numerical problems - check the range of the latent
> variables using pl=TRUE).
>
>
> model2<-MCMCglmm(Toes~1,random=~animal,family="ordinal",pedigree=ped,data=data,prior=prior2,nitt=500000,thin=500,burnin=200000,verbose=FALSE)
>
>>
>>
>> followed by the same model but with the prior from the Wilson paper
>> erratum:
>>
>> prior2.1 <- list(G = list(G1 = list(V = 1, nu = 0.002)), R = list(V =
>> 1,nu = 0.002))
>>
>>
>> the plots for these models look ok, and they both give a heritability
>> estimate of ~0.55, with very similar 95% HDP intervals, although the
>> auto correlation is quite high ~0.6
>
> Is this calculated as Vanimal/(Vanimal+Vresidual+1) or
> Vanimal/(Vanimal+Vresidual)?
>
>>
>>
>> Next I tried a categorical model, treating Toes as a binary response
>> variable (deforemed or not deformed), using
>>
>> p.var<-var(data$ToesBinary,na.rm=TRUE)
>>
>> prior3= list(R = list(V = 1, n = 0, fix = 1), G = list(G1 = list(V =
>> 1,n = 0)))
>>
>> model3<-
>> MCMCglmm(ToesBinary~1,random=~animal,family="categorical",data=data,pedigree=ped,prior=prior3,nitt=200000,thin=250,burnin=100000,verbose
>> = FALSE)
>>
>>
>> Heritability came out higher with this model at about ~0.85.
>
> Is this calculated as Vanimal/(Vanimal+Vresidual+pi^2/3) or
> Vanimal/(Vanimal+Vresidual)?
>>
>>
>> And then last but not least a multinomial2 model:
>>
>> prior4=list(R=list(V=1,nu=0.002),G=list(G1=list(V=1,nu=1,alpha.mu=0,alpha.V=100)))
>>
>>
>> model4<-MCMCglmm(cbind(Deformed,Fine)~1,random=~animal,family="multinomial2",data=data,pedigree=ped,prior=prior4,nitt=200000,thin=250,burnin=100000,verbose=FALSE)
>>
>>
>> The variance plots for this one didn't look very pretty and the
>> autocorrelation values were still quite high, ~0.8, but it gave me a
>> pretty heritability plot and an estimate similar to my ordinal models.
>>
>>
>> As you can probably gather from all of this, I am something of an
>> indecisive person...., but at this stage I am inlcined to stick with
>> the ordinal models and believe that there is a heritable element to
>> the number of deformed toes (which seems rather intuitive after all
>> that!). Given this decision, can you recommend anything I might do to
>> reduce the autocorrelation in these models. Also, do I have to run
>> the same model many times to confirm its conclusions?
>
> Re: autocorrelation. Use parameter expansion and/or run for longer
> (particularly for binary animal models and/or if there is little
> overdispersion) If the absolute value of the latent variables under
> the logit link exceed 20, or under the probit link exceed 7, then at
> the moment I would consider another program unless the latent
> variables are associated with particular levels of a fixed effect.
> Then flatish priors on the probability scale may overcome the problem.
>
>
>>
>>
>> Finally finally, I will have to do another zipoisson model in the
>> future, can you direct me to where I might find an example of the
>> code used for estimating heritability in these cases please?
>
> I think it makes most sense to think about the heritabilities of each
> process separately - the first as a Poisson trait and the second as
> a threshold (binary) trait.
>
>>
>>
>> Many many thanks once more,
>> Helen
>>
>>
>>
>>
>> On 08/07/2012 12:03, Jarrod Hadfield wrote:
>>> Hi,
>>>
>>> I agree with David that an ordinal model may behave better for this
>>> type of data (or even a 0/1 response with deformed or not) but in
>>> case you want to pursue a ZIP model...
>>>
>>> The model you have specified is probably not what you had in mind.
>>> You should think of the ZIP response as being two responses
>>> (traits): the Poisson counts (trait 1) and the probability that a
>>> zero comes from the zero-inflation process (trait 2). At the moment
>>> you have a common intercept for both and therefore assume that the
>>> probability that a zero is from the zero-inflation process (on the
>>> logit scale) is equal to the mean Poisson count (on the log scale).
>>> There is probably little justification for this (but see ZAP
>>> models). Moreover, you have a single animal term, and therefore
>>> assume that both processes have the same genetic variance and that
>>> the genetic correlation between them is 1. You also assume that the
>>> residual variance for each process is equal (but the residual
>>> correlation is 0), although this is not actually that bigger deal
>>> because residual variation in the zero-inflation process (and
>>> residual corrletaion) is not identifiable from the data and so the
>>> value it could take is essentially arbitrary. However, I like to
>>> recognise this non-identifiability by fixing the variance at 1 and
>>> the covariance at zero (but see ZAP models where your specification
>>> is useful).
>>>
>>> So, a more appropriate model would be:
>>>
>>> model1.1<-MCMCglmm(Toes~trait,random=~us(trait):animal,family="zipoisson",rcov=~idh(trait):units,pedigree=ped,data=data,prior=prior1.1,nitt=500000,thin=500,burnin=200000,verbose=FALSE)
>>> or perhaps
>>>
>>> model1.2<-MCMCglmm(Toes~trait,random=~idh(trait):animal,family="zipoisson",rcov=~idh(trait):units,pedigree=ped,data=data,prior=prior1.1,nitt=500000,thin=500,burnin=200000,verbose=FALSE)
>>> The difference between these two is that the genetic correlation
>>> between the zero-inflation and the Poisson processes is estimated,
>>> and in the latter it is set to 0.
>>>
>>> I'm not sure about how much data you have, but I would think that if
>>> the incidence of deformed toes is very low, you probably have little
>>> information for some of the parameters. In this case you have to be
>>> VERY careful with priors. The prior recommendations in the MCMCglmm
>>> tutorial in Wilson et al. 2011 will generally be quite informative
>>> for small or even moderate sized datasets (the authors have kindly
>>> published an erratum as many people seem to be using their
>>> recommendation without realising its implications).
>>>
>>> The prior could look like:
>>>
>>> prior=list(R=list(V=diag(2), nu=0, fix=2), G=list(G1=G1))
>>>
>>> which fixes the residual variance in the zero-inflation process to 1
>>> with a flat improper prior on the Poisson over-dispersion. For the
>>> prior on the genetic (co)variances (G1) I generally use parameter
>>> expanded priors of the form:
>>>
>>> G1=list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000)
>>>
>>> or
>>>
>>> G1=list(V=diag(2), nu=1, alpha.mu=c(0,0), alpha.V=diag(2)*1000)
>>>
>>> for model.2.
>>>
>>> These generally give posterior modes for the variance components
>>> that are close to (RE)ML estimates, although it is unclear to me how
>>> they behave for the covariances and functions of variances (e.g.
>>> heritability). At least for binary traits (e.g. the zero-inflation
>>> bit) I have seen a chi-square prior recommended because of its
>>> prior properties on the heritability (de Villemereuil in MEE). This
>>> type of prior can also be obtained using parameter expansion:
>>>
>>> (V=1, nu=1000, alpha.mu=0, alpha.V=1)
>>>
>>> and I would guess that the multivariate analogue would be
>>> (V=diag(2), nu=1001, alpha.mu=c(0,0), alpha.V=diag(2)). However, a
>>> chi-square prior probably has poor properties for the Poisson
>>> variance and in MCMCglmm you could only mix these different priors
>>> for model.2 using a different syntax:
>>>
>>> model1.2<-MCMCglmm(Toes~trait,random=~idh(at(trait,1)):animal+idh(at(trait,2)):animal,family="zipoisson",rcov=~idh(trait):units,pedigree=ped,data=data,prior=prior1.1,nitt=500000,thin=500,burnin=200000,verbose=FALSE)
>>> G=list(
>>> G1=list(V=1, nu=1, alpha.mu=0, alpha.V=1000),
>>> G2=list(V=1, nu=1000, alpha.mu=0, alpha.V=1)
>>> )
>>>
>>>
>>> Cheers,
>>>
>>> Jarrod
>>>
>>>
>>>
>>>
>>> Quoting Helen Ward on Wed, 04 Jul 2012
>>> 16:24:36 +0100:
>>>
>>>> Dear list,
>>>>
>>>> I work on bats, some of which have deformed toes. I have just started
>>>> trying to use MCMCglmm to investigate whether there is hertiable
>>>> variation for the number of deformed toes a bat has. Since most bats
>>>> have 0 deformed toes I am using a zero-inflated Poisson model.
>>>>
>>>> I have put together some code (see below) based on the MCMCglmm
>>>> tutorial
>>>> from the Wilson et al. 2011 'Ecologist's guide to the animal model'
>>>> paper, the MCMCglmm Course Notes and trial and error. This code
>>>> runs and
>>>> I get a model and a heritability estimate. However, the heritability
>>>> estimate I get is extremely high and this, teamed with the
>>>> observations
>>>> that my posterior distributions of the variance componants are not
>>>> pretty as the example ones I've seen and my autocorrelation values are
>>>> high (~0.4), makes me suspicious that I am not modelling this trait
>>>> very
>>>> well. I suspect I have over simplified the model. I based my prior
>>>> on a
>>>> model done in ASReml modelling 'Toes' as a Poisson distribution, which
>>>> can an heritability estimate of 0.35.
>>>>
>>>> In addition, when I try and add a second random variable - year of
>>>> birth
>>>> (Born) - into a second model, both plots and auto correlation get
>>>> worse
>>>> and my heritability estimate for number of deformed toes (Toes) goes
>>>> from about 0.8 to 0.05!
>>>>
>>>> Finally, I suspect you're not, but I would like to know if you are
>>>> allowed to compare DIC values from models with different priors
>>>> please.
>>>>
>>>> I have pasted my code below. If anyone can offer constructive
>>>> criticism
>>>> it would be much appreciated.
>>>>
>>>> Many thanks,
>>>> Helen
>>>>
>>>> Model 1: With number of deformed toes (Toes) as the only random
>>>> variable
>>>>
>>>>
>>>> p.var<-var(data$Toes,na.rm=TRUE)
>>>> prior1.1<-list(G=list(G1=list(V=matrix(p.var*0.35),n=1)),
>>>> R=list(V=matrix(p.var*0.65),n=1))
>>>>
>>>> model1.1<-MCMCglmm(Toes~1,random=~animal,family="zipoisson",rcov=~trait:units,pedigree=ped,data=data,prior=prior1.1,nitt=500000,thin=500,burnin=200000,verbose=FALSE)
>>>> posterior.heritability1.1<-model1.1$VCV[,"animal"]/(model1.1$VCV[,"animal"]+model1.1$VCV[,"trait:units"])
>>>> HPDinterval(posterior.heritability1.1,0.95)
>>>> posterior.mode(posterior.heritability1.1)
>>>>
>>>>
>>>> Model 2: With number of toes (Toes) and year of birth (Born) as random
>>>> variables
>>>>
>>>>
>>>> p.var<-var(data$Toes,na.rm=TRUE)
>>>> prior1.2<-list(G=list(G1=list(V=matrix(p.var/3),n=1),G2=list(V=matrix(p.var/3),n=1)),R=list(V=matrix(p.var/3),n=1))
>>>> model1.2<-
>>>> MCMCglmm(Toes~1,random=~animal+Born,family="zipoisson",rcov=~trait:units,pedigree=ped,data=data,prior=prior1.2,nitt=1000000,thin=1000,burnin=500000,verbose=FALSE)
>>>> posterior.heritability1.2<-model1.2$VCV[,"animal"]/(model1.2$VCV[,"animal"]+model1.2$VCV[,"Born"]+model1.2$VCV[,"trait:units"])
>>>> HPDinterval(posterior.heritability1.2,0.95)
>>>> posterior.mode(posterior.heritability1.2)
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models@r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>
>>>>
>>>
>>>
>>>
>>
>>
>>
>
>
>
[[alternative HTML version deleted]]