R/vcv.phylo.R

   1 ## vcv.phylo.R (2012-02-21)
   2
   3 ##   Phylogenetic Variance-Covariance or Correlation Matrix
   4
   5 ## Copyright 2002-2012 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     pp <- prop.part(phy)
  18
  19     phy <- reorder(phy, "pruningwise")
  20
  21     n <- length(phy$tip.label)
  22     e1 <- phy$edge[, 1]
  23     e2 <- phy$edge[, 2]
  24     EL <- phy$edge.length
  25
  26     ## xx: vecteur donnant la distance d'un noeud
  27     ##     ou d'un tip à partir de la racine
  28     ## (same than in is.ultrametric)
  29     xx <- numeric(n + phy$Nnode)
  30
  31     vcv <- matrix(0, n, n)
  32
  33     ## the loop below starts from the bottom of the edge matrix, so
  34     ## from the root
  35
  36     for (i in length(e1):1) {
  37         var.cur.node <- xx[e1[i]]
  38         xx[e2[i]] <- var.cur.node + EL[i] # sets the variance
  39         j <- i - 1L
  40         while (e1[j] == e1[i] && j > 0) {
  41             left <- if (e2[j] > n) pp[[e2[j] - n]] else e2[j]
  42             right <- if (e2[i] > n) pp[[e2[i] - n]] else e2[i]
  43             vcv[left, right] <- vcv[right, left] <- var.cur.node # sets the covariance
  44             j <- j - 1L
  45         }
  46     }
  47
  48     diag.elts <- 1 + 0:(n - 1)*(n + 1)
  49     vcv[diag.elts] <- xx[1:n]
  50
  51     if (corr) {
  52         ## This is inspired from the code of cov2cor (2005-09-08):
  53         Is <- sqrt(1 / vcv[diag.elts])
  54         ## below 'vcv[] <- ...' has been changed to 'vcv <- ...'
  55         ## which seems to be twice faster for n = 1000 and
  56         ## respects the additional attributes (2012-02-21):
  57         vcv <- Is * vcv * rep(Is, each = n)
  58         vcv[diag.elts] <- 1
  59     }
  60
  61     dimnames(vcv)[1:2] <- list(phy$tip.label)
  62     vcv
  63 }
  64
  65 vcv.corPhyl <- function(phy, corr = FALSE, ...)
  66 {
  67     labels <- attr(phy, "tree")$tip.label
  68     dummy.df <- data.frame(1:length(labels), row.names = labels)
  69     res <- nlme::corMatrix(Initialize.corPhyl(phy, dummy.df), corr = corr)
  70     dimnames(res) <- list(labels, labels)
  71     res
  72 }