[R-sig-ME] Fwd: Priors for and estimating heritability from an ordinal and a Poisson animal model

[David Duffy]

On Tue, 25 Sep 2012, Helen Ward wrote:

> Please can someone also tell me whether the heritability estimate from
> the ordinal model means the same thing as a heritability estimated from
> a model in which the residual variance is not fixed?

Looks OK from here, except maybe nu in prior2.  It is probably 
worthwhile simulating some data from what you think is a sensible model 
and running these jobs, as well as tinkering with the priors, and seeing 
what effect that has.  I don't know, for instance, why age at first 
reproduction should be poisson (you could show your conclusions are robust 
to misspecification if the true model is something else).

As a check for "toes", at least, you can calculate the intralitter 
polychoric correlation (I presume you have litters), by producing all
possible sibling pairs, and using the polycor package eg

library(polycor)
litter <- paste(data$sire, data$dam, sep="x")
litter[litter == "NAxNA"] <- NA
# phenotype and sibship indicator
litter.data <- data.frame(litter=litter, trait=data$trait)
litter.data <- litter.data[complete.cases(litter.data),]
# keep those with at least two members
litter.sizes <- as.data.frame(table(litter.data$litter))
useful <- litter.sizes$Var1[litter.sizes$Freq > 1]
litter.data <- litter.data[litter.data$litter %in% useful, ]
# set up for combn() to extract all pairs of siblings
litter.data$litter <- factor(litter.data$litter)
litter.data$trait <- as.character(litter.data$trait)
allpairs <- sapply(split(litter.data$trait, litter.data$litter),
                    function(x) combn(x,2))
allpairs <- as.data.frame(matrix(unlist(unlist(allpairs)), nc=2, byrow=T))
# finally call polychor() to estimate sibling intraclass correlation
names(allpairs) <- c("sib1", "sib2")
polychor(c(allpairs$sib1, allpairs$sib2), c(allpairs$sib2, allpairs$sib1))


You could extend this to other classes of relative and make it bivariate.

Just 2c, David Duffy.