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