R/rTrait.R

   1 ## rTrait.R (2010-01-11)
   2
   3 ##   Trait Evolution
   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 rTraitDisc <-
  11     function(phy, model = "ER", k = if (is.matrix(model)) ncol(model) else 2,
  12              rate = 0.1, states = LETTERS[1:k], freq = rep(1/k, k),
  13              ancestor = FALSE, root.value = 1)
  14 {
  15     if (is.null(phy$edge.length))
  16         stop("tree has no branch length")
  17     if (any(phy$edge.length < 0))
  18         stop("at least one branch length negative")
  19
  20     if (is.character(model)) {
  21         switch(model, "ER" = {
  22                    if (length(rate) != 1)
  23                        stop("`rate' must have one element")
  24                    Q <- matrix(rate, k, k)
  25                }, "ARD" = {
  26                    if (length(rate) != k*(k - 1))
  27                        stop("`rate' must have k(k - 1) elements")
  28                    Q <- matrix(0, k, k)
  29                    Q[col(Q) != row(Q)] <- rate
  30                }, "SYM" = {
  31                    if (length(rate) != k*(k - 1)/2)
  32                        stop("`rate' must have k(k - 1)/2 elements")
  33                    Q <- matrix(0, k, k)
  34                    sel <- col(Q) < row(Q)
  35                    Q[sel] <- rate
  36                    Q <- t(Q)
  37                    Q[sel] <- rate
  38                })
  39     }
  40     if (is.matrix(model)) {
  41         Q <- model
  42         if (ncol(Q) != nrow(Q))
  43             stop("the matrix given as `model' must be square")
  44     }
  45
  46     phy <- reorder(phy, "pruningwise")
  47     n <- length(phy$tip.label)
  48     N <- dim(phy$edge)[1]
  49     ROOT <- n + 1L
  50     x <- integer(n + phy$Nnode)
  51     x[ROOT] <- as.integer(root.value)
  52
  53     anc <- phy$edge[, 1]
  54     des <- phy$edge[, 2]
  55     el <- phy$edge.length
  56
  57     if (is.function(model)) {
  58         environment(model) <- environment() # to find 'k'
  59         for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
  60     } else {
  61         diag(Q) <- -rowSums(Q)
  62         for (i in N:1) {
  63             p <- matexpo(Q * freq * el[i])[x[anc[i]], ]
  64             x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
  65         }
  66     }
  67
  68     if (ancestor) {
  69         if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
  70         names(x) <- c(phy$tip.label, phy$node.label)
  71     } else {
  72         x <- x[1:n]
  73         names(x) <- phy$tip.label
  74     }
  75     class(x) <- "factor"
  76     levels(x) <- states
  77     x
  78 }
  79
  80 rTraitCont <-
  81     function(phy, model = "BM", sigma = 0.1, alpha = 1,
  82              theta = 0, ancestor = FALSE, root.value = 0)
  83 {
  84     if (is.null(phy$edge.length))
  85         stop("tree has no branch length")
  86     if (any(phy$edge.length < 0))
  87         stop("at least one branch length negative")
  88
  89     phy <- reorder(phy, "pruningwise")
  90     n <- length(phy$tip.label)
  91     N <- dim(phy$edge)[1]
  92     ROOT <- n + 1L
  93     x <- numeric(n + phy$Nnode)
  94     x[ROOT] <- root.value
  95
  96     anc <- phy$edge[, 1]
  97     des <- phy$edge[, 2]
  98     el <- phy$edge.length
  99
 100     if (is.function(model)) {
 101         environment(model) <- environment()
 102         for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
 103     } else {
 104         model <- pmatch(model, c("BM", "OU"))
 105         if (length(sigma) == 1) sigma <- rep(sigma, N)
 106         else if (length(sigma) != N)
 107             stop("'sigma' must have one or Nedge(phy) elements")
 108         if (model == 2) { # "OU"
 109             if (length(alpha) == 1) alpha <- rep(alpha, N)
 110             else if (length(alpha) != N)
 111                 stop("'alpha' must have one or Nedge(phy) elements")
 112             if (length(theta) == 1) theta <- rep(theta, N)
 113             else if (length(theta) != N)
 114                 stop("'theta' must have one or Nedge(phy) elements")
 115         }
 116         .C("rTraitCont", as.integer(model), as.integer(N), as.integer(anc - 1L), as.integer(des - 1L), el, sigma, alpha, theta, x, DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
 117     }
 118
 119     if (ancestor) {
 120         if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
 121         names(x) <- c(phy$tip.label, phy$node.label)
 122     } else {
 123         x <- x[1:n]
 124         names(x) <- phy$tip.label
 125     }
 126     x
 127 }