3 ## Phylogenetically Independent Contrasts
5 ## Copyright 2002-2011 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")
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("missing data in 'x': 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.
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),
50 ##for (i in seq(from = 1, by = 2, length.out = nb.node)) {
52 ## anc <- phy$edge[i, 1]
53 ## des1 <- phy$edge[i, 2]
54 ## des2 <- phy$edge[j, 2]
55 ## sumbl <- bl[i] + bl[j]
57 ## contr[ic] <- phenotype[des1] - phenotype[des2]
58 ## if (scaled) contr[ic] <- contr[ic]/sqrt(sumbl)
59 ## if (var.contrasts) var.con[ic] <- sumbl
60 ## phenotype[anc] <- (phenotype[des1]*bl[j] + phenotype[des2]*bl[i])/sumbl
61 ## k <- which(phy$edge[, 2] == anc)
62 ## bl[k] <- bl[k] + bl[i]*bl[j]/sumbl
67 if (is.null(phy$node.label)) as.character(1:nb.node + nb.tip)
70 contr <- cbind(contr, ans[[8]])
71 dimnames(contr) <- list(lbls, c("contrasts", "variance"))
72 } else names(contr) <- lbls
76 pic.ortho <- function(x, phy, var.contrasts = FALSE, intra = FALSE)
79 m <- n - 1L # number of nodes
80 phy <- reorder(phy, "pruningwise")
81 xx <- unlist(lapply(x, mean)) # 'x' in Felsenstein's paper
82 xx <- c(xx, numeric(m))
83 delta.v <- numeric(n + m)
84 s <- 1/unlist(lapply(x, length))
86 contrast <- var.cont <- numeric(m)
91 d2 <- phy$edge[i + 1L, 2]
93 tmp1 <- 1/(phy$edge.length[i] + delta.v[d1])
94 tmp2 <- 1/(phy$edge.length[i + 1L] + delta.v[d2])
95 xx[a] <- (tmp1 * xx[d1] + tmp2 * xx[d2])/(tmp1 + tmp2)
96 delta.v[a] <- 1/(tmp1 + tmp2)
97 f1 <- tmp1/(tmp1 + tmp2)
98 f2 <- tmp2/(tmp1 + tmp2)
99 s[a] <- f1*f1 * s[d1] + f2*f2 * s[d2]
100 tmp <- 1/(s[d1] + s[d2])
101 contrast[a - n] <- (xx[d1] - xx[d2]) * sqrt(tmp)
102 var.cont[a - n] <- (1/tmp1 + 1/tmp2) * tmp
107 if (is.null(phy$node.label)) as.character(1:m + n)
111 contrast <- cbind(contrast, var.cont)
112 dimnames(contrast) <- list(lbls, c("contrasts", "variance"))
113 } else names(contrast) <- lbls
116 intraspe.ctr <- function(x) {
120 ctr[1L] <- x[1L] - x[2L]
123 ctr[i] <- x[i + 1L] - mean(x[1:i])
124 sqrt((1:k)/(1:k + 1)) * ctr
126 tmp <- lapply(x, intraspe.ctr)
127 names(tmp) <- phy$tip.label
128 attr(contrast, "intra") <- tmp
134 varCompPhylip <- function(x, phy, exec = NULL)
137 if (is.vector(x)) x <- as.list(x)
138 if (is.matrix(x) || is.data.frame(x)) {
139 tmpx <- vector("list", n)
140 for (i in 1:n) tmpx[[i]] <- x[i, , drop = FALSE]
141 names(tmpx) <- rownames(x)
144 p <- if (is.vector(x[[1]])) 1L else ncol(x[[1]])
145 if (!is.null(names(x))) x <- x[phy$tip.label]
147 phy <- makeLabel(phy, len = 10)
150 ni <- sapply(x, function(xx) if (is.vector(xx)) 1L else nrow(xx))
153 write.tree(phy, file = paste(pfx, "intree", sep = "/"))
154 infile <- paste(pfx, "infile", sep = "/")
156 cat(n, " ", p, "\n", sep = "", file = infile, append = TRUE)
158 cat(lbs[i], file = infile, append = TRUE)
159 ## can surely be better but OK for the moment:
160 cat(paste(rep(" ", 11 - nchar(lbs[i])), collapse = ""),
161 file = infile, append = TRUE)
162 cat(ni[i], "\n", sep = "", file = infile, append = TRUE)
164 cat(x[[i]], sep = " ", file = infile, append = TRUE)
165 cat("\n", file = infile, append = TRUE)
166 } else write(t(x[[i]]), file = infile, ncolumns = p, append = TRUE)
171 if (.Platform$OS.type == "unix") "phylip contrast"
176 if (file.exists("outfile")) unlink("outfile")
177 system(exec, intern = TRUE, input = c("W", "A", "Y"))
178 varA <- scan("outfile", skip = 7, nlines = p, quiet = TRUE)
179 varE <- scan("outfile", skip = 11 + p, nlines = p, quiet = TRUE)
181 varA <- matrix(varA, p, p, byrow = TRUE)
182 varE <- matrix(varE, p, p, byrow = TRUE)
184 list(varA = varA, varE = varE)