sqrt_OU_covariance: (inverse) square root of the phylogenetic covariance

Description Usage Arguments Value References Examples

View source: R/sqrt_OU_covariance.R

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

1
2
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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 July 30, 2018, 12:05 p.m.