doSimulation: Discrete-Time Character Simulation

In bomeara/treevo: Using ABC to Understand Trait Evolution

Discrete-Time Character Simulation

Description

The function doSimulation evolves continuous characters under a discrete time process. These functions are mainly used as internal components, generating simulations within ABC analyses using the doRun functions. See Note below.

Usage

doSimulation(
  phy = NULL,
  intrinsicFn,
  extrinsicFn,
  startingValues,
  intrinsicValues,
  extrinsicValues,
  generation.time = 1000,
  TreeYears = max(branching.times(phy)) * 1e+06,
  timeStep = NULL,
  returnAll = FALSE,
  taxonDF = NULL,
  checkTimeStep = TRUE
)

Arguments

phy	A phylogenetic tree, in package ape's phylo format.
intrinsicFn	Name of (previously-defined) function that governs how traits evolve within a lineage, regardless of the trait values of other taxa.
extrinsicFn	Name of (previously-defined) function that governs how traits evolve within a lineage, based on their own ('internal') trait vlaue and the trait values of other taxa.
startingValues	State at the root.
intrinsicValues	Vector of values corresponding to the parameters of the intrinsic model.
extrinsicValues	Vector of values corresponding to the parameters of the extrinsic model.
generation.time	The number of years per generation. This sets the coarseness of the simulation; if it's set to 1000, for example, the population's trait values change every 1000 calendar years. Note that this is in calendar years (see description for argument TreeYears), and not in millions of years (as is typical for dated trees in macroevolutionary studies). Thus, if a branch is 1 million-year time-unit long, and a user applies the default generation.time = 1000, then 1000 evolutionary changes will be simulated along that branch. See documentation for doSimulation for further details.
TreeYears	The amount of calendar time from the root to the furthest tip. Most trees in macroevolutionary studies are dated with branch lengths in units of millions of years, and thus the default for this argument is max(branching.times(phy)) * 1e6. If your tree has the most recent tip at time zero (i.e., the modern day), this would be the same as the root age of the tree. If your branch lengths are not in millions of years, you should alter this argument. Otherwise, leave this argument alone. See documentation for doSimulation for further details.
timeStep	This value corresponds to the length of intervals between discrete evolutionary events ('generations') simulated along branches, relative to a rescaled tree where the root to furthest tip distance is 1. For example, timeStep = 0.01 of would mean 100 (i.e., 1 / 0.01) evolutionary changes would be expected to occur from the root to the furthest tip. (Note that the real number simulated will be much less, because simulations start over at each branching node.) Ideally, timeStep (or its effective value, via other arguments) should be as short as is computationally possible. Typically NULL by default and determined internally as follows: timeStep = generation.time / TreeYears. Can be provided a value as an alternative to using arguments generation.time and TreeYears, which would then be overridden. See documentation for doSimulation for further details.
returnAll	If TRUE, the output returned is a data.frame containing the values at each node from the simulation.
taxonDF	A data.frame containing data on nodes (both tips and internal nodes) output by various internal functions. Can be supplied as input to spead up repeated calculations, but by default is NULL, which instead forces a calculation from input phy.
checkTimeStep	If TRUE, warnings will be issued if TimeStep is too short.

Details

The phylogenetic tree used is rescaled such that the distance from the root to the furthest tip is rescaled to equal 1 time-unit, and it is this rescaled edge lengths to with arguments timeStep refers to. Typically, this will be determined though as a ratio of TreeYears (which is the number of calendar years constituing the root-to-furthest-tip distance, and is determined by default as if the user had supplied a tree with edge lengths in time-units of 1 time-unit = 1 million years), and generation.time, which gives the length of timeSteps in calendar years (e.g. generation.time = 1000 means an evolutionary change in trait values every 1000 years). Note that the real number of trait change events simulated may be less less, because simulations start over at each branching node, but intervals between changes should be so fine that this should be negligible (related, the results should be independent of your choice for generation.time or timeStep). We recommend that the effective timeStep should be as short as is computationally possible. SaveRealParams is useful for tracking the real true values if simulating data to test the performance of ABC analyses. It is not useful for ABC analyses of empirical data.

Value

If returnAll = FALSE (the default), this function returns a data frame of species character (tip) values in the tree, as a single column named 'states' with rownames reflecting the taxon names given in phy$tip.label. If returnAll = TRUE, the raw data.frame from the simulation will instead be returned.

Note

The simulateWithPriors functions are effectively the engine that powers the doRun functions, while the doSimulation functions are the pistons within the simulateWithPriors engine. In general, most users will just drive the car - they will just use doRun, but some users may want to use simulateWithPriors or doSimulation functions to do various simulations.

Author(s)

Brian O'Meara and Barb Banbury

Examples

set.seed(1)
tree <- rcoal(20)
# get realistic edge lengths
tree$edge.length <- tree$edge.length*20

#Simple Brownian motion

char <- doSimulation(
    phy = tree, 
    generation.time = 100000, 
    intrinsicFn = brownianIntrinsic, 
    extrinsicFn = nullExtrinsic, 
    startingValues = c(10), #root state
    intrinsicValues = c(0.01), 
    extrinsicValues = c(0))

#Character displacement model with minimum bound

char <- doSimulation(
    phy = tree, 
    generation.time = 100000, 
    intrinsicFn = boundaryMinIntrinsic, 
    extrinsicFn = ExponentiallyDecayingPushExtrinsic, 
    startingValues = c(10), #root state
    intrinsicValues = c(0.05, 10, 0.01), 
    extrinsicValues = c(0, .1, .25))

bomeara/treevo documentation built on Aug. 19, 2023, 6:52 p.m.