R/pic.R

   1 ## pic.R (2012-11-20)
   2
   3 ##   Phylogenetically Independent Contrasts
   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 pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE, rescaled.tree = 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")
  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("missing data in 'x': you may consider removing the species with missing data from your tree with the function 'drop.tip'.")
  24
  25     phy <- reorder(phy, "postorder")
  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
  39     ## No need to copy the branch lengths: they are rescaled
  40     ## in the C code, so it's important to leave the default
  41     ## `DUP = TRUE' of .C.
  42     ans <- .C("pic", as.integer(nb.tip), as.integer(nb.node),
  43               as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
  44               as.double(phy$edge.length), as.double(phenotype),
  45               double(nb.node), double(nb.node),
  46               as.integer(var.contrasts), as.integer(scaled),
  47               PACKAGE = "ape")
  48
  49     contr <- ans[[7]]
  50     lbls <-
  51         if (is.null(phy$node.label)) as.character(1:nb.node + nb.tip)
  52         else phy$node.label
  53     if (var.contrasts) {
  54         contr <- cbind(contr, ans[[8]])
  55         dimnames(contr) <- list(lbls, c("contrasts", "variance"))
  56     } else names(contr) <- lbls
  57     if (rescaled.tree) {
  58         phy$edge.length <- ans[[5]]
  59         contr <- list(contr = contr, rescaled.tree = phy)
  60     }
  61     contr
  62 }
  63
  64 pic.ortho <- function(x, phy, var.contrasts = FALSE, intra = FALSE)
  65 {
  66     n <- length(x)
  67     m <- n - 1L # number of nodes
  68     phy <- reorder(phy, "pruningwise")
  69     xx <- unlist(lapply(x, mean)) # 'x' in Felsenstein's paper
  70     xx <- c(xx, numeric(m))
  71     delta.v <- numeric(n + m)
  72     s <- 1/unlist(lapply(x, length))
  73     s <- c(s, numeric(m))
  74     contrast <- var.cont <- numeric(m)
  75
  76     i <- 1L
  77     while (i < m + n) {
  78         d1 <- phy$edge[i, 2]
  79         d2 <- phy$edge[i + 1L, 2]
  80         a <- phy$edge[i, 1]
  81         tmp1 <- 1/(phy$edge.length[i] + delta.v[d1])
  82         tmp2 <- 1/(phy$edge.length[i + 1L] + delta.v[d2])
  83         xx[a] <- (tmp1 * xx[d1] + tmp2 * xx[d2])/(tmp1 + tmp2)
  84         delta.v[a] <- 1/(tmp1 + tmp2)
  85         f1 <- tmp1/(tmp1 + tmp2)
  86         f2 <- tmp2/(tmp1 + tmp2)
  87         s[a] <- f1*f1 * s[d1] + f2*f2 * s[d2]
  88         tmp <- 1/(s[d1] + s[d2])
  89         contrast[a - n] <- (xx[d1] - xx[d2]) * sqrt(tmp)
  90         var.cont[a - n] <- (1/tmp1 + 1/tmp2) * tmp
  91         i <- i + 2L
  92     }
  93
  94     lbls <-
  95         if (is.null(phy$node.label)) as.character(1:m + n)
  96         else phy$node.label
  97
  98     if (var.contrasts) {
  99         contrast <- cbind(contrast, var.cont)
 100         dimnames(contrast) <- list(lbls, c("contrasts", "variance"))
 101     } else names(contrast) <- lbls
 102
 103     if (intra) {
 104         intraspe.ctr <- function(x) {
 105             k <- length(x) - 1L
 106             if (!k) return(NULL)
 107             ctr <- numeric(k)
 108             ctr[1L] <- x[1L] - x[2L]
 109             if (k > 1)
 110                 for (i in 2:k)
 111                     ctr[i] <- x[i + 1L] - mean(x[1:i])
 112             sqrt((1:k)/(1:k + 1)) * ctr
 113         }
 114         tmp <- lapply(x, intraspe.ctr)
 115         names(tmp) <- phy$tip.label
 116         attr(contrast, "intra") <- tmp
 117     }
 118
 119     contrast
 120 }
 121
 122 varCompPhylip <- function(x, phy, exec = NULL)
 123 {
 124     n <- Ntip(phy)
 125     if (is.vector(x)) x <- as.list(x)
 126     if (is.matrix(x) || is.data.frame(x)) {
 127         tmpx <- vector("list", n)
 128         for (i in 1:n) tmpx[[i]] <- x[i, , drop = FALSE]
 129         names(tmpx) <- rownames(x)
 130         x <- tmpx
 131     }
 132     p <- if (is.vector(x[[1]])) 1L else ncol(x[[1]])
 133     if (!is.null(names(x))) x <- x[phy$tip.label]
 134
 135     phy <- makeLabel(phy, len = 10)
 136     lbs <- phy$tip.label
 137
 138     ni <- sapply(x, function(xx) if (is.vector(xx)) 1L else nrow(xx))
 139
 140     pfx <- tempdir()
 141     write.tree(phy, file = paste(pfx, "intree", sep = "/"))
 142     infile <- paste(pfx, "infile", sep = "/")
 143     file.create(infile)
 144     cat(n, " ", p, "\n", sep = "", file = infile, append = TRUE)
 145     for (i in 1:n) {
 146         cat(lbs[i], file = infile, append = TRUE)
 147         ## can surely be better but OK for the moment:
 148         cat(paste(rep(" ", 11 - nchar(lbs[i])), collapse = ""),
 149             file = infile, append = TRUE)
 150         cat(ni[i], "\n", sep = "", file = infile, append = TRUE)
 151         if (ni[i] == 1) {
 152             cat(x[[i]], sep = " ", file = infile, append = TRUE)
 153             cat("\n", file = infile, append = TRUE)
 154         } else write(t(x[[i]]), file = infile, ncolumns = p, append = TRUE)
 155     }
 156
 157     if (is.null(exec))
 158         exec <-
 159             if (.Platform$OS.type == "unix") "phylip contrast"
 160             else "contrast"
 161
 162     odir <- setwd(pfx)
 163     on.exit(setwd(odir))
 164     if (file.exists("outfile")) unlink("outfile")
 165     system(exec, intern = TRUE, input = c("W", "A", "Y"))
 166     varA <- scan("outfile", skip = 7, nlines = p, quiet = TRUE)
 167     varE <- scan("outfile", skip = 11 + p, nlines = p, quiet = TRUE)
 168     if (p > 1) {
 169         varA <- matrix(varA, p, p, byrow = TRUE)
 170         varE <- matrix(varE, p, p, byrow = TRUE)
 171     }
 172     list(varA = varA, varE = varE)
 173 }