3 ## Phylogenetically Independent Contrasts
5 ## Copyright 2002-2009 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE)
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)
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")
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'.")
25 phy <- reorder(phy, "pruningwise")
26 phenotype <- numeric(nb.tip + nb.node)
28 if (is.null(names(x))) {
29 phenotype[1:nb.tip] <- x
31 if (all(names(x) %in% phy$tip.label))
32 phenotype[1:nb.tip] <- x[phy$tip.label]
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.')
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)
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),
51 ##for (i in seq(from = 1, by = 2, length.out = nb.node)) {
53 ## anc <- phy$edge[i, 1]
54 ## des1 <- phy$edge[i, 2]
55 ## des2 <- phy$edge[j, 2]
56 ## sumbl <- bl[i] + bl[j]
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
68 contr <- cbind(contr, ans[[8]])
69 dimnames(contr) <- list(1:nb.node + nb.tip, c("contrasts", "variance"))
70 } else names(contr) <- 1:nb.node + nb.tip