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

David Duffy David.Duffy at qimr.edu.au
Wed Sep 26 01:52:24 CEST 2012

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

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.

More information about the R-sig-mixed-models mailing list