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, such as the rho(residual) of model_ar1(), 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(), model_mt() or model_ar1()
group: covariance group; by default the group of the direct genetic term: the first pedigree or genomic term that is not maternal() or indirect(), else the first kernel() with an estimated variance over the animals (most of its levels are animals of the pedigree, or there is no pedigree term), else the first maternal() term, else the first other kernel() with an estimated variance, 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, required when the model has an indirect() term 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

With an indirect() term the phenotypic variance depends on the group size: each record carries the indirect effects of its n - 1 group mates. The denominator is then the one of Bijma, Muir and Van Arendonk (2007), sigma2_AD + c_V [1 + (n - 2) r] sigma2_AS + 2 r c_D sigma_ADS plus the other components, with r the average relationship between group mates and, under the dilution d of the term (Bijma 2010), c_D = (n - 1)^(1 - d) and c_V = (n - 1)^(1 - 2d). That is why n is required there. t2() gives the total heritable variance of the same model.

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, and t2() returns it for the indirect-effect model.

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, metafounders or a declared kernel() K: 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 BLUPF90, ASReml and WOMBAT report from the components, and the same likelihood gives the same number: animal(id) and kernel(id, K = A) with the A of the same pedigree, or animal(id) with genotypes = and kernel(id, K = H), give the same h2, in h2(), t2(), the share column of print(fit) and summary(fit), and the h2 of summary() of a gibbs() fit. Two warnings give the other reading when the base of a matrix is far from the records:

a declared K whose mean diagonal over the records that entered the fit is more than 10% away from 1 (a partial multibreed matrix, whose diagonal is about the breed fraction, or an unscaled kernel) gives a warning with the mean diagonal and the value the ratio would take with each such K at its mean diagonal over the records; scaling K to mean diagonal 1 (K / mean(diag(K))) makes the two agree;
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.

Two kinds of kernel() are not relationships but covariances of the records, and enter at the variance they put on a record:

a kernel() that is the associative_matrix() of the indirect() term of the same fit (recognised by its diagonal, (n_i - 1)^(1 - 2d) on every record, and, when that diagonal is 1 everywhere, with pens of 2 or d = 0.5, by the attribute br_associative that associative_matrix() sets): at (n - 1)^(1 - 2d) times its component for the group size n asked for, or at its mean diagonal over the records without n (the share column);
a kernel(fixed = v) term, the known covariance of a two-stage analysis (the sampling variances of the first-stage estimates): at v times the mean diagonal of its K over the records, the mean sampling variance, and not at v.

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). A fit made by indirect_residual() reads its residual at the group size: 1 + (n - 1)^(1 - 2d) k times var(residual). weight = on a fit without weights is refused.

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() and model_ar1()) 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.

Bijma, P., Muir, W.M. & Van Arendonk, J.A.M. (2007). Multilevel selection 1: quantitative genetics of inheritance and response to selection. Genetics 175:277-288.

Bijma, P. (2010). Multilevel selection 4: modeling the relationship of indirect genetic effects and group size. Genetics 186:1029-1031.