Hi
I am estimating the heritability of an ordinal trait using MCMCglmm and 
have come across two problems: one regarding heritability and one 
specific to ordinal data sets.

_My data_
Trait1 = ordinal score from 0 - 4
Colony = factor with 2 options
ID = repeated measures per individual (818 individuals)

/prior2.1 <- list(G = list(G1=list(V=1, nu=1000, alpha.mu=0, alpha.V=1),/
/G2=list(V=1, nu=1000, alpha.mu=0, alpha.V=1)), /
/R = list(V = 1, fix = 1)) /

The simplest model just estimates Va (animal) and Vi (ID):
/model2.1<- MCMCglmm(Trait1~ Colony , random =~animal+ID, pedigree = 
Ped3, data = Data, prior = prior2.1, family='ordinal',burnin = 20000, 
nitt = 500000, thin = 200, pr=TRUE)/

When I extract the variance components:
/posterior.mode(model2.1$VCV)/
/HPDinterval(model2.1$VCV)/

/animal = //0.89 (0.05-1.69/
/ID = 0.61 (0.00-1.60)/

So the heritability:

/posterior.heritability1.1 <- model2.1$VCV[, "animal"]/(model2.1$VCV[, 
"animal"]+ model2.1$VCV[, "ID"] + model2.1$VCV[, "units"]+1)/
/posterior.mode(posterior.heritability1.1)/
/h^2 = 0.21 (0.05- 0.45)/

/Now if I fit the full model (prior structure changed accordingly): /

/model2.2<-MCMCglmm(Trait1~ Colony , random =~animal + BYEAR + MOTHER + 
ID, pedigree = Ped3, data = Data, prior = prior2.1, 
family='ordinal',burnin = 20000, nitt = 500000, thin = 200, pr=TRUE)/

And then extract the variance components and heritability as above:

/animal = //0.00 (0.00-0.37)///

/BYEAR = //0.00 (0.00-1.47)///

/MOTHER=//0.00 (0.00-0.25)///

/ID=//0.00 (0.00-1.45)/

/h^2 = /0.00 (0.00-0.37)

I have read about using DIC for model comparisons for non guassian data 
and I think I understand why it is problematic 
(https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q2/018096.html).

However, for Gaussian models, using MCMCglmm, to estimate heritability, 
model selection using DIC is the norm (e.g. Charmantier et al 2011, J 
Evol Bio; Taylor et al 2012, J Evol Bio).I have seen some posts that 
suggest that model selection is unnecessary for Bayesian models but as 
you can see from the output above, including all possible random effects 
in this model pushes the variance components to zero, and dramatically 
changes the interpretation of the results.I was just wondering if anyone 
has any thoughts on how to chose the best model or does it simple 
involve deciding the best biological reasoning, and fitting this model?I 
am relatively new to these analyses but the dramatic change in variance 
estimates seems worrying to me.

My second question is specific to ordinal analyses:As I understand, for 
any fixed effect, I can estimate a probability matrix, such that I get a 
probability of having a specific score:

https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003671.html

So if I fit the model:
/model2.1<- MCMCglmm(Trait1~ ID , random =~ year, pedigree = Ped3, data 
= Data, prior = prior2.1, family='ordinal',burnin = 20000, nitt = 
500000, thin = 200, pr=TRUE)/

And there are 800 individuals (ID) and 5 ordinal levels of trait 1, I 
think I will have a 800 x 5 matrix, with a probability of each 
individual having each value of trait 1.  As in example the first line 
may read:

Individual01234

10.90.050.0500

However, I need to extract one score per individual and I wondered if 
anyone knows if there is any methods for doing this?I can fit ID as a 
random effect but I am not sure this changes anything, and is associated 
with the curse of BLUPS. I have been looking for weeks for specific 
information on this but I can not find anyone discussing it.

I can happily provide data if it helps.

Many Thanks
Sam //



-- 

Dr Samantha Patrick
Post Doctoral Fellow
Centre d'Etudes Biologiques de Chizé - CNRS
79360 Villiers-en-Bois
France
T:+33 5 49 09 78 46
M:+33 6 75 06 34 51
Skype: sammy_patrick
http://www.cebc.cnrs.fr/Fidentite/patrick/patrick.htm
http://www.researchgate.net/profile/Samantha_Patrick/


	[[alternative HTML version deleted]]

