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


Paradis, E. (2006) Analyses of Phylogenetics and Evolution with R. New York: Springer.

See Also

rTraitCont, rTraitMult, ace

  • rTraitDisc
### 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])
# }
Documentation reproduced from package ape, version 5.3, License: GPL (>= 2)

Community examples

Looks like there are no examples yet.