ape (version 2.6-2)

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 (2006, p. 151) for a brief introduction.

Usage

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

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.
linear
a logical indicating which parameterisation of the OU model to use (see details).

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 ``pruningwise'' order: to see the order of the branches you can use: tr <- reorder(tr, "p"); 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 more interesting for the last one to model varying phenotypic optima.

By default the following formula is used:

$$x_{t''} = x_{t'} - \alpha l (x_{t'} - \theta) + \sigma l \epsilon$$

where $l (= t'' - t')$ is the branch length, and $\epsilon \sim N(0, 1)$. If $\alpha > 1$, this may lead to chaotic oscillations. Thus an alternative parameterisation is used if linear = FALSE:

$$x_{t''} = x_{t'} - (1 - exp(-\alpha l)) * (x_{t'} - \theta) + \sigma l \epsilon$$}

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

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

See Also

rTraitDisc, ace

Examples

Run this code
data(bird.orders)
rTraitCont(bird.orders) # BM with sigma = 0.1
### OU model with two optima:
tr <- reorder(bird.orders, "p")
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)

Run the code above in your browser using DataCamp Workspace