-## rTrait.R (2010-01-11)
+## rTrait.R (2011-06-15)
## Trait Evolution
-## Copyright 2010 Emmanuel Paradis
+## Copyright 2010-2011 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
rTraitDisc <-
function(phy, model = "ER", k = if (is.matrix(model)) ncol(model) else 2,
rate = 0.1, states = LETTERS[1:k], freq = rep(1/k, k),
- ancestor = FALSE, root.value = 1)
+ ancestor = FALSE, root.value = 1, ...)
{
if (is.null(phy$edge.length))
stop("tree has no branch length")
stop("at least one branch length negative")
if (is.character(model)) {
- switch(model, "ER" = {
+ switch(toupper(model), "ER" = {
if (length(rate) != 1)
stop("`rate' must have one element")
Q <- matrix(rate, k, k)
if (is.function(model)) {
environment(model) <- environment() # to find 'k'
- for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
+ for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i], ...)
} else {
+ freq <- rep(freq, each = k)
+ Q <- Q * freq
+ diag(Q) <- 0
diag(Q) <- -rowSums(Q)
for (i in N:1) {
- p <- matexpo(Q * freq * el[i])[x[anc[i]], ]
+ p <- matexpo(Q * el[i])[x[anc[i]], ]
x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
}
}
}
rTraitCont <-
- function(phy, model = "BM", sigma = 0.1, alpha = 1,
- theta = 0, ancestor = FALSE, root.value = 0)
+ function(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0,
+ ancestor = FALSE, root.value = 0, ...)
{
if (is.null(phy$edge.length))
stop("tree has no branch length")
if (is.function(model)) {
environment(model) <- environment()
- for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
+ for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i], ...)
} else {
- model <- pmatch(model, c("BM", "OU"))
+ model <- pmatch(toupper(model), c("BM", "OU"))
if (length(sigma) == 1) sigma <- rep(sigma, N)
else if (length(sigma) != N)
stop("'sigma' must have one or Nedge(phy) elements")
else if (length(theta) != N)
stop("'theta' must have one or Nedge(phy) elements")
}
- .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")
+ .C("rTraitCont", as.integer(model), as.integer(N),
+ as.integer(anc - 1L), as.integer(des - 1L), el,
+ as.double(sigma), as.double(alpha), as.double(theta), x,
+ DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
}
if (ancestor) {
}
x
}
+
+rTraitMult <-
+ function(phy, model, p = 1, root.value = rep(0, p), ancestor = FALSE,
+ asFactor = NULL, trait.labels = paste("x", 1:p, sep = ""), ...)
+{
+ phy <- reorder(phy, "pruningwise")
+ n <- length(phy$tip.label)
+ m <- phy$Nnode
+ N <- dim(phy$edge)[1]
+ ROOT <- n + 1L
+
+ x <- matrix(0, n + m, p)
+ x[ROOT, ] <- root.value
+
+ anc <- phy$edge[, 1]
+ des <- phy$edge[, 2]
+
+ el <- phy$edge.length
+ if (is.null(el)) el <- numeric(N)
+
+ environment(model) <- environment() # to find 'p'
+
+ for (i in N:1) x[des[i], ] <- model(x[anc[i], ], el[i], ...)
+
+ if (ancestor) {
+ if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
+ rownames(x) <- c(phy$tip.label, phy$node.label)
+ } else {
+ x <- x[1:n, , drop = FALSE]
+ rownames(x) <- phy$tip.label
+ }
+ x <- as.data.frame(x)
+ names(x) <- trait.labels
+ if (!is.null(asFactor)) {
+ for (i in asFactor) {
+ y <- x[, i]
+ x[, i] <- factor(y, labels = LETTERS[1:length(unique(y))])
+ }
+ }
+ x
+}