Title:	Heritability Estimation from Mixed Models
Version:	0.1.0
Description:	Reporting heritability estimates is an important to quantitative genetics studies and breeding experiments. Here we provide functions to calculate various broad-sense heritabilities from 'asreml' and 'lme4' model objects. All methods we have implemented in this package have extensively discussed in the article by Schmidt et al. (2019) <doi:10.1534/genetics.119.302134>.
License:	GPL (≥ 3)
Imports:	cli, emmeans, Matrix, methods, stringr, vctrs
Suggests:	testthat (≥ 3.0.0), agridat, knitr, rmarkdown, lme4, pbkrtest, dplyr, ggplot2, tidyr, purrr, here
SystemRequirements:	asreml (4.2.0)
Config/testthat/edition:	3
Encoding:	UTF-8
RoxygenNote:	7.3.2
Depends:	R (≥ 4.1.0)
LazyData:	true
VignetteBuilder:	knitr
URL:	https://anu-aagi.github.io/heritable/
Language:	en
NeedsCompilation:	no
Packaged:	2026-01-08 06:34:59 UTC; fontikar
Author:	Fonti Kar [aut, cre], Yidi Deng [aut], Emi Tanaka [aut, cph]
Maintainer:	Fonti Kar <fonti.kar@anu.edu.au>
Repository:	CRAN
Date/Publication:	2026-01-09 19:00:02 UTC

Calculate heritability from mixed model object

Description

A case-specific wrapper for calculating heritability.

• The upper case prefix H2_ refers to the wrapper or subfunctions e.g. H2_Delta() for calculating broad sense heritability

Usage

H2(model, target, method = c("Cullis", "Oakey", "Delta", "Piepho", "Standard"), options)

Arguments

Details

The following methods are currently implemented for broad-sense heritability H2(method = "XX"):

• "Cullis":

  H^2_{Cullis} = 1 - \frac{PEV^{BLUP}_{\overline\Delta ij}}{2\sigma^2_g}

• "Oakey":

  H^2_{Oakey} = \frac{\sum_{i = n_z+1}^{n_g} \lambda_i}{\sum_{n_g}^{\lambda_i\neq 0}}

• "Delta":

  H^2_{\Delta ..} = 1 - \frac{PEV^{BLUP}_{\overline\Delta ..}}{2\sigma^2_g}

• "Piepho":

  H^2_{Piepho} = \frac{\sigma^2_g}{\sigma^2_g + \overline{PEV_{BLUE_g}} / 2}

• "Standard":

  H^2_{Standard} = \frac{\sigma^2_g}{\sigma^2_g + \frac{1}{n_g}\sum_{n_g}^{i=1} \sigma^2_p / n_{gi}}

For further details of a specific method - take a look at help file for each subfunctions ?H2_Cullis

Value

A named numeric vector, length matching number of methods supplied

References

• Cullis, B. R., Smith, A. B., & Coombes, N. E. (2006). On the design of early generation variety trials with correlated data. Journal of Agricultural, Biological, and Environmental Statistics, 11(4), 381–393. https://doi.org/10.1198/108571106X154443

• Oakey, H., Verbyla, A., Pitchford, W., Cullis, B., & Kuchel, H. (2006). Joint modeling of additive and non-additive genetic line effects in single field trials. Theoretical and Applied Genetics, 113(5), 809–819. https://doi.org/10.1007/s00122-006-0333-z

• Schmidt, P., Hartung, J., Rath, J., & Piepho, H.-P. (2019). Estimating Broad-Sense Heritability with Unbalanced Data from Agricultural Cultivar Trials. Crop Science, 59(2), 525–536. https://doi.org/10.2135/cropsci2018.06.0376

• Piepho, H.-P., & Möhring, J. (2007). Computing Heritability and Selection Response From Unbalanced Plant Breeding Trials. Genetics, 177(3), 1881–1888. https://doi.org/10.1534/genetics.107.074229

• Falconer, D. S., & Mackay, T. F. C. (1996). Introduction to quantitative genetics (4th ed.). Longman.

See Also

H2_Cullis(), H2_Oakey(), H2_Delta(), H2_Piepho(), H2_Standard()

Examples

# lme4 model
lettuce_subset <- lettuce_phenotypes |> subset(loc == "L2")
lettuce_lme4 <- lme4::lmer(y ~ rep + (1 | gen), data = lettuce_subset)
H2(lettuce_lme4, target = "gen", method = c("Standard", "Delta"))

# asreml model (Requires license)
## Not run: 
lettuce_asreml <- asreml::asreml(fixed = y ~ rep,
                                 random = ~ gen,
                                 data = lettuce_subset,
                                 trace = FALSE
                                 )

