sqrt_OU_covariance: (inverse) square root of the phylogenetic covariance

In khabbazian/l1ou: Tools for detecting past changes in the expected mean trait values and studying trait evolution from comparative data

(inverse) square root of the phylogenetic covariance

Description

Computes an inverse square root and square root of the phylogenetic covariance matrix, under the Brownian motion (BM) or the Ornstein-Uhlenbeck (OU) model. The algorithm traverses the tree only once, hence the algorithm is very fast and can be applied to very big trees.

Usage

sqrt_OU_covariance(tree, alpha = 0, root.model = c("OUfixedRoot",
  "OUrandomRoot"), check.order = TRUE, check.ultrametric = TRUE)

Arguments

tree	tree of class phylo with branch lengths. If alpha>0, i.e. under the OU model, the tree has to be ultrametric.
alpha	adaptation rate for the OU model. The default is 0, which corresponds to the BM mode with a fixed ancestral state at the root.
root.model	ancestral state model at the root.
check.order	logical. If TRUE, the order will be checked to be in postorder traversal.
check.ultrametric	logical. If TRUE, the tree will be checked to ultrametric.

Value

sqrtInvSigma	inverse square root of the phylogenetic covariance matrix.
sqrtSigma	square root of the phylogenetic covariance matrix.

References

Mohammad Khabbazian, Ricardo Kriebel, Karl Rohe, and Cécile Ané (2016). "Fast and accurate detection of evolutionary shifts in Ornstein-Uhlenbeck models". Methods in Ecology and Evolution. doi:10.1111/2041-210X.12534

Eric A. Stone. 2011. "Why the phylogenetic regression appears robust to tree misspecification". Systematic Biology, 60(3):245-260.

Examples

data(lizard.tree)
res <- sqrt_OU_covariance(lizard.tree) # alpha not provided: so BM model.
Sigma <- vcv(lizard.tree)
dimnames(Sigma) <- NULL
all.equal(res$sqrtSigma %*% t(res$sqrtSigma), Sigma) # TRUE
all.equal(res$sqrtInvSigma %*% t(res$sqrtInvSigma), solve(Sigma)) # TRUE


##Here's the example from "Eric A. Stone. 2011." (See references)

tr <-  read.tree(text="((((Homo:.21,Pongo:.21):.28,Macaca:.49):.13,Ateles:.62):.38,Galago:1);") 
RE <- sqrt_OU_covariance(tr) 
B <- round( RE$sqrtSigma, digits=3)
D <- round( RE$sqrtInvSigma, digits=3)
print(B)
print(D)


##Here is the examples on how to get the contrasts using sqrt_OU_covariance
data(lizard.tree, lizard.traits)
lizard <- adjust_data(lizard.tree, lizard.traits[,1])
eModel <- estimate_shift_configuration(lizard$tree, lizard$Y)
theta <- eModel$intercept + l1ou:::convert_shifts2regions(eModel$tree,
                             eModel$shift.configuration, eModel$shift.values)
REf <- sqrt_OU_covariance(eModel$tree, alpha=eModel$alpha,
                                         root.model = "OUfixedRoot",
                                         check.order=FALSE, check.ultrametric=FALSE)
 covInverseSqrtf  <- t(REf$sqrtInvSigma)
 covSqrtf   <- REf$sqrtSigma
# `covInverseSqrtf` represents the transpose of square root of  the inverse matrix of covariance for FixedRoot model.
# `covSqrtf` represents the square root of the covariance matrix for FixedRoot model.
 Y  <- rTraitCont(eModel$tree, "OU", theta=theta, 
                                     alpha=eModel$alpha, 
                                     sigma=eModel$sigma, root.value=eModel$intercept)
 contrast    <-  covInverseSqrtf%*%(Y - eModel$mu)
 

khabbazian/l1ou documentation built on Aug. 10, 2022, 7:41 p.m.