Heritability from the estimated components

h2 = var(group) / phenotypic variance, where the phenotypic variance is the sum of the variance AND covariance components of the trait. The denominator is worth stating because it is where this kind of function usually goes quietly wrong: in a direct-maternal model the phenotypic variance is sigma2_a + sigma2_m + sigma_am + sigma2_e (Willham 1972), so the covariance BELONGS in it, and dropping it inflates the ratio. The coefficient 1 on sigma_am is the one in the worked example of OPTION se_covar_function in the BLUPF90 documentation; part of the literature writes 2 sigma_am, which is the variance of the sum a + m of the SAME animal and not the phenotypic variance of a record. Parameters that are not (co)variances stay out of the denominator.

Usage

h2(
  fit,
  group = NULL,
  trait = NULL,
  n = NULL,
  r = 0,
  weight = NULL,
  scale = FALSE
)

Arguments

fit: result of model() or model_mt()
group: covariance group; by default the group of the direct genetic term: the first pedigree or genomic term that is not maternal(), else the first maternal() term, else the first group; the same default in ebv(), solutions(), accuracy(), h2() and h2_curve(). A named group of several terms gives the h2 of its direct term, also when the group has the name of one of its terms; the name of a term outside such a group gives that term's ratio
trait: for a multi-trait fit, which trait; all of them by default
n: group size, read only by a fit with indirect genetic effects and ignored otherwise; a vector gives one value per size
r: average additive relationship between group mates, used with n
weight: for a fit with weights, the weight at which var(residual) is read; by default the harmonic mean of the weights of the records that entered the fit. weight = 1 is the ratio of the components as estimated
scale: FALSE (default) reads every component as estimated, the pedigree components of a fit with metafounders on the base of Gamma; TRUE multiplies those by var_scale, the usual base. Without metafounders it changes nothing. It does not touch the residual, which is read at weight either way: weight = 1 is the ratio of the raw components

Value

a single number, or one per trait in the multi-trait case, named by trait; with a vector n, one per size

Details

Two cases are refused instead of answered. With a reaction norm the heritability is a function of the gradient and not a number, so this sends the caller to h2_curve(). On the observed scale of a threshold trait the ratio is not the liability heritability either; h2_observed() and h2_liability() convert between the two.

No standard error travels with the number here; se_function() gives it by the delta method.

THE MARKER OF THE TERM decides what the ratio is. In a sire model the sire variance is a quarter of the additive variance, so the heritability is 4 var(sire) / Vp (Falconer and Mackay 1996; Mrode and Pocrnic 2023, ch. 3), and that is what is returned for a sire() term. In the sire and maternal-grandsire model, sire(sire, mgs = "mgs"), a record with a known grandsire carries s + m / 2 and one with the grandsire unknown carries s alone, so with the sire and the grandsire unrelated the mean phenotypic variance of the records is (1 + f / 4) var(sire) + var(residual) + others, with f the fraction of the records that entered the fit with a known grandsire (stored in fit$mgs_known), and h2 = 4 var(sire) over that. The relationship between sire and grandsire is ignored there. A term with no relationship matrix (pe(), random()) is refused: var / Vp there is a variance ratio, not a heritability, and the share column of print(fit) gives it.

THE SCALE OF THE COMPONENTS. One rule holds for every term with a relationship matrix, pedigree, genomic, H or metafounders: it enters with its component as estimated, the variance on the base of its matrix (diagonal 1 for a non-inbred base animal). That is the ratio that BLUPF90, ASReml and WOMBAT report from the components. The same number is in h2(), the share column of print(fit) and summary(fit), and the h2 of summary() of a gibbs() fit. A warning gives the other reading when the base of the matrix is far from the records. With metafounders, the pedigree components are on the base of Gamma. When var_scale = 1 + mean(diag(Gamma)) / 2 - mean(Gamma) (Legarra et al. 2015, the factor estimate_gamma() returns) is more than 10% away from 1 and scale is not given, a warning gives the ratio on the usual base, each pedigree component times var_scale; scale = TRUE returns that number (then the h2 of the fit with metafounders is close to the h2 of the fit without them), and scale = FALSE keeps the default and silences the warning. A covariance between two terms takes the geometric mean of the two factors.

WEIGHTS. With weights in model(), var(residual) is the residual variance at weight 1, and a record of weight w has var(residual) / w; weights w and 10 w give the same breeding values and a different var(residual). BLUPF90 (airemlf90), ASReml and WOMBAT report the components, that is, the residual at weight 1, and the ratio from them depends on the scale of the weights. Here the ratio is read at weight, by default the harmonic mean of the weights of the records that entered the fit: var(residual) times mean(1 / w) is the mean residual variance of those records, so the answer does not depend on the scale of the weights. The arithmetic mean weight would put var(residual) / mean(w) in the denominator, which is the residual variance of no record and is below the mean one. weight = 1 gives the ratio of the components as estimated, the number airemlf90 gives (the per-observation heritability when the weights are the numbers of observations averaged in each record). h2() refuses weight = on a fit without weights.

A FIT THAT DID NOT CONVERGE is read at its last iterate, with a warning. When the fit reports the information matrix SINGULAR and the variance of the numerator is in a set of components the data do not separate, the ratio moves along that ridge and is returned as NA, with a warning naming the components; a set that only touches the denominator (permanent environment and residual with one record per animal) leaves the number, with a warning. A fit with no estimates (theta NA) gives NA, saying so. A fit at given components (start with maxiter = 0 and n_em = 0 in model(), start with maxiter = 0 in model_mt()) is read at those components without a word; with n_em > 0 the EM rounds moved the start, and without start the point is the default start, so either fit is read as one that did not converge. In a multi-trait fit one warning covers every trait.

References

Willham, R.L. (1972). The role of maternal effects in animal breeding: III. Biometrical aspects of maternal effects in animals. Journal of Animal Science 35:1288-1293.

Falconer, D.S. & Mackay, T.F.C. (1996). Introduction to Quantitative Genetics, 4th ed. Longman.

Legarra, A., Christensen, O.F., Vitezica, Z.G., Aguilar, I. & Misztal, I. (2015). Ancestral relationships using metafounders: finite ancestral populations and across population relationships. Genetics 200:455-468.