The Standard Error of Heritability

Exactly what does ASReml compute?

The expressions used for calculating the heritability and its standard error are described in chapter 13 of the User Guide. The detail that is implied is the calculation of variances of the components.

Consider an example sent by another user for a bivariate analysis of a repeated trait so there is genetic variance and betweed and within animals. The variance components reported in the .asr file were
 Source                Model  terms     Gamma     Component    Comp/SE   % C
 Residual            UnStruct     1   1.26730       1.26730     108.05   0 P
 Residual            UnStruct     1  0.610439      0.610439      33.58   0 P
 Residual            UnStruct     2   4.77802       4.77802     110.16   0 P
 Residual            UnStruct     1  0.394353E-01  0.394353E-01   3.88   0 P
 Residual            UnStruct     2 -0.274366     -0.274366     -12.47   0 P
 Residual            UnStruct     3   1.84899       1.84899     108.08   0 P
 Trait.id            UnStruct     1  0.303319E-01  0.303319E-01   4.05   0 P
 Trait.id            UnStruct     1  0.277637E-01  0.277637E-01   2.50   0 P
 Trait.id            UnStruct     2  0.158267      0.158267       5.25   0 P
 Trait.id            UnStruct     1  0.349000E-02  0.349000E-02   0.27   0 P
 Trait.id            UnStruct     2 -0.905091E-01 -0.905091E-01  -3.51   0 P
 Trait.id            UnStruct     3  0.417714      0.417714      10.82   0 P
 Trait.ide(id)       UnStruct     1  0.562876E-01  0.562876E-01   6.72   0 P
 Trait.ide(id)       UnStruct     1  0.248613E-01  0.248613E-01   2.08   0 P
 Trait.ide(id)       UnStruct     2  0.242284      0.242284       7.56   0 P
 Trait.ide(id)       UnStruct     1 -0.735900E-02 -0.735900E-02  -0.70   0 P
 Trait.ide(id)       UnStruct     2 -0.439076E-01 -0.439076E-01  -2.13   0 P
 Trait.ide(id)       UnStruct     3  0.224123      0.224123       8.64   0 P
 Covariance/Variance/Correlation Matrix UnStructured
   1.267     0.2481     0.2576E-01
  0.6104      4.778    -0.9231E-01
  0.3944E-01-0.2744      1.849
 Covariance/Variance/Correlation Matrix UnStructured
  0.3033E-01 0.4007     0.3101E-01
  0.2776E-01 0.1583    -0.3520
  0.3490E-02-0.9051E-01 0.4177
 Covariance/Variance/Correlation Matrix UnStructured
  0.5629E-01 0.2129    -0.6552E-01
  0.2486E-01 0.2423    -0.1884
 -0.7359E-02-0.4391E-01 0.2241
The variance matrix for these 18 parameters is printed in the .vvp file. This variance matrix is simply the inverse of the Average Information matrix.

To calculate the heritability for the first trait, we first need the phenotypic variance. This is produced in the .pvc file by ASReml using a .pin file containing
 F Phen 1 7 13
 H H2a  7 19
The components are for the first trait are
   1 Residual                 1.26730
   7 Trait.id                0.303319E-01
  13 Trait.ide(id)           0.562876E-01
The SUM is 1.35392.

The variances for the components (from rows 1, 7 and 13 of the vvp file are)
    0.137567E-03
    0.265106E-06     0.561911E-04
   -0.278960E-04    -0.381871E-04      0.700629E-04
Making it square and adjusting the power
    1.375670E-04     0.002651E-04     -0.278960E-04
    0.002651E-04     0.561911E-04     -0.381871E-04
   -0.278960E-04    -0.381871E-04      0.700629E-04
Now to calculate the heritability, we must first compute the phenotypic variance (the sum of the three components 1.26730 + 0.030332 + 0.056288). The variance of this value is obtained by adding the rows, then the columns gives of the variance matrix to obtain.
    1.375670E-04     0.002651E-04     -0.278960E-04     1.099361E-04
    0.002651E-04     0.561911E-04     -0.381871E-04     0.182691E-04
   -0.278960E-04    -0.381871E-04      0.700629E-04     0.039798E-04
    1.099361E-04     0.182691E-04      0.039798E-04     1.321850E-04
This gives the SE of the phenotypic variance as 1.14972E-2 which agrees with
      19  Phen  1       1.354      0.1150E-01
The heritability is 0.03033/1.35392 = 0.0224

Now var(heritability) is
Var(Sn/Sd) = (Sn/Sd)2 (Var(Sn)/(Sn)2 + Var(Sd)/(Sd)2 -2Cov(Sn,Sd)/Sn/Sd)
= ( 0.0224)2 ( 0.561911E-04/(0.03033)2 + 1.321850E-04/(1.35392)2
- 2 (0.182691E-04)/0.03033/1.35392)
= .00050176 (.061083 + 0.000072 - 0.000890)
= 0.000029796
SE(h2)=0.0054586 which again agrees with what was in the .pvc file.
          H2a         = Trait.id   7/ Phen  1  19=          0.0224    0.0055

How good is this estimate?

This estimate is an approximation and appears to be a little inflated. In the component table above is the field headed Comp/SE. This contains the ratio of the component with its standard error (the square root of the variance value reported in the .vvp file. A value of 1 for this ratio is often significant when tested with a likelihood ratio test (LRT) (the exact value correstonding to significant usually ranges from 0.7 to 1.4 depending on the problem). Therefore, I use LRT for actually testing variance parameters. This however implies that the SE based on the AI matrix is inflated. Nevertheless, it is useful as an indicator of precision of the estimate.

Univariate models.

In a univariate analysis, the analysis is performed using the variance ratio but is reported as a variance. In this case, the reported variance of the variance parameter has been derived from the inverse AI matrix using standard first order expressions for the variance of a product (ratio times residual).

Boundary models.

When constraints apply to the variance parameter (e.g. positive definite) but the AI algorithm would produce a value outside the parameter space, the corresponding rows of the AI matrix are zeroed and standard errors will either be not available or conditional on the corresponding parameters being 'known' (i.e. no variance).

Arthur Gilmour 7 Nov 2008

See Also