![[Experimental]](figures/lifecycle-experimental.svg)
Details
heritability() is part of the planned v0.1 fitted-object contract. It
works for hsquared_fit objects that contain a Julia result.
What is h2 here? The denominator Vp is the sum of the random-effect
variances plus residual. Fixed-effect variance (sex, age, year, block) is
not in Vp, so the reported h2 is conditional on the fixed effects.
Default animal model:
h2 = Va / VpwithVp = Va + Ve.Repeatability (
target = "repeatability"):h2 = Va / VpwithVp = Va + Vpe + Ve. Permanent environment is in the denominator.repeatability()on the same fit is(Va + Vpe) / Vp.Multi-effect (
target = "multi_effect"):h2for the animal block isVa / Vp, whereVpisVaplus every other random-effect variance plus residual.
Falconer fence for the opt-in two-effect model (target = "two_effect"):
the reported number is the narrow-sense direct heritability
h2 = sigma_a2 / (sigma_a2 + sigma_2 + sigma_e2) within that model (the
additive-genetic variance is divided by the total phenotypic variance,
which now includes the second component sigma_2). The second component's
variance ratio (common-environment c2 or maternal m2) is not a
heritability and is returned separately by common_env_proportion() /
maternal_proportion().
Willham fence for the opt-in direct-maternal correlated model
(target = "direct_maternal"): heritability() returns the labelled
triple – direct h2_d, maternal m2, Willham total h2_T, and r_am – as a
data frame (Willham 1963, 1972). sigma^2_P = sigma_ad + sigma_am + sigma_dm + sigma_e2 = Var(y_i) (coefficient 1 on sigma_dm). A warning is issued
because h2 is denominator-dependent under maternal effects and
h2_T < h2_d is expected when r_am < 0. Use direct_heritability() or
total_heritability() for targeted accessors without the warning.
Reaction-norm fence for the opt-in random-regression model
(target = "random_regression"): heritability is a trajectory h2(t)
over the covariate, not a scalar, so heritability() errors on such a
fit and names the implemented accessor. Use rr_heritability() for the
h2(t) curve.
The bounded genetic GLLVM route does not define heritability; this extractor errors because its G and trait genetic modes are on the link scale.
Non-Gaussian responses (target = "nongaussian", poisson(log) or
binomial(logit)): heritability() returns a table of up to three scales,
not a single number. Those rows follow de Villemereuil, Schielzeth,
Nakagawa and Morrissey (2016, Genetics 204:1281-1294). h2_latent is
Eq 4, V_A / (V_A + V_RE + V_O), and is defined for every family,
including Poisson. The log-link liability variance V_link is 0; that
is not the latent residual V_O, and V_link does not appear in Eq 4.
For logit, report h2_liability = Eq 24,
V_A / (V_A + V_RE + V_O + pi^2/3) (link variance pi^2/3). That
liability scale is not the latent scale. h2_latent equals 1 when
V_RE and V_O are zero (the usual intercept-only animal() model).
Poisson reports h2_latent and a count-scale h2_observation only;
there is no liability scale, and h2_latent is not the
Nakagawa-Schielzeth ln(1 + 1/lambda) convention. binomial(link = "probit") is not available in R. See the Fitting models article for
which row to report.
A non-converged fit still returns the engine number so you can inspect it,
but warns: that number is not an estimate. A near-zero value from a
failed fit is not evidence that heritability is zero. Use
fit_diagnostics() before reading any number. logLik() already refuses
a non-converged fit; heritability() keeps the value and shouts instead.
The same warn-and-return applies to heritability_interval(),
heritability_standard_error(), and
variance_component_standard_errors(). print() says the heritability
is not reportable when the fit did not converge. summary() omits the
printed SE/CI block when the fit did not converge or is at a
variance-component boundary; it does not print those numbers as results.
Testing V_A = 0 and boundary fits: additive variance sits on the edge of
the parameter space (V_A >= 0). The Wald SE and the logit-delta CI are
not valid when V_A is near 0; the usual normal approximation is an
interior-point result. The likelihood-ratio test of V_A = 0 uses a
50:50 mixture of chi^2_0 and chi^2_1 (Self and Liang 1987; Stram and
Lee 1994), so a plain chi^2_1 p-value is twice too large. R does not
run that test: stats::anova() stops, and logLik() / AIC() are
exposed so a reader can build a naive chi^2_1 LRT by hand – do not.
The Julia engine already implements the mixture as
HSquared.nested_lrt(loglik_constrained, loglik_full; df = 1, boundary_df = 1).
That names an existing helper; it is not a new R test statistic.