H2(lettuce_asreml, target = "gen", method = c("Standard", "Delta"))

## End(Not run)

Calculate Cullis' heritability from model object

Description

Compute "generalised heritability" for unbalanced experimental designs. See Cullis, Smith and Coombes (2006) for derivation.

Usage

H2_Cullis(model, target, options)

Arguments

Details

The equation for Cullis heritability is as follow

H^2_{Cullis} = 1 - \frac{PEV^{BLUP}_{\overline\Delta ij}}{2\sigma^2_g}

where:

• PEV is the prediction error variance matrix of the pairwise differences among BLUPS

• \sigma^2 is the variance attributed to differences between genotype

Value

Numeric value

References

Cullis, B. R., Smith, A. B., & Coombes, N. E. (2006). On the design of early generation variety trials with correlated data. Journal of Agricultural, Biological, and Environmental Statistics, 11(4), 381–393. https://doi.org/10.1198/108571106X154443

Examples

# lme4 model
lettuce_subset <- lettuce_phenotypes |> subset(loc == "L2")
lettuce_lme4 <- lme4::lmer(y ~ rep + (1 | gen), data = lettuce_subset)
H2_Cullis(lettuce_lme4, target = "gen")

# asreml model (Requires license)
## Not run: 
lettuce_asreml <- asreml::asreml(fixed = y ~ rep,
                                 random = ~ gen,
                                 data = lettuce_subset,
                                 trace = FALSE
                                 )

H2_Cullis(lettuce_asreml, target = "gen")

## End(Not run)

Calculate Cullis heritability using variance parameters

Description

Compute the Cullis heritability for genotype means using the average variance of pairwise differences of best linear unbiased predictors (BLUPs).

Usage

H2_Cullis_parameters(vd_BLUP_avg, vc_g)

Arguments

Details

The equation for Cullis heritability is as follow

H^2_{Cullis} = 1 - \frac{PEV^{BLUP}_{\overline\Delta ij}}{2\sigma^2_g}

where:

• PEV is the prediction error variance matrix of the pairwise differences among BLUPS

• \sigma^2 is the variance attributed to differences between genotype

Value

Numeric value

References

Cullis, B. R., Smith, A. B., & Coombes, N. E. (2006). On the design of early generation variety trials with correlated data. Journal of Agricultural, Biological, and Environmental Statistics, 11(4), 381–393. https://doi.org/10.1198/108571106X154443

Examples

H2_Cullis_parameters(vd_BLUP_avg = 0.25, vc_g = 0.8)

Calculate average heritability of differences between genotypes from model object

Description

Instead of computing heritability on a "entry-mean" basis, this method calculates heritability using "entry-differences". Entry here is referring to the genotype, line or variety of interest. See reference for origin and interpretation of H2_Delta and it's variants

Usage

H2_Delta(model,
         target,
         type = c("BLUP", "BLUE"),
         aggregate = c("arithmetic", "harmonic"),
         options
         )

Arguments

Details

The heritability of differences between genotypes is given by:

H^2_{\Delta ..} = 1 - \frac{PEV^{BLUP}_{\overline\Delta ..}}{2\sigma^2_g}

where:

• PEV^{BLUP}_{\overline\Delta ..} is the mean of the prediction error variance matrix for the pairwise differences among BLUPs (BLUEs if method = "BLUE") across all genotypes

• \sigma^2 is the variance attributed to differences between genotype

See reference page 995 - 997 for full derivation of this heritability measure and related variants

Value

Numeric

References

Schmidt, P., Hartung, J., Rath, J., & Piepho, H.-P. (2019). Estimating Broad-Sense Heritability with Unbalanced Data from Agricultural Cultivar Trials. Crop Science, 59(2), 525–536. https://doi.org/10.2135/cropsci2018.06.0376

See Also

H2_Delta_by_genotype(), H2_Delta_pairwise()

Examples

# lme4 model
lettuce_subset <- lettuce_phenotypes |> subset(loc == "L2")
lettuce_lme4 <- lme4::lmer(y ~ rep + (1 | gen), data = lettuce_subset)
H2_Delta(lettuce_lme4, target = "gen", type = "BLUP")

# asreml model (Requires license)
## Not run: 
lettuce_asreml <- asreml::asreml(fixed = y ~ rep,
                                 random = ~ gen,
                                 data = lettuce_subset,
                                 trace = FALSE
                                 )

H2_Delta(lettuce_asreml, target = "gen", type = "BLUP")

## End(Not run)

Calculate heritability of differences for a given genotype from model object

