Discrete Character Simulation
This function simulates the evolution of a discrete character along a
phylogeny. If model is a character or a matrix, evolution is
simulated with a Markovian model; the transition probabilities are
calculated for each branch with P = e^{Qt} where Q is the
rate matrix given by model and t is the branch length.
The calculation is done recursively from the root. See Paradis (2006,
p. 101) for a general introduction applied to evolution.
rTraitDisc(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, ...)phy |
an object of class |
model |
a character, a square numeric matrix, or a function specifying the model (see details). |
k |
the number of states of the character. |
rate |
the rate of change used if |
states |
the labels used for the states; by default “A”, “B”, ... |
freq |
a numeric vector giving the equilibrium relative frequencies of each state; by default the frequencies are equal. |
ancestor |
a logical value specifying whether to return the values at the nodes as well (by default, only the values at the tips are returned). |
root.value |
an integer giving the value at the root (by default,
it's the first state). To have a random value, use |
... |
further arguments passed to |
There are three possibilities to specify model:
A matrix:it must be a numeric square matrix; the diagonal is
always ignored. The arguments k and rate are ignored.
A character:these are the same short-cuts than in the function
ace: "ER" is an equal-rates model, "ARD"
is an all-rates-different model, and "SYM" is a symmetrical
model. Note that the argument rate must be of the appropriate
length, i.e., 1, k(k - 1), or k(k - 1)/2 for the three models,
respectively. The rate matrix Q is then filled column-wise.
A function:it must be of the form foo(x, l) where
x is the trait of the ancestor and l is the branch
length. It must return the value of the descendant as an integer.
A factor with names taken from the tip labels of phy. If
ancestor = TRUE, the node labels are used if present,
otherwise, “Node1”, “Node2”, etc.
Emmanuel Paradis
Paradis, E. (2006) Analyses of Phylogenetics and Evolution with R. New York: Springer.
data(bird.orders)
### the two followings are the same:
rTraitDisc(bird.orders)
rTraitDisc(bird.orders, model = matrix(c(0, 0.1, 0.1, 0), 2))
### two-state model with irreversibility:
rTraitDisc(bird.orders, model = matrix(c(0, 0, 0.1, 0), 2))
### simple two-state model:
tr <- rcoal(n <- 40, br = runif)
x <- rTraitDisc(tr, ancestor = TRUE)
plot(tr, show.tip.label = FALSE)
nodelabels(pch = 19, col = x[-(1:n)])
tiplabels(pch = 19, col = x[1:n])
### an imaginary model with stasis 0.5 time unit after a node, then
### random evolution:
foo <- function(x, l) {
if (l < 0.5) return(x)
sample(2, size = 1)
}
tr <- rcoal(20, br = runif)
x <- rTraitDisc(tr, foo, ancestor = TRUE)
plot(tr, show.tip.label = FALSE)
co <- c("blue", "yellow")
cot <- c("white", "black")
Y <- x[1:20]
A <- x[-(1:20)]
nodelabels(A, bg = co[A], col = cot[A])
tiplabels(Y, bg = co[Y], col = cot[Y])Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.