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.