R/ltt.plot.R

   1 ## ltt.plot.R (2012-03-05)
   2
   3 ##    Lineages Through Time Plot
   4
   5 ## Copyright 2002-2012 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 ltt.plot.coords <- function(phy, backward = TRUE, tol = 1e-6)
  11 {
  12     if (is.ultrametric(phy, tol)) {
  13         if (is.binary.tree(phy)) {
  14             N <- numeric(phy$Nnode + 1)
  15             N[] <- 1
  16         } else {
  17             node.order <- tabulate(phy$edge[, 1])
  18             N <- node.order[-(1:length(phy$tip.label))] - 1
  19         }
  20         bt <- branching.times(phy)
  21         names(bt) <- NULL
  22         o <- order(bt, decreasing = TRUE)
  23         time <- c(-bt[o], 0)
  24         if (!is.binary.tree(phy)) N <- c(1, N[o])
  25     } else {
  26         if (!is.binary.tree(phy)) phy <- multi2di(phy)
  27         n <- Ntip(phy)
  28         m <- phy$Nnode
  29         ROOT <- n + 1L
  30         event <- time.event <- numeric(n + m)
  31
  32         time.event[ROOT] <- 0
  33         phy <- reorder(phy)
  34
  35         for (i in 1:nrow(phy$edge))
  36             time.event[phy$edge[i, 2]] <- time.event[phy$edge[i, 1]] + phy$edge.length[i]
  37
  38         present <- max(time.event)
  39         event[1:n] <- -1
  40         event[ROOT:(n + m)] <- 1
  41
  42         ## delete the events that are too close to present:
  43         past.event <- present - time.event > tol
  44         event <- event[past.event]
  45         time.event <- time.event[past.event]
  46
  47         ## reorder wrt time:
  48         o <- order(time.event)
  49         time.event <- time.event[o]
  50         event <- event[o]
  51
  52         time <- c(time.event - present, 0)
  53         N <- c(1, event)
  54     }
  55     N <- cumsum(N)
  56     if (!is.null(phy$root.edge)) {
  57         time <- c(time[1] - phy$root.edge, time)
  58         N <- c(1, N)
  59     }
  60     if (!backward) time <- time - time[1]
  61     cbind(time, N)
  62 }
  63
  64 ltt.plot <- function(phy, xlab = "Time", ylab = "N",
  65                      backward = TRUE, tol = 1e-6, ...)
  66 {
  67     if (!inherits(phy, "phylo"))
  68         stop("object \"phy\" is not of class \"phylo\"")
  69
  70     xy <- ltt.plot.coords(phy, backward, tol)
  71
  72     plot.default(xy, xlab = xlab, ylab = ylab, xaxs = "r",
  73                  yaxs = "r", type = "S", ...)
  74 }
  75
  76 ltt.lines <- function(phy, backward = TRUE, tol = 1e-6, ...)
  77 {
  78     xy <- ltt.plot.coords(phy, backward, tol)
  79     lines(xy, type = "S", ...)
  80 }
  81
  82 mltt.plot <-
  83     function(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE,
  84              xlab = "Time", ylab = "N", log = "", backward = TRUE, tol = 1e-6)
  85 {
  86     if (inherits(phy, "phylo")) { # if a tree of class "phylo"
  87         TREES <- list(ltt.plot.coords(phy, backward, tol))
  88         names(TREES) <- deparse(substitute(phy))
  89     } else { # a list of trees
  90         TREES <- lapply(phy, ltt.plot.coords, backward = backward, tol = tol)
  91         names(TREES) <- names(phy)
  92         if (is.null(names(TREES)))
  93             names(TREES) <-
  94                 paste(deparse(substitute(phy)), "-", 1:length(TREES))
  95     }
  96     dts <- list(...)
  97     n <- length(dts)
  98     if (n) {
  99         mc <- as.character(match.call())[-(1:2)]
 100         nms <- mc[1:n]
 101         for (i in 1:n) {
 102             if (inherits(dts[[i]], "phylo")) {
 103                 a <- list(ltt.plot.coords(dts[[i]], backward, tol))
 104                 names(a) <- nms[i]
 105             } else { # a list of trees
 106                 a <- lapply(dts[[i]], ltt.plot.coords, backward = backward, tol = tol)
 107                 names(a) <- names(dts[[i]])
 108                 if (is.null(names(a)))
 109                     names(a) <- paste(deparse(substitute(phy)), "-", seq_along(a))
 110             }
 111             TREES <- c(TREES, a)
 112         }
 113     }
 114     n <- length(TREES)
 115     range.each.tree <- sapply(TREES, function(x) range(x[, 1]))
 116     xl <- range(range.each.tree)
 117     yl <- c(1, max(sapply(TREES, function(x) max(x[, 2]))))
 118
 119     ## if backward is FALSE, we have to rescale the time scales of each tree:
 120     if (!backward) {
 121         for (i in seq_along(TREES)) {
 122             tmp <- TREES[[i]]
 123             tmp[, 1] <- tmp[, 1] + xl[2] - range.each.tree[2, i]
 124             TREES[[i]] <- tmp
 125         }
 126     }
 127
 128     plot.default(NA, type = "n", xlim = xl, ylim = yl, xaxs = "r",
 129                  yaxs = "r", xlab = xlab, ylab = ylab, log = log)
 130
 131     lty <- if (!dlty) rep(1, n) else 1:n
 132     col <- if (!dcol) rep(1, n) else topo.colors(n)
 133
 134     for (i in 1:n)
 135         lines(TREES[[i]], col = col[i], lty = lty[i], type = "S")
 136
 137     if (legend)
 138         legend(xl[1], yl[2], legend = names(TREES),
 139                lty = lty, col = col, bty = "n")
 140 }
 141
 142 ltt.coplot <- function(phy, backward = TRUE, ...)
 143 {
 144     layout(matrix(1:2, 2))
 145     par(mar = c(0, 3, 0.5, 0.5))
 146     o <- plot(phy, root.edge = TRUE, ...)
 147     par(mar = c(3, 3, 0, 0.5))
 148     ltt.plot(phy, xlim = o$x.lim, backward = FALSE, xaxt = "n")
 149     if (backward) axisPhylo() else axis(1)
 150 }