View source: R/sim.trait.values.R
sim.trait.values | R Documentation |
Fossil recovery rates or other parameter values can be simulated for a phylo (tree
) or taxonomy (taxonomy
) object.
Under the autocorrelated
model, trait values evolve along lineages according to a Brownian motion process, where the strength of the relationship between ancestor and descendant values is determined by the parameter ν (v
).
If ν is small values will be more similar between ancestor and descendants, and if ν is zero all trait values will be equal.
For a given species i with ancestor j, a new trait value κ_i is drawn from a lognormal distribution with
κ_i ~ LN( ln([κ_j] - (σ^2/2), σ)
where σ = ν * t_i and t_i is the lineage duration of the species.
This fossil recovery model is described in Heath et al. (2014) and is equivalent to the autocorrelated relaxed clock model described in Kishino et al. (2001).
Under the BM
and OU
models traits are simulated under a standard Brownian motion or Ornstein-Uhlenbeck process with rate parameter ν (v
).
The OU model has the additional parameter alpha
, which determines the strength with which trait values are attracted to the mean. Note the init
argument will specify both the value at the root and the mean of the process under the OU model.
Under the independent
model a new trait value is drawn for each species from any valid user-specified distribution (dist
).
change.pr
is the probability that a trait value will change at each speciation event.
If change.pr = 1
trait values will be updated at every speciation events.
Finally, traits can be simulated under the standard Lewis Mk model (Mk
), with symmetric rates of change. The rate is specified using v
and number of states using k
.
sim.trait.values( init = 1, tree = NULL, taxonomy = NULL, root.edge = TRUE, model = "autocorrelated", v = 0.01, alpha = 0.1, min.value = -Inf, max.value = Inf, dist = function() { runif(1, 0, 2) }, change.pr = 1, k = 2 )
init |
Initial value at the origin or root of the phylo or taxonomy object. Default = 1. |
tree |
Phylo object. |
taxonomy |
Taxonomy object. |
root.edge |
If TRUE include the root edge. Default = TRUE. |
model |
Model used to simulate rate variation across lineages. Options include "autocorrelated" (default), "BM" (Brownian motion), "OU" (Ornstein-Uhlenbeck), "independent" or the Lewis "Mk" model. |
v |
Brownian motion parameter v used in the autocorrelated, BM and OU models. Or rate change under the Mk model. Default = 0.01. |
alpha |
Ornstein-Uhlenbeck parameter alpha. Determines the strength with which trait values are pulled back towards the mean. |
min.value |
Min trait value allowed under the BM and OU models. Default = -Inf. |
max.value |
Max trait value allowed under the BM and OU models. Default = Inf. |
dist |
Distribution of trait values used to draw new values under the "independent" model. This parameter is ignored if |
change.pr |
Probability that trait values change at speciation events. Default = 1. |
k |
Number of states used for the Mk model. Default = 2. |
A vector of parameter values.
Values are output for each species in the order in which they appear in the taxonomy object (if taxonomy was provided) or for each edge in the order in which they appear in the tree object.
If the tree object has a root edge (root.edge
), the first entry in the vector will correspond to this edge.
Heath et al. 2014. The fossilized birth-death process for coherent calibration of divergence-time estimates. PNAS 111:E2957-E2966.
Kishino et al. 2001. Performance of a divergence time estimation method under a probabilistic model of rate evolution MBE 18:352-361.
# simulate tree t = ape::rtree(6) # simulate taxonomy s = sim.taxonomy(t, 0.5, 1, 0.5) # simulate rates under the autocorrelated trait values model rate = 2 rates = sim.trait.values(rate, taxonomy = s, v = 1) f = sim.fossils.poisson(rates, taxonomy = s) plot(f, t) # simulate rates under the independent trait values model dist = function() { rlnorm(1, log(rate), 1) } rates = sim.trait.values(rate, taxonomy = s, model = "independent", dist = dist) f = sim.fossils.poisson(rates, taxonomy = s) plot(f, t) # simulate rates under the independent trait values model with infrequent changes rates = sim.trait.values(rate, taxonomy = s, model = "independent", dist = dist, change.pr = 0.1) f = sim.fossils.poisson(rates, taxonomy = s) plot(f, t) # simulate traits under Brownian motion and convert into rates traits = sim.trait.values(0, taxonomy = s, model = "BM", v = 2) # function for translating states into rates translate.states = function(traits, low, high) sapply(traits, function(t) if(t < 0) low else high) # sampling rates low = 0.1 high = 2 rates = translate.states(traits, low, high) f = sim.fossils.poisson(rates, taxonomy = s) plot(f, tree = t)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.