-## vcv.phylo.R (2006-10-04)
+## vcv.phylo.R (2012-02-21)
## Phylogenetic Variance-Covariance or Correlation Matrix
-## Copyright 2002-2006 Emmanuel Paradis
+## Copyright 2002-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-vcv.phylo <- function(phy, model = "Brownian", cor = FALSE)
+vcv <- function(phy, ...) UseMethod("vcv")
+
+vcv.phylo <- function(phy, model = "Brownian", corr = FALSE, ...)
{
- if (class(phy) != "phylo")
- stop('object "phy" is not of class "phylo"')
if (is.null(phy$edge.length))
- stop("the tree has no branch lengths")
-
- foo <- function(node, var, endofclade) {
- ## First, get the extent of clade descending
- ## from `node' in the matrix `edge':
- from <- which(phy$edge[, 1] == node)
- to <- c(from[-1] - 1, endofclade)
- ## Get the #'s of the descendants of `node':
- desc <- phy$edge[from, 2]
- ## The variance of each of these is easy:
- vcv[desc, desc] <<- var + phy$edge.length[from]
- nd <- length(desc)
- ## The variance of `node' is equal to the covariance of
- ## each possible pair among its descendant clades.
- for (i in 1:(nd - 1))
- for (j in (i + 1):nd)
- for (k in phy$edge[from[i]:to[i], 2])
- for (l in phy$edge[from[j]:to[j], 2])
- vcv[k, l] <<- vcv[l, k] <<- var
- for (i in 1:nd) {
- if (desc[i] <= n) next
- foo(desc[i], vcv[desc[i], desc[i]], to[i])
+ stop("the tree has no branch lengths")
+
+ pp <- prop.part(phy)
+
+ phy <- reorder(phy, "pruningwise")
+
+ n <- length(phy$tip.label)
+ e1 <- phy$edge[, 1]
+ e2 <- phy$edge[, 2]
+ EL <- phy$edge.length
+
+ ## xx: vecteur donnant la distance d'un noeud
+ ## ou d'un tip à partir de la racine
+ ## (same than in is.ultrametric)
+ xx <- numeric(n + phy$Nnode)
+
+ vcv <- matrix(0, n, n)
+
+ ## the loop below starts from the bottom of the edge matrix, so
+ ## from the root
+
+ for (i in length(e1):1) {
+ var.cur.node <- xx[e1[i]]
+ xx[e2[i]] <- var.cur.node + EL[i] # sets the variance
+ j <- i - 1L
+ while (e1[j] == e1[i] && j > 0) {
+ left <- if (e2[j] > n) pp[[e2[j] - n]] else e2[j]
+ right <- if (e2[i] > n) pp[[e2[i] - n]] else e2[i]
+ vcv[left, right] <- vcv[right, left] <- var.cur.node # sets the covariance
+ j <- j - 1L
}
}
- n <- length(phy$tip.label)
- n.node <- phy$Nnode
- N <- n.node + n
- vcv <- matrix(0, N, N)
- foo(n + 1, 0, dim(phy$edge)[1])
-
- vcv <- vcv[1:n, 1:n]
- if (cor) {
- ## This is inspired from the code of `cov2cor' (2005-09-08):
- M <- vcv
- Is <- sqrt(1/M[1 + 0:(n - 1)*(n + 1)])
- vcv[] <- Is * M * rep(Is, each = n)
- vcv[1 + 0:(n - 1)*(n + 1)] <- 1
+ diag.elts <- 1 + 0:(n - 1)*(n + 1)
+ vcv[diag.elts] <- xx[1:n]
+
+ if (corr) {
+ ## This is inspired from the code of cov2cor (2005-09-08):
+ Is <- sqrt(1 / vcv[diag.elts])
+ ## below 'vcv[] <- ...' has been changed to 'vcv <- ...'
+ ## which seems to be twice faster for n = 1000 and
+ ## respects the additional attributes (2012-02-21):
+ vcv <- Is * vcv * rep(Is, each = n)
+ vcv[diag.elts] <- 1
}
- rownames(vcv) <- colnames(vcv) <- phy$tip.label
+
+ dimnames(vcv)[1:2] <- list(phy$tip.label)
vcv
}
+
+vcv.corPhyl <- function(phy, corr = FALSE, ...)
+{
+ labels <- attr(phy, "tree")$tip.label
+ dummy.df <- data.frame(1:length(labels), row.names = labels)
+ res <- nlme::corMatrix(Initialize.corPhyl(phy, dummy.df), corr = corr)
+ dimnames(res) <- list(labels, labels)
+ res
+}