+
+pic.ortho <- function(x, phy, var.contrasts = FALSE, intra = FALSE)
+{
+ n <- length(x)
+ m <- n - 1L # number of nodes
+ phy <- reorder(phy, "pruningwise")
+ xx <- unlist(lapply(x, mean)) # 'x' in Felsenstein's paper
+ xx <- c(xx, numeric(m))
+ delta.v <- numeric(n + m)
+ s <- 1/unlist(lapply(x, length))
+ s <- c(s, numeric(m))
+ contrast <- var.cont <- numeric(m)
+
+ i <- 1L
+ while (i < m + n) {
+ d1 <- phy$edge[i, 2]
+ d2 <- phy$edge[i + 1L, 2]
+ a <- phy$edge[i, 1]
+ tmp1 <- 1/(phy$edge.length[i] + delta.v[d1])
+ tmp2 <- 1/(phy$edge.length[i + 1L] + delta.v[d2])
+ xx[a] <- (tmp1 * xx[d1] + tmp2 * xx[d2])/(tmp1 + tmp2)
+ delta.v[a] <- 1/(tmp1 + tmp2)
+ f1 <- tmp1/(tmp1 + tmp2)
+ f2 <- tmp2/(tmp1 + tmp2)
+ s[a] <- f1*f1 * s[d1] + f2*f2 * s[d2]
+ tmp <- 1/(s[d1] + s[d2])
+ contrast[a - n] <- (xx[d1] - xx[d2]) * sqrt(tmp)
+ var.cont[a - n] <- (1/tmp1 + 1/tmp2) * tmp
+ i <- i + 2L
+ }
+
+ lbls <-
+ if (is.null(phy$node.label)) as.character(1:m + n)
+ else phy$node.label
+
+ if (var.contrasts) {
+ contrast <- cbind(contrast, var.cont)
+ dimnames(contrast) <- list(lbls, c("contrasts", "variance"))
+ } else names(contrast) <- lbls
+
+ if (intra) {
+ intraspe.ctr <- function(x) {
+ k <- length(x) - 1L
+ if (!k) return(NULL)
+ ctr <- numeric(k)
+ ctr[1L] <- x[1L] - x[2L]
+ if (k > 1)
+ for (i in 2:k)
+ ctr[i] <- x[i + 1L] - mean(x[1:i])
+ sqrt((1:k)/(1:k + 1)) * ctr
+ }
+ tmp <- lapply(x, intraspe.ctr)
+ names(tmp) <- phy$tip.label
+ attr(contrast, "intra") <- tmp
+ }
+
+ contrast
+}
+
+varCompPhylip <- function(x, phy, exec = NULL)
+{
+ n <- Ntip(phy)
+ if (is.vector(x)) x <- as.list(x)
+ if (is.matrix(x) || is.data.frame(x)) {
+ tmpx <- vector("list", n)
+ for (i in 1:n) tmpx[[i]] <- x[i, , drop = FALSE]
+ names(tmpx) <- rownames(x)
+ x <- tmpx
+ }
+ p <- if (is.vector(x[[1]])) 1L else ncol(x[[1]])
+ if (!is.null(names(x))) x <- x[phy$tip.label]
+
+ phy <- makeLabel(phy, len = 10)
+ lbs <- phy$tip.label
+
+ ni <- sapply(x, function(xx) if (is.vector(xx)) 1L else nrow(xx))
+
+ pfx <- tempdir()
+ write.tree(phy, file = paste(pfx, "intree", sep = "/"))
+ infile <- paste(pfx, "infile", sep = "/")
+ file.create(infile)
+ cat(n, " ", p, "\n", sep = "", file = infile, append = TRUE)
+ for (i in 1:n) {
+ cat(lbs[i], file = infile, append = TRUE)
+ ## can surely be better but OK for the moment:
+ cat(paste(rep(" ", 11 - nchar(lbs[i])), collapse = ""),
+ file = infile, append = TRUE)
+ cat(ni[i], "\n", sep = "", file = infile, append = TRUE)
+ if (ni[i] == 1) {
+ cat(x[[i]], sep = " ", file = infile, append = TRUE)
+ cat("\n", file = infile, append = TRUE)
+ } else write(t(x[[i]]), file = infile, ncolumns = p, append = TRUE)
+ }
+
+ if (is.null(exec))
+ exec <-
+ if (.Platform$OS.type == "unix") "phylip contrast"
+ else "contrast"
+
+ odir <- setwd(pfx)
+ on.exit(setwd(odir))
+ if (file.exists("outfile")) unlink("outfile")
+ system("phylip contrast", intern = TRUE, input = c("W", "A", "Y"))
+ varA <- scan("outfile", skip = 7, nlines = p, quiet = TRUE)
+ varE <- scan("outfile", skip = 11 + p, nlines = p, quiet = TRUE)
+ if (p > 1) {
+ varA <- matrix(varA, p, p, byrow = TRUE)
+ varE <- matrix(varE, p, p, byrow = TRUE)
+ }
+ list(varA = varA, varE = varE)
+}