Estimate the heritability of a continuous trait in a phylogenetically linked population.

This function estimates the broad-sense heritability of a trait from trait measurements at the tips of a phylogenetic tree connecting the carriers of this trait. By default the function performs a phylogenetic pair (PP) analysis, a phylogenetic mixed model fit (PMM) and a phylogenetic Ornstein-Uhlenbeck mixed model fit (POUMM) and returns an object of class H2Analysis. Use the generic functions summary and plot to see the heritability estimates.

estimateH2(z, tree, methods = list(PP = TRUE, PMM = TRUE, POUMM = TRUE),
  verbose = FALSE)

Arguments

z	A numeric vector of length equal to the number of tips in tree matching their order.
tree	A phylo object.
methods	A named list of logical values and/or lists, specifying each estimating method. The names of the elements should start with PP, POUMM, or PMM. For example, specifying `methods = list(PP = TRUE, POUMM_1 = TRUE, POUMM_2 = list(nSamplesMCMC = 1e6), PMM = TRUE)` will cause the methods PP and PMM to be executed with default settings, while the method POUMM will be executed with a million iterations in each MCMC chain. The values in the list methods should be either logical or lists. To disable the fit of a method, simply exclude it from the list or set it to FALSE.
verbose	Logical: should messages be written on the screen at runtime.

Value

an S3 object of class H2Analysis. Use the generic functions summary and plot to display the results.

Details

The order of elements in z should match the order of the tips in tree. By default the function will run six MCMC chains of length 100,000. Without parallelization of the MCMC chains and the likelihood calculation, this may take more than 20 minutes to run on a tree of 1000 tips.

References

Mitov, V. and Stadler, T. (2016). The heritability of pathogen traits - definitions and estimators. bioRxiv, 058503 doi: https://doi.org/10.1101/058503

Examples

# NOT RUN {
# Please, read the package vignette for more detailed examples.
N <- 200
tr <- ape::rtree(N)
z <- POUMM::rVNodesGivenTreePOUMM(tr, 0, 2, 3, 1, 1)[1:N]

# Enable parallel execution of the POUMM and PMM chains
cluster <- parallel::makeCluster(parallel::detectCores(logical = FALSE))
doParallel::registerDoParallel(cluster)

H2Analysis <- estimateH2(z, tree)

parallel::stopCluster(cluster)

summary(H2Analysis)
# }