est_herit: Estimate heritability with a linear mixed model

Description

Estimate the heritability of a set of traits via a linear mixed model, with possible allowance for covariates.

Usage

est_herit(
  pheno,
  kinship,
  addcovar = NULL,
  weights = NULL,
  reml = TRUE,
  cores = 1,
  ...
)

Value

A vector of estimated heritabilities, corresponding to the columns in pheno. The result has attributes "sample_size", "log10lik" and "resid_sd".

Arguments

pheno: A numeric matrix of phenotypes, individuals x phenotypes.
kinship: A kinship matrix.
addcovar: An optional numeric matrix of additive covariates.
weights: An optional numeric vector of positive weights for the individuals. As with the other inputs, it must have names for individual identifiers.
reml: If true, use REML; otherwise, use maximimum likelihood.
cores: Number of CPU cores to use, for parallel calculations. (If 0, use parallel::detectCores().) Alternatively, this can be links to a set of cluster sockets, as produced by parallel::makeCluster().
...: Additional control parameters (see details).

Details

We fit the model $y = X β + ϵ$ where $ϵ$ is multivariate normal with mean 0 and covariance matrix $σ^{2} [h^{2} (2 K) + I]$ where $K$ is the kinship matrix and $I$ is the identity matrix.

If weights are provided, the covariance matrix becomes $σ^{2} [h^{2} (2 K) + D]$ where $D$ is a diagonal matrix with the reciprocal of the weights.

For each of the inputs, the row names are used as individual identifiers, to align individuals.

If reml=TRUE, restricted maximum likelihood (reml) is used to estimate the heritability, separately for each phenotype.

Additional control parameters include tol for the tolerance for convergence, quiet for controlling whether messages will be display, max_batch for the maximum number of phenotypes in a batch, and check_boundary for whether the 0 and 1 boundary values for the estimated heritability will be checked explicitly.

Examples

Run this code

# read data
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
iron <- iron[,c("19","X")] # subset to two chromosomes

# insert pseudomarkers into map
map <- insert_pseudomarkers(iron$gmap, step=1)

# calculate genotype probabilities
probs <- calc_genoprob(iron, map, error_prob=0.002)

# kinship matrix
kinship <- calc_kinship(probs)

# grab phenotypes and covariates; ensure that covariates have names attribute
pheno <- iron$pheno
covar <- match(iron$covar$sex, c("f", "m")) # make numeric
names(covar) <- rownames(iron$covar)

# perform genome scan
hsq <- est_herit(pheno, kinship, covar)