Description

Instead of computing heritability on a "entry-mean" basis, this method calculates heritability using "entry-differences". Entry here is referring to the genotype, line or variety of interest. See reference for origin and interpretation of h2/H2_Delta_by_genotype and it's variants

Usage

H2_Delta_by_genotype(model, target, type = c("BLUE", "BLUP"), options)

Arguments

Details

The heritability of differences for a given genotype is given by:

H^2_{\Delta i.} = 1 - \frac{PEV^{BLUP}_{\overline\Delta i.}}{2\sigma^2_g}

where:

• PEV^{BLUP}_{\overline\Delta i.} is the arithmetic mean of the prediction error variance matrix for pairwise differences among BLUPs (or BLUEs if method = "BLUE") for a given genotype

• \sigma^2 is the variance attributed to differences between genotype

See reference page 995 - 997 for full derivation of this heritability measure and related variants

Value

Numeric

Named list, with each element containing a named numeric vector

References

Schmidt, P., Hartung, J., Rath, J., & Piepho, H.-P. (2019). Estimating Broad-Sense Heritability with Unbalanced Data from Agricultural Cultivar Trials. Crop Science, 59(2), 525–536. https://doi.org/10.2135/cropsci2018.06.0376

See Also

H2_Delta(), H2_Delta_pairwise()

Examples

# lme4 model
lettuce_subset <- lettuce_phenotypes |> subset(loc == "L2")
lettuce_lme4 <- lme4::lmer(y ~ rep + (1 | gen), data = lettuce_subset)
H2_Delta_by_genotype(lettuce_lme4, target = "gen", type = "BLUP")

# asreml model (Requires license)
## Not run: 
lettuce_asreml <- asreml::asreml(fixed = y ~ rep,
                                 random = ~ gen,
                                 data = lettuce_subset,
                                 trace = FALSE
                                 )

H2_Delta_by_genotype(lettuce_asreml, target = "gen", type = "BLUP")

## End(Not run)

Calculate pairwise heritability of differences between genotypes from model object

Description

Instead of computing heritability on a "entry-mean" basis, this method calculates heritability using "entry-differences". Entry here is referring to the genotype, line or variety of interest. See reference for origin and interpretation of h2/H2_Delta_pairwise and it's variants

Usage

H2_Delta_pairwise(model, target, type = c("BLUE", "BLUP"), options)

Arguments

Value

A dspMatrix

References

Schmidt, P., Hartung, J., Rath, J., & Piepho, H.-P. (2019). Estimating Broad-Sense Heritability with Unbalanced Data from Agricultural Cultivar Trials. Crop Science, 59(2), 525–536. https://doi.org/10.2135/cropsci2018.06.0376

See Also

H2_Delta_by_genotype(), H2_Delta()

Examples

# lme4 model
lettuce_subset <- lettuce_phenotypes |> subset(loc == "L2")
lettuce_lme4 <- lme4::lmer(y ~ rep + (1 | gen), data = lettuce_subset)
H2_Delta_pairwise(lettuce_lme4, target = "gen", type = "BLUP")

# asreml model (Requires license)
## Not run: 
lettuce_asreml <- asreml::asreml(fixed = y ~ rep,
                                 random = ~ gen,
                                 data = lettuce_subset,
                                 trace = FALSE
                                 )

H2_Delta_pairwise(lettuce_asreml, target = "gen", type = "BLUP")

## End(Not run)

Calculate heritability of pairwise differences using variance parameters

Description

Compute broad-sense heritability of differences using the variance of differences between two BLUPs/BLUEs

Usage

h2_Delta_parameters(G_g, vd_matrix, type)

H2_Delta_parameters(vc_g, vd_matrix, type)

Arguments

Details

See H2_Delta() and reference for full derivation and equation for heritability Delta

Value

Matrix of pairwise heritability of differences among BLUES or BLUPs

References

Schmidt, P., Hartung, J., Rath, J., & Piepho, H.-P. (2019). Estimating Broad-Sense Heritability with Unbalanced Data from Agricultural Cultivar Trials. Crop Science, 59(2), 525–536. https://doi.org/10.2135/cropsci2018.06.0376

Examples

h2_Delta_parameters(G_g = diag(0.15, 2, 2), vd_matrix = matrix(c(NA,0.2,0.2,NA),2,2), type = "BLUP")

H2_Delta_parameters(vc_g = 0.01, vd_matrix = matrix(c(NA,0.2,0.2,NA),2,2), "BLUE")

Calculate Oakey's heritability from model object

Description

Compute heritability for genotype means using the variance–covariance matrix of the genotype BLUPs as described by Oakey et al. (2006).

