3 ## Most Parsimonious Reconstruction
5 ## Copyright 2010 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 MPR <- function(x, phy, outgroup)
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)
21 if (!is.null(names(x))) {
22 if (all(names(x) %in% 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.")
27 n <- length(phy$tip.label)
28 if (is.null(phy$node.label))
29 phy$node.label <- n + 1:(phy$Nnode)
31 phy <- drop.tip(root(phy, outgroup), outgroup)
34 phy <- reorder(phy, "pruningwise")
36 root.state <- x[outgroup]
37 I <- as.integer(x[-outgroup])
39 I <- cbind(I, I) # interval map
47 s <- seq(from = 1, by = 2, length.out = m)
49 des <- matrix(phy$edge[, 2], ncol = 2, byrow = TRUE)
50 for (i in 1:m) I[anc[i], ] <- med(I[des[i, ], ])
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) {
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)))
66 rownames(out) <- phy$node.label