R/rTrait.R

   1 ## rTrait.R (2011-06-15)
   2
   3 ##   Trait Evolution
   4
   5 ## Copyright 2010-2011 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(toupper(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         freq <- rep(freq, each = k)
  62         Q <- Q * freq
  63         diag(Q) <- 0
  64         diag(Q) <- -rowSums(Q)
  65         for (i in N:1) {
  66             p <- matexpo(Q * el[i])[x[anc[i]], ]
  67             x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
  68         }
  69     }
  70
  71     if (ancestor) {
  72         if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
  73         names(x) <- c(phy$tip.label, phy$node.label)
  74     } else {
  75         x <- x[1:n]
  76         names(x) <- phy$tip.label
  77     }
  78     class(x) <- "factor"
  79     levels(x) <- states
  80     x
  81 }
  82
  83 rTraitCont <-
  84     function(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0,
  85              ancestor = FALSE, root.value = 0, ...)
  86 {
  87     if (is.null(phy$edge.length))
  88         stop("tree has no branch length")
  89     if (any(phy$edge.length < 0))
  90         stop("at least one branch length negative")
  91
  92     phy <- reorder(phy, "pruningwise")
  93     n <- length(phy$tip.label)
  94     N <- dim(phy$edge)[1]
  95     ROOT <- n + 1L
  96     x <- numeric(n + phy$Nnode)
  97     x[ROOT] <- root.value
  98
  99     anc <- phy$edge[, 1]
 100     des <- phy$edge[, 2]
 101     el <- phy$edge.length
 102
 103     if (is.function(model)) {
 104         environment(model) <- environment()
 105         for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i], ...)
 106     } else {
 107         model <- pmatch(toupper(model), c("BM", "OU"))
 108         if (length(sigma) == 1) sigma <- rep(sigma, N)
 109         else if (length(sigma) != N)
 110             stop("'sigma' must have one or Nedge(phy) elements")
 111         if (model == 2) { # "OU"
 112             if (length(alpha) == 1) alpha <- rep(alpha, N)
 113             else if (length(alpha) != N)
 114                 stop("'alpha' must have one or Nedge(phy) elements")
 115             if (length(theta) == 1) theta <- rep(theta, N)
 116             else if (length(theta) != N)
 117                 stop("'theta' must have one or Nedge(phy) elements")
 118         }
 119         .C("rTraitCont", as.integer(model), as.integer(N),
 120            as.integer(anc - 1L), as.integer(des - 1L), el,
 121            as.double(sigma), as.double(alpha), as.double(theta), x,
 122            DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
 123     }
 124
 125     if (ancestor) {
 126         if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
 127         names(x) <- c(phy$tip.label, phy$node.label)
 128     } else {
 129         x <- x[1:n]
 130         names(x) <- phy$tip.label
 131     }
 132     x
 133 }
 134
 135 rTraitMult <-
 136     function(phy, model, p = 1, root.value = rep(0, p), ancestor = FALSE,
 137              asFactor = NULL, trait.labels = paste("x", 1:p, sep = ""), ...)
 138 {
 139     phy <- reorder(phy, "pruningwise")
 140     n <- length(phy$tip.label)
 141     m <- phy$Nnode
 142     N <- dim(phy$edge)[1]
 143     ROOT <- n + 1L
 144
 145     x <- matrix(0, n + m, p)
 146     x[ROOT, ] <- root.value
 147
 148     anc <- phy$edge[, 1]
 149     des <- phy$edge[, 2]
 150
 151     el <- phy$edge.length
 152     if (is.null(el)) el <- numeric(N)
 153
 154     environment(model) <- environment() # to find 'p'
 155
 156     for (i in N:1) x[des[i], ] <- model(x[anc[i], ], el[i], ...)
 157
 158     if (ancestor) {
 159         if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
 160         rownames(x) <- c(phy$tip.label, phy$node.label)
 161     } else {
 162         x <- x[1:n, , drop = FALSE]
 163         rownames(x) <- phy$tip.label
 164     }
 165     x <- as.data.frame(x)
 166     names(x) <- trait.labels
 167     if (!is.null(asFactor)) {
 168         for (i in asFactor) {
 169             y <- x[, i]
 170             x[, i] <- factor(y, labels = LETTERS[1:length(unique(y))])
 171         }
 172     }
 173     x
 174 }