est_herit {qtl2}

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)