Cullis et al general heritability in asreml-r (long)

Tena koutou

I’m wishing to calculate the heritability following the general method proposed Cullis et al. (2006, JABES 11:381-393) using R and asreml-r. At a workshop in NZ earlier this year, Brian supplied the following code (I’ve simplified the model a little) to do this based on a bivariate analysis of two traits:

parents <- unique(acpedred$Plant)[16:21]
# Identify parents since it is the family-based heritability that is
# of interest

App <- round(Ared[16:21,16:21], 6)
# App is the numerator relationship matrix (A) for the parents

asr4 <- asreml(cbind(var04,var05)~trait,
               random=~ped(Plant):corh(trait),
               rcov=~units:corh(trait), ...)

Czz3 <- predict(asr4, classify='trait:Plant', maxiter=1,
                levels=list('Plant'=parents, trait='var04'),
                only='ped(Plant):trait', vcov=T)$predictions$vpv
# Using predict to get a variance. levels restricts the predict
# set; only overwrites defaults to give only BLUPs and not general
# mean plus BLUPs; vcov=T will save PEV (prediction error variance
# matrix) and then this is extracted with the $vpv
Czz4 <- predict(asr4, classify='trait:Plant', maxiter=1,
                levels=list('Plant'=parents, trait='var05'),
                only='ped(Plant):trait', vcov=T)$predictions$vpv

1-sum(diag(solve(asr4$gammas[1]*App)%*%Czz3))/6
1-sum(diag(solve(asr4$gammas[2]*App)%*%Czz4))/6

This works fine. If var04 and var05 are the same trait measured in two separate years (on the same perennial plants), it is not unreasonable to want the overall heritability for this trait. Am I correct in assuming this can be obtained simply by taking the average of the variances:

1-sum(diag(solve(mean(asr4$gammas[1:2])*App)%*%Czz3))/6

In my particular example, I have used the as.multivariate way of specifying the multivariate analysis (I have three years of data), and allowed for some spatial auto-correlation:

ASR <- asreml(myVar~year, random=~corh(year):ped(genotype),
              as.multivariate='year',
              rcov=~us(year):id(spRow):ar1(spPlant), ...)
# genotype is equivalent to Plant above

Czz <- predict(ASR, classify='genotype', maxiter=1,
               levels=list('genotype'=parents),
               only='year:ped(genotype)', vcov=T)$predictions$vpv

1-sum(diag(solve(mean(ASR$gammas[2:4])*AParents)%*%Czz))/nrow(Czz)
# Overall heritability. Aparents is equivalent to App above
1-sum(diag(solve(ASR$gammas[2]*AParents)%*%Czz))/nrow(Czz)
# Heritability for the first year
1-sum(diag(solve(ASR$gammas[3]*AParents)%*%Czz))/nrow(Czz)
# Heritability for the second year
1-sum(diag(solve(ASR$gammas[4]*AParents)%*%Czz))/nrow(Czz)
# Heritability for the third year

An alternative way of getting the heritabilities for each year would be to do three separate univariate analyses:

ASR <- asreml(myVar~1, random=~ped(genotype),
              rcov=~id(spRow):ar1(spPlant), ...)
# I’m using version 3 so ped(genotype, var=T) is unnecessary

Czz <- predict(ASR, classify='genotype', maxiter=1,
               levels=list('genotype'=parents),
               only='ped(genotype)', vcov=T)$predictions$vpv

1-sum(diag(solve(summary(ASR)$varcomp[1,2]*AParents)%*%Czz))/nrow(Czz)
# gammas not equivalent to variances in this case and hence it was
# necessary to specify the variances explicitly

When I do this I get the following values:

From multi-variate analysis:
        Overall 0.66
        First year -1.70
        Second year 0.78
        Third year 0.76
From the three univariate analyses:
        First year 0.00
        Second year 0.44
        Third year 0.50

Some difference between the estimates for the individual years is to be expected for a variety of reasons. For example, the multivariate model fits the same autocorrelation between plants within a row for all three years (am I correct that one cannot do otherwise due to lack of separability?). Also, it is a well known phenomenon that the first year of harvest is not well correlated with subsequent years and thus corh(year) might be a poor choice rather than corg(year). However, restricting the multivariate analysis to just the second and third years gave an overall heritability estimate of 0.73, and 0.72 and 0.73 for the second and third years respectively.

So, finally, my main question: Does the magnitude of the difference suggest there is something wrong with my model (or code), or is it not unreasonable? Also, what does the negative estimate for the first year in the multi-variate analysis indicate? (Note: its variance component was not constrained to the boundary, although it was an order of magnitude less than the other two years.)

Thanks to those who have read this far. I’d be interested in your thoughts.

Hei kona ra …

Peter Alspach.

-------------------- m2f --------------------

Sent using Mail2Forum (http://www.mail2forum.com).

Read this topic online here:
http://www.vsni.co.uk/forum/viewtopic.php?p=544#544

-------------------- m2f --------------------
Received on Wed May 19 2009 - 03:47:16 EST