Usage

H2_Oakey(model, target, options)

Arguments

Details

H^2_{Oakey} = \frac{\sum_{i = n_z+1}^{n_g} \lambda_i}{\sum_{n_g}^{\lambda_i\neq 0}}

where:

• n_g is the number of genotypes

• n_z is the number of zero eigenvalues

• \lambda_i is the ith eigenvalue of the matrix I_{m} - G^{-1}C^{gg}

• \sigma^2 is the variance attributed to differences between genotype

See pages 813 and 818 of the reference for full derivation and explanation for Oakey's heritability

Value

Numeric

References

Oakey, H., Verbyla, A., Pitchford, W., Cullis, B., & Kuchel, H. (2006). Joint modeling of additive and non-additive genetic line effects in single field trials. Theoretical and Applied Genetics, 113(5), 809–819. https://doi.org/10.1007/s00122-006-0333-z

Examples

# lme4 model
lettuce_subset <- lettuce_phenotypes |> subset(loc == "L2")
lettuce_lme4 <- lme4::lmer(y ~ rep + (1 | gen), data = lettuce_subset)
H2_Oakey(lettuce_lme4, target = "gen")

# asreml model (Requires license)
## Not run: 
lettuce_asreml <- asreml::asreml(fixed = y ~ rep,
                                 random = ~ gen,
                                 data = lettuce_subset,
                                 trace = FALSE
                                 )

H2_Oakey(lettuce_asreml, target = "gen")

## End(Not run)

Calculate Oakey's heritability using variance parameters

Description

Rather than providing a model object, supply the necessary components to compute this heritability measure.

Usage

H2_Oakey_parameters(Gg_inv, C_gg)

Arguments

Value

Numeric value

Examples

Gg_inv = diag(1/0.15, 3, 3)
C_gg <- matrix(
  c(
    0.08, 0.01, 0.00,
    0.01, 0.07, 0.01,
    0.00, 0.01, 0.09
  ),
  nrow = 3, byrow = TRUE
)
H2_Oakey_parameters(Gg_inv, C_gg)

Calculate Piepho's heritability from model object Compute Piepho's heritability using variance differences between genotype BLUEs

Description

Calculate Piepho's heritability from model object Compute Piepho's heritability using variance differences between genotype BLUEs

Usage

H2_Piepho(model, target, options)

Arguments

Details

The equation for Piepho's heritability is as follows:

H^2_{Piepho} = \frac{\sigma^2_g}{\sigma^2_g + \overline{PEV_{BLUE_g}} / 2}

where:

• \overline{PEV_{BLUE_g}} is the prediction error variance matrix for genotype BLUEs

• \sigma^2_g is the variance attributed to differences between genotype

See reference for full derivation and details.

Value

Numeric

References

Piepho, H.-P., & Möhring, J. (2007). Computing Heritability and Selection Response From Unbalanced Plant Breeding Trials. Genetics, 177(3), 1881–1888. https://doi.org/10.1534/genetics.107.074229

Examples

# lme4 model
lettuce_subset <- lettuce_phenotypes |> subset(loc == "L2")
lettuce_lme4 <- lme4::lmer(y ~ rep + (1 | gen), data = lettuce_subset)
H2_Piepho(lettuce_lme4, target = "gen")

# asreml model (Requires license)
## Not run: 
lettuce_asreml <- asreml::asreml(fixed = y ~ rep,
                                 random = ~ gen,
                                 data = lettuce_subset,
                                 trace = FALSE
                                 )

H2_Piepho(lettuce_asreml, target = "gen")

## End(Not run)

Calculate Piepho's heritability using variance parameters

Description

Compute Piepho's heritability using the variance of differences between two BLUES.

Usage

H2_Piepho_parameters(vc_g, vd_BLUE_avg)

Arguments

Details

The equation for Piepho's heritability is as follows:

H^2_{Piepho} = \frac{\sigma^2_g}{\sigma^2_g + \overline{PEV_{BLUE_g}} / 2}

where:

• \overline{PEV_{BLUE_g}} is the prediction error variance matrix for genotype BLUEs

• \sigma^2_g is the variance attributed to differences between genotype

Value

Numeric value

References

Piepho, H.-P., & Möhring, J. (2007). Computing Heritability and Selection Response From Unbalanced Plant Breeding Trials. Genetics, 177(3), 1881–1888. https://doi.org/10.1534/genetics.107.074229

Examples

H2_Piepho_parameters(vc_g = 0.25, vd_BLUE_avg = 0.68)

Calculate standard heritability from model object

Description

