R/pic.R

   1 ## pic.R (2010-08-18)
   2
   3 ##   Phylogenetically Independent Contrasts
   4
   5 ## Copyright 2002-2010 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE)
  11 {
  12     if (!inherits(phy, "phylo"))
  13       stop("object 'phy' is not of class \"phylo\"")
  14     if (is.null(phy$edge.length))
  15       stop("your tree has no branch lengths: you may consider setting them equal to one, or using the function `compute.brlen'.")
  16     nb.tip <- length(phy$tip.label)
  17     nb.node <- phy$Nnode
  18     if (nb.node != nb.tip - 1)
  19       stop("'phy' is not rooted and fully dichotomous")
  20     if (length(x) != nb.tip)
  21       stop("length of phenotypic and of phylogenetic data do not match")
  22     if (any(is.na(x)))
  23       stop("the present method cannot (yet) be used directly with missing data: you may consider removing the species with missing data from your tree with the function `drop.tip'.")
  24
  25     phy <- reorder(phy, "pruningwise")
  26     phenotype <- numeric(nb.tip + nb.node)
  27
  28     if (is.null(names(x))) {
  29         phenotype[1:nb.tip] <- x
  30     } else {
  31         if (all(names(x) %in% phy$tip.label))
  32           phenotype[1:nb.tip] <- x[phy$tip.label]
  33         else {
  34             phenotype[1:nb.tip] <- x
  35             warning('the names of argument "x" and the tip labels of the tree did not match: the former were ignored in the analysis.')
  36         }
  37     }
  38     ## No need to copy the branch lengths: they are rescaled
  39     ## in the C code, so it's important to leave the default
  40     ## `DUP = TRUE' of .C.
  41     contr <- var.con <- numeric(nb.node)
  42
  43     ans <- .C("pic", as.integer(nb.tip), as.integer(nb.node),
  44               as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
  45               as.double(phy$edge.length), as.double(phenotype),
  46               as.double(contr), as.double(var.con),
  47               as.integer(var.contrasts), as.integer(scaled),
  48               PACKAGE = "ape")
  49
  50     ## The "old" R code:
  51     ##for (i in seq(from = 1, by = 2, length.out = nb.node)) {
  52     ##    j <- i + 1
  53     ##    anc <- phy$edge[i, 1]
  54     ##    des1 <- phy$edge[i, 2]
  55     ##    des2 <- phy$edge[j, 2]
  56     ##    sumbl <- bl[i] + bl[j]
  57     ##    ic <- anc - nb.tip
  58     ##    contr[ic] <- phenotype[des1] - phenotype[des2]
  59     ##    if (scaled) contr[ic] <- contr[ic]/sqrt(sumbl)
  60     ##    if (var.contrasts) var.con[ic] <- sumbl
  61     ##    phenotype[anc] <- (phenotype[des1]*bl[j] + phenotype[des2]*bl[i])/sumbl
  62     ##    k <- which(phy$edge[, 2] == anc)
  63     ##    bl[k] <- bl[k] + bl[i]*bl[j]/sumbl
  64     ##
  65     ##}
  66     contr <- ans[[7]]
  67     lbls <-
  68         if (is.null(phy$node.label)) as.character(1:nb.node + nb.tip)
  69         else phy$node.label
  70     if (var.contrasts) {
  71         contr <- cbind(contr, ans[[8]])
  72         dimnames(contr) <- list(lbls, c("contrasts", "variance"))
  73     } else names(contr) <- lbls
  74     contr
  75 }