rTraitCont {ape}R Documentation

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.

Details

There are three possibilities to specify model:

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.

Author(s)

Emmanuel Paradis

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

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)

[Package ape version 5.5 Index]