est_herit: Estimate heritability with a linear mixed model

In qtl2: Quantitative Trait Locus Mapping in Experimental Crosses

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,
  ...
)

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 \beta + \epsilon where \epsilon is multivariate normal with mean 0 and covariance matrix \sigma^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 \sigma^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.

Value

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

Examples

# read data
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))


# 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)

qtl2 documentation built on May 29, 2024, 11:46 a.m.