Compute standard heritability using the classic ratio method of genotypic and phenotypic variance. See Falconer & Mackay (1996)

Usage

H2_Standard(model, target, options)

Arguments

Details

The equation used to calculate standard heritability is:

H^2_{Standard} = \frac{\sigma^2_g}{\sigma^2_g + \frac{1}{n_g}\sum_{n_g}^{i=1} \sigma^2_p / n_{gi}}

where:

• n_g is the number of genotypes

• n_{gi} is the number of replicate for a given genotype i

• \sigma_g is the variance attributed to genotype differences

• \sigma_p is the variance attributed to phenotypic differences

Value

Numeric value

References

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

Examples

# lme4 model
lettuce_subset <- lettuce_phenotypes |> subset(loc == "L2")
lettuce_lme4 <- lme4::lmer(y ~ rep + (1 | gen), data = lettuce_subset)
H2_Standard(lettuce_lme4, target = "gen")

# asreml model (Requires license)
## Not run: 
lettuce_asreml <- asreml::asreml(fixed = y ~ rep,
                                 random = ~ gen,
                                 data = lettuce_subset,
                                 trace = FALSE
                                 )

H2_Standard(lettuce_asreml, target = "gen")

## End(Not run)

Calculate Standard heritability using variance parameters

Description

Compute Standard heritability for genotype means using the variance components of genotype and residuals.

Usage

H2_Standard_parameters(vc_g, vc_e, n_r = 1)

Arguments

Details

The equation for Standard heritability is as follows:

H^2_{Standard} = \frac{\sigma^2_g}{\sigma^2_g + \frac{1}{n_g}\sum_{n_g}^{i=1} \sigma^2_p / n_{gi}}

where:

• n_g is the number of genotypes

• n_{gi} is the number of replicate for a given genotype i

• \sigma_g is the variance attributed to genotype differences

• \sigma_p is the variance attributed to phenotypic differences

Value

Numeric value

References

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

Examples

H2_Standard_parameters(vc_g = 0.25, vc_e = 0.8)

Molecular marker data and genomic relatedness matrix of 89 lettuce varieties

Description

Molecular marker data and genomic relatedness matrix of 89 lettuce varieties

Usage

lettuce_markers

lettuce_GRM

Format

lettuce_markers

A data frame with 89 rows and 301 columns:

• gen genotype identifier

• 300 genetic markers scored as -1, 0, 1 (see Details)

lettuce_GRM

A matrix array with 89 rows and 89 columns where each row/column represents a genotype

Details

The varieties were genotyped with a total of 300 markers (i.e. 95 single nucleotide polymorphisms and 205 amplified fragment length polymorphism markers, see Hayes et al. (2014) for more details of marker matrix. The biallelic marker M_iw for the ith genotype and the wth marker with alleles A_1 (i.e. the reference allele) and A_2 was coded as:

• 1 for A_1 A_1,

• -1 for A_2 A_2

• 0 for A_1 A_2 and A_2 A_1

Source

https://figshare.com/articles/dataset/Lettuce_trial_phenotypic_and_marker_data_/8299493

References

Hadasch, S., Simko, I., Hayes, R.J., Ogutu, J.O. and Piepho, H.-P. (2016), Comparing the Predictive Abilities of Phenotypic and Marker-Assisted Selection Methods in a Biparental Lettuce Population. The Plant Genome, 9: plantgenome2015.03.0014. doi:10.3835/plantgenome2015.03.0014

Hayes, R. J., Galeano, C. H., Luo, Y., Antonise, R., & Simko, I. (2014). Inheritance of Decay of Fresh-cut Lettuce in a Recombinant Inbred Line Population from ‘Salinas 88’ × ‘La Brillante’. Journal of the American Society for Horticultural Science, 139(4), 388–398. doi:10.21273/JASHS.139.4.388

Phenotypic of 89 lettuce varieties

Description

89 lettuce varieties tested at three environments, each laid out as a randomized complete block design. The measured trait was resistance to downy mildew scored on a scale ranging from 0 to 5.

Usage

lettuce_phenotypes

Format

lettuce_phenotypes

A data frame with 703 rows and 4 columns:

• loc environment identifier

• gen genotype identifier

• rep replicate identifier

• y resistance to downy mildew scored on a scale ranging from 0 to 5

Pull fixed and random terms from a model formula

Description

Extract the labels of fixed and random terms from a model object that exposes a formula with fixed and random components (for example objects produced by asreml::asreml). The function returns a named list containing two character vectors: fixed and random.

Usage

## S3 method for class 'asreml'
pull_terms(model)

Arguments

Value

A named list with components: