1 ## rTrait.R (2011-06-15)
5 ## Copyright 2010-2011 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
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, ...)
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")
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)
26 if (length(rate) != k*(k - 1))
27 stop("`rate' must have k(k - 1) elements")
29 Q[col(Q) != row(Q)] <- rate
31 if (length(rate) != k*(k - 1)/2)
32 stop("`rate' must have k(k - 1)/2 elements")
34 sel <- col(Q) < row(Q)
40 if (is.matrix(model)) {
42 if (ncol(Q) != nrow(Q))
43 stop("the matrix given as `model' must be square")
46 phy <- reorder(phy, "pruningwise")
47 n <- length(phy$tip.label)
50 x <- integer(n + phy$Nnode)
51 x[ROOT] <- as.integer(root.value)
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], ...)
61 freq <- rep(freq, each = k)
64 diag(Q) <- -rowSums(Q)
66 p <- matexpo(Q * el[i])[x[anc[i]], ]
67 x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
72 if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
73 names(x) <- c(phy$tip.label, phy$node.label)
76 names(x) <- phy$tip.label
84 function(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0,
85 ancestor = FALSE, root.value = 0, ...)
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")
92 phy <- reorder(phy, "pruningwise")
93 n <- length(phy$tip.label)
96 x <- numeric(n + phy$Nnode)
101 el <- phy$edge.length
103 if (is.function(model)) {
104 environment(model) <- environment()
105 for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i], ...)
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")
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")
126 if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
127 names(x) <- c(phy$tip.label, phy$node.label)
130 names(x) <- phy$tip.label
136 function(phy, model, p = 1, root.value = rep(0, p), ancestor = FALSE,
137 asFactor = NULL, trait.labels = paste("x", 1:p, sep = ""), ...)
139 phy <- reorder(phy, "pruningwise")
140 n <- length(phy$tip.label)
142 N <- dim(phy$edge)[1]
145 x <- matrix(0, n + m, p)
146 x[ROOT, ] <- root.value
151 el <- phy$edge.length
152 if (is.null(el)) el <- numeric(N)
154 environment(model) <- environment() # to find 'p'
156 for (i in N:1) x[des[i], ] <- model(x[anc[i], ], el[i], ...)
159 if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
160 rownames(x) <- c(phy$tip.label, phy$node.label)
162 x <- x[1:n, , drop = FALSE]
163 rownames(x) <- phy$tip.label
165 x <- as.data.frame(x)
166 names(x) <- trait.labels
167 if (!is.null(asFactor)) {
168 for (i in asFactor) {
170 x[, i] <- factor(y, labels = LETTERS[1:length(unique(y))])