This is the first release of paleobuddy, an R package
dedicated to flexible simulations of diversification, fossil records,
and phylogenetic trees. Below we list current features, and above
sections will be filled as new features and fixes are implemented.
rexp.var generalizes rexp to take any
function of time as a rate. Also allows for a shape
parameter, in which case it similarly generalizes
rweibull.bd.sim simulates species diversification with high
flexibility in allowed speciation and extinction rate scenarios.
Produces a sim object.sample.clade simulates fossil sampling. Similar
flexibility to bd.sim, though even more so in the case of
age-dependent rates.make.phylo creates a bifurcating phylogenetic tree as
aphylo object (see APE) from a
sim object. Can take fossils to be added as length
0 branches.draw.sim plots a sim by drawing species
durations and relationships, and optionally adding fossils as time
points or ranges.find.lineages creates subsets of a sim
defined as the clades descended from one or more species present in the
simulation. Needed to e.g. generate phylogenetic trees from simulations
with more than one starting species.make.rate creates a purely time-dependent (or constant)
rate based on optional inputs. Used internally to allow for users to
define rate scenarios easily in bd.sim.phylo.to.sim creates a sim object from a
phylo object, provided the user makes choices to solve
ambiguities on bifurcating phylogenetic trees.var.rate.div calculates expected diversity for a given
diversification rate and set of times. Useful for testing
bd.sim and planning rate scenarios.sim a class returned by bd.sim and used as
an input for many functions in the package. Formally, it is a named list
of vectors recording speciation time, extinction time, status (extant or
extinct), and parent information for each lineage in the simulation. It
contains the following methods: ** print gives some quick
details about number of extant and total species, and the first few
members of each vector. ** head and tail
return the sim object containing only a given number of
species from the beginning and end of its vectors, respectively. **
summary gives quantitative details, e.g. quantiles of
durations and speciation waiting times. ** plot plots
lineages through time (LTT) plots for births, deaths, and diversity. **
sim.counts counts numbers of births, deaths, and diversity
for some given time t. ** is.sim checks the
object is a valid sim object. Used internally for error
checking.temp temperature data during the Cenozoic. Modified
from RPANDA.co2 CO2 data during the Jurassic. Modified from RPANDA.overview gives a reasonably in depth look at the main
features of the package, including examples of workflows using most
available rate scenarios, and examples of applications.The question of how to structure time came up a lot during
development of the package. Most of the literature in macroevolution and
paleontology considers absolute geological time, i.e. t = 0
at the present and t = 5 five million years ago. It becomes
challenging, however, to visualize and program complex rates going
backwards. As such, the code is structured such that all functions are
considered to go from 0 to the maximum simulation time
tMax—i.e. the inverse of absolute geological time. There is
only one exception to this rule, the adFun parameter
describing age-dependent distribution of fossil occurrences in
sample.clade. In any case, all returned objects in the
package are set to follow absolute geological time, so as to conform to
the literature.
The integrate function occasionally fails when given
functions that vary suddenly and rashly, usually happening in the case
of environmentally-dependent rates. I have tracked this error down to
numerical problems in integrate, and testing seems to
indicate the error does not prevent integrate from finding
the correct result. As such this is not currently something I intend to
fix, though if issues are found that indicate this could be a
paleobuddy problem, not an integrate problem,
that could change.
paleobuddy is the first package to implement
time-dependent parameters for Weibull-distributed waiting times. Since
the authors currently are not aware of an analytical solution to
important quantities in the BD process in this case, it is challenging
to test exactly. Simulation tests indicate pretty strongly that the
algorithm works, however, with one exception—in cases where shape is
time-dependent and varies dramatically, especially when close to
0, rexp.var seems to have a hard time finding
the correct waiting time distribution. When maintained within levels
generally accepted as sensible throughout the literature—around 0.8 to
3, say—, and even a reasonable amount outside of that range, tests
indicate the algorithm functions as it should. A testthat
routine will be implemented in the future to formalize these claims, and
this issue is one I plan to work on soon, especially if users report it
as more prevalent than I thought.