1 ## vcv.phylo.R (2006-10-04)
3 ## Phylogenetic Variance-Covariance or Correlation Matrix
5 ## Copyright 2002-2006 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 vcv.phylo <- function(phy, model = "Brownian", cor = FALSE)
12 if (class(phy) != "phylo")
13 stop('object "phy" is not of class "phylo"')
14 if (is.null(phy$edge.length))
15 stop("the tree has no branch lengths")
17 foo <- function(node, var, endofclade) {
18 ## First, get the extent of clade descending
19 ## from `node' in the matrix `edge':
20 from <- which(phy$edge[, 1] == node)
21 to <- c(from[-1] - 1, endofclade)
22 ## Get the #'s of the descendants of `node':
23 desc <- phy$edge[from, 2]
24 ## The variance of each of these is easy:
25 vcv[desc, desc] <<- var + phy$edge.length[from]
27 ## The variance of `node' is equal to the covariance of
28 ## each possible pair among its descendant clades.
31 for (k in phy$edge[from[i]:to[i], 2])
32 for (l in phy$edge[from[j]:to[j], 2])
33 vcv[k, l] <<- vcv[l, k] <<- var
35 if (desc[i] <= n) next
36 foo(desc[i], vcv[desc[i], desc[i]], to[i])
40 n <- length(phy$tip.label)
43 vcv <- matrix(0, N, N)
44 foo(n + 1, 0, dim(phy$edge)[1])
48 ## This is inspired from the code of `cov2cor' (2005-09-08):
50 Is <- sqrt(1/M[1 + 0:(n - 1)*(n + 1)])
51 vcv[] <- Is * M * rep(Is, each = n)
52 vcv[1 + 0:(n - 1)*(n + 1)] <- 1
54 rownames(vcv) <- colnames(vcv) <- phy$tip.label