R/vcv.phylo.R

   1 ## vcv.phylo.R (2009-11-19)
   2
   3 ##   Phylogenetic Variance-Covariance or Correlation Matrix
   4
   5 ## Copyright 2002-2009 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 vcv <- function(phy, ...) UseMethod("vcv")
  11
  12 vcv.phylo <- function(phy, model = "Brownian", corr = FALSE, ...)
  13 {
  14     if (is.null(phy$edge.length))
  15       stop("the tree has no branch lengths")
  16
  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]
  26         nd <- length(desc)
  27         ## The variance of `node' is equal to the covariance of
  28         ## each possible pair among its descendant clades.
  29         for (i in 1:(nd - 1))
  30           for (j in (i + 1):nd)
  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
  34         for (i in 1:nd) {
  35             if (desc[i] <= n) next
  36             foo(desc[i], vcv[desc[i], desc[i]], to[i])
  37         }
  38     }
  39
  40     phy <- reorder(phy)
  41     n <- length(phy$tip.label)
  42     n.node <- phy$Nnode
  43     N <- n.node + n
  44     vcv <- matrix(0, N, N)
  45     foo(n + 1, 0, dim(phy$edge)[1])
  46
  47     vcv <- vcv[1:n, 1:n]
  48     if (corr) {
  49         ## This is inspired from the code of `cov2cor' (2005-09-08):
  50         M <- vcv
  51         Is <- sqrt(1/M[1 + 0:(n - 1)*(n + 1)])
  52         vcv[] <- Is * M * rep(Is, each = n)
  53         vcv[1 + 0:(n - 1)*(n + 1)] <- 1
  54     }
  55     rownames(vcv) <- colnames(vcv) <- phy$tip.label
  56     vcv
  57 }
  58
  59 vcv.corPhyl <- function(phy, corr = FALSE, ...)
  60 {
  61     labels <- attr(phy, "tree")$tip.label
  62     dummy.df <- data.frame(1:length(labels), row.names = labels)
  63     res <- nlme::corMatrix(Initialize.corPhyl(phy, dummy.df), corr = corr)
  64     dimnames(res) <- list(labels, labels)
  65     res
  66 }