R/chronoMPL.R

   1 ## chronoMPL.R (2007-08-29)
   2
   3 ##   Molecular Dating with Mean Path Lengths
   4
   5 ## Copyright 2007 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 chronoMPL <- function(phy, se = TRUE, test = TRUE)
  11 {
  12     if (!is.binary.tree(phy)) stop("the tree is not dichotomous.")
  13     n <- length(phy$tip.label)
  14     m <- phy$Nnode
  15     N <- dim(phy$edge)[1]
  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,
  20                 PACKAGE = "ape")[[6]]
  21     s <- numeric(n + m) # sum of path lengths
  22     if (se) ss <- s
  23     if (test) Pval <- numeric(m)
  24     for (i in seq(1, N - 1, 2)) {
  25         j <- i + 1
  26         a <- obj$edge[i, 2]
  27         b <- obj$edge[j, 2]
  28         o <- obj$edge[i, 1]
  29         A <- s[a] + ndesc[a]*obj$edge.length[i]
  30         B <- s[b] + ndesc[b]*obj$edge.length[j]
  31         s[o] <- A + B
  32         if (se)
  33           ss[o] <- ss[a] + ndesc[a]^2 * obj$edge.length[i] + ss[b] +
  34             ndesc[b]^2 * obj$edge.length[j]
  35         if (test) {
  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
  39             z <- z/sqrt(tmp)
  40             Pval[o - n] <- 2*pnorm(z, lower.tail = FALSE)
  41         }
  42     }
  43     node.age <- s/ndesc
  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
  47     phy
  48 }