1 ## chronoMPL.R (2007-08-29)
3 ## Molecular Dating with Mean Path Lengths
5 ## Copyright 2007 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 chronoMPL <- function(phy, se = TRUE, test = TRUE)
12 if (!is.binary.tree(phy)) stop("the tree is not dichotomous.")
13 n <- length(phy$tip.label)
16 obj <- reorder(phy, "pruningwise")
17 ndesc <- .C("node_depth", as.integer(n), as.integer(m),
18 as.integer(obj$edge[, 1]), as.integer(obj$edge[, 2]),
19 as.integer(N), double(n + m), DUP = FALSE,
21 s <- numeric(n + m) # sum of path lengths
23 if (test) Pval <- numeric(m)
24 for (i in seq(1, N - 1, 2)) {
29 A <- s[a] + ndesc[a]*obj$edge.length[i]
30 B <- s[b] + ndesc[b]*obj$edge.length[j]
33 ss[o] <- ss[a] + ndesc[a]^2 * obj$edge.length[i] + ss[b] +
34 ndesc[b]^2 * obj$edge.length[j]
36 z <- abs(A/ndesc[a] - B/ndesc[b])
37 tmp <- (ss[a] + ndesc[a]^2 * obj$edge.length[i])/ndesc[a]^2
38 tmp <- tmp + (ss[b] + ndesc[b]^2 * obj$edge.length[j])/ndesc[b]^2
40 Pval[o - n] <- 2*pnorm(z, lower.tail = FALSE)
44 phy$edge.length <- node.age[phy$edge[, 1]] - node.age[phy$edge[, 2]]
45 if (se) attr(phy, "stderr") <- sqrt(ss[-(1:n)]/ndesc[-(1:n)]^2)
46 if (test) attr(phy, "Pval") <- Pval