ape (version 4.0)

rTraitDisc: Discrete Character Simulation

Description

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.

Usage

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

Arguments

phy

an object of class "phylo".

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 model is a character; it is not recycled if model = "ARD" of model = "SYM".

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 root.value = sample(k).

...

further arguments passed to model if it is a function.

Value

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.

Details

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.

References

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

See Also

rTraitCont, rTraitMult, ace

Examples

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])