R/MPR.R

   1 ## MPR.R (2010-08-10)
   2
   3 ##   Most Parsimonious Reconstruction
   4
   5 ## Copyright 2010 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 MPR <- function(x, phy, outgroup)
  11 {
  12     if (is.rooted(phy))
  13         stop("the tree must be unrooted")
  14     if (!is.binary.tree(phy))
  15         stop("the tree must be fully dichotomous")
  16     if (length(outgroup) > 1L)
  17         stop("outgroup must be a single tip")
  18     if (is.character(outgroup))
  19         outgroup <- which(phy$tip.label == outgroup)
  20
  21     if (!is.null(names(x))) {
  22         if (all(names(x) %in% phy$tip.label))
  23             x <- x[phy$tip.label]
  24         else warning("the names of 'x' and the tip labels of the tree do not match: the former were ignored in the analysis.")
  25     }
  26
  27     n <- length(phy$tip.label)
  28     if (is.null(phy$node.label))
  29         phy$node.label <- n + 1:(phy$Nnode)
  30
  31     phy <- drop.tip(root(phy, outgroup), outgroup)
  32     n <- n - 1L
  33     m <- phy$Nnode
  34     phy <- reorder(phy, "pruningwise")
  35
  36     root.state <- x[outgroup]
  37     I <- as.integer(x[-outgroup])
  38     I[n + 1:m] <- NA
  39     I <- cbind(I, I) # interval map
  40
  41     med <- function(x) {
  42         i <- length(x)/2
  43         sort(x)[c(i, i + 1L)]
  44     }
  45
  46     ## 1st pass
  47     s <- seq(from = 1, by = 2, length.out = m)
  48     anc <- phy$edge[s, 1]
  49     des <- matrix(phy$edge[, 2], ncol = 2, byrow = TRUE)
  50     for (i in 1:m) I[anc[i], ] <- med(I[des[i, ], ])
  51
  52     ## 2nd pass
  53     out <- matrix(NA, m, 2)
  54     colnames(out) <- c("lower", "upper")
  55     ## do the most basal node before looping
  56     Iw <- as.vector(I[des[m, ], ]) # interval maps of the descendants
  57     out[anc[m] - n, ] <- range(med(c(root.state, root.state, Iw)))
  58     for (i in (m - 1):1) {
  59         j <- anc[i]
  60         Iw <- as.vector(I[des[i, ], ]) # interval maps of the descendants
  61         k <- which(phy$edge[, 2] == j) # find the ancestor
  62         tmp <- out[phy$edge[k, 1] - n, ]
  63         out[j - n, 1] <- min(med(c(tmp[1], tmp[1], Iw)))
  64         out[j - n, 2] <- max(med(c(tmp[2], tmp[2], Iw)))
  65     }
  66     rownames(out) <- phy$node.label
  67     out
  68 }