ape (version 5.1)

rTraitCont: Continuous Character Simulation

Description

This function simulates the evolution of a continuous character along a phylogeny. The calculation is done recursively from the root. See Paradis (2012, pp. 232 and 324) for an introduction.

Usage

rTraitCont(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0,
           ancestor = FALSE, root.value = 0, ...)

Arguments

phy

an object of class "phylo".

model

a character (either "BM" or "OU") or a function specifying the model (see details).

sigma

a numeric vector giving the standard-deviation of the random component for each branch (can be a single value).

alpha

if model = "OU", a numeric vector giving the strength of the selective constraint for each branch (can be a single value).

theta

if model = "OU", a numeric vector giving the optimum for each branch (can be a single value).

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

a numeric giving the value at the root.

…

further arguments passed to model if it is a function.

Value

A numeric vector 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:

• "BM":a Browian motion model is used. If the arguments sigma has more than one value, its length must be equal to the the branches of the tree. This allows to specify a model with variable rates of evolution. You must be careful that branch numbering is done with the tree in ``postorder'' order: to see the order of the branches you can use: tr <- reorder(tr, "po"); plor(tr); edgelabels(). The arguments alpha and theta are ignored.

• "OU":an Ornstein-Uhlenbeck model is used. The above indexing rule is used for the three parameters sigma, alpha, and theta. This may be interesting for the last one to model varying phenotypic optima. The exact updating formula from Gillespie (1996) are used which are reduced to BM formula if alpha = 0.

• 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. The arguments sigma, alpha, and theta are ignored.

References

Gillespie, D. T. (1996) Exact numerical simulation of the Ornstein-Uhlenbeck process and its integral. Physical Review E, 54, 2084--2091.

Paradis, E. (2012) Analysis of Phylogenetics and Evolution with R (Second Edition). New York: Springer.

See Also

rTraitDisc, rTraitMult, ace

Examples

# NOT RUN {
data(bird.orders)
rTraitCont(bird.orders) # BM with sigma = 0.1
### OU model with two optima:
tr <- reorder(bird.orders, "postorder")
plot(tr)
edgelabels()
theta <- rep(0, Nedge(tr))
theta[c(1:4, 15:16, 23:24)] <- 2
## sensitive to 'alpha' and 'sigma':
rTraitCont(tr, "OU", theta = theta, alpha=.1, sigma=.01)
### an imaginary model with stasis 0.5 time unit after a node, then
### BM evolution with sigma = 0.1:
foo <- function(x, l) {
    if (l <= 0.5) return(x)
    x + (l - 0.5)*rnorm(1, 0, 0.1)
}
tr <- rcoal(20, br = runif)
rTraitCont(tr, foo, ancestor = TRUE)
### a cumulative Poisson process:
bar <- function(x, l) x + rpois(1, l)
(x <- rTraitCont(tr, bar, ancestor = TRUE))
plot(tr, show.tip.label = FALSE)
Y <- x[1:20]
A <- x[-(1:20)]
nodelabels(A)
tiplabels(Y)
# }