-## pic.R (2012-09-11)
+## pic.R (2013-02-18)
## Phylogenetically Independent Contrasts
-## Copyright 2002-2012 Emmanuel Paradis
+## Copyright 2002-2013 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE)
+pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE, rescaled.tree = FALSE)
{
- if (!inherits(phy, "phylo"))
- stop("object 'phy' is not of class \"phylo\"")
- if (is.null(phy$edge.length))
- stop("your tree has no branch lengths")
+ if (!inherits(phy, "phylo")) stop("object 'phy' is not of class \"phylo\"")
+ if (is.null(phy$edge.length)) stop("your tree has no branch lengths")
nb.tip <- length(phy$tip.label)
nb.node <- phy$Nnode
if (nb.node != nb.tip - 1)
- stop("'phy' is not rooted and fully dichotomous")
+ stop("'phy' is not rooted and fully dichotomous")
if (length(x) != nb.tip)
- stop("length of phenotypic and of phylogenetic data do not match")
+ stop("length of phenotypic and of phylogenetic data do not match")
if (any(is.na(x)))
- stop("missing data in 'x': you may consider removing the species with missing data from your tree with the function 'drop.tip'.")
+ stop("missing data in 'x': you may consider removing the species with missing data from your tree with the function 'drop.tip'.")
phy <- reorder(phy, "postorder")
phenotype <- numeric(nb.tip + nb.node)
warning("the names of argument 'x' and the tip labels of the tree did not match: the former were ignored in the analysis.")
}
}
+
## No need to copy the branch lengths: they are rescaled
## in the C code, so it's important to leave the default
## `DUP = TRUE' of .C.
-
ans <- .C("pic", as.integer(nb.tip), as.integer(nb.node),
as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
as.double(phy$edge.length), as.double(phenotype),
as.integer(var.contrasts), as.integer(scaled),
PACKAGE = "ape")
- ## The "old" R code:
- ##for (i in seq(from = 1, by = 2, length.out = nb.node)) {
- ## j <- i + 1
- ## anc <- phy$edge[i, 1]
- ## des1 <- phy$edge[i, 2]
- ## des2 <- phy$edge[j, 2]
- ## sumbl <- bl[i] + bl[j]
- ## ic <- anc - nb.tip
- ## contr[ic] <- phenotype[des1] - phenotype[des2]
- ## if (scaled) contr[ic] <- contr[ic]/sqrt(sumbl)
- ## if (var.contrasts) var.con[ic] <- sumbl
- ## phenotype[anc] <- (phenotype[des1]*bl[j] + phenotype[des2]*bl[i])/sumbl
- ## k <- which(phy$edge[, 2] == anc)
- ## bl[k] <- bl[k] + bl[i]*bl[j]/sumbl
- ##
- ##}
contr <- ans[[7]]
lbls <-
if (is.null(phy$node.label)) as.character(1:nb.node + nb.tip)
contr <- cbind(contr, ans[[8]])
dimnames(contr) <- list(lbls, c("contrasts", "variance"))
} else names(contr) <- lbls
+ if (rescaled.tree) {
+ phy$edge.length <- ans[[5]]
+ contr <- list(contr = contr, rescaled.tree = phy)
+ }
contr
}
{
n <- length(x)
m <- n - 1L # number of nodes
- phy <- reorder(phy, "pruningwise")
+ phy <- reorder(phy, "postorder")
xx <- unlist(lapply(x, mean)) # 'x' in Felsenstein's paper
xx <- c(xx, numeric(m))
delta.v <- numeric(n + m)