rTraitDisc {ape}R Documentation

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, ...)



an object of class "phylo".


a character, a square numeric matrix, or a function specifying the model (see details).


the number of states of the character.


the rate of change used if model is a character; it is not recycled if model = "ARD" of model = "SYM".


the labels used for the states; by default “A”, “B”, ...


a numeric vector giving the equilibrium relative frequencies of each state; by default the frequencies are equal.


a logical value specifying whether to return the values at the nodes as well (by default, only the values at the tips are returned).


an integer giving the value at the root (by default, it's the first state). To have a random value, use root.value = sample(k).


further arguments passed to model if it is a function.


There are three possibilities to specify model:


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.

See Also

rTraitCont, rTraitMult, ace


### the two followings are the same:
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])

[Package ape version 5.5 Index]