R/paleobuddy.R

In paleobuddy: Simulating Diversification Dynamics

#' paleobuddy: Simulating diversification dynamics
#'
#' \code{paleobuddy} provides users with flexible scenarios for species 
#' birth-death simulations. It also provides the possibility of generating 
#' phylogenetic trees (with extinct and extant species) and fossil records 
#' (with a number of preservation scenarios) from the same underlying process. 
#' 
#' @section Birth-death simulation:
#' Users have access to a large array of scenarios to use and combine for 
#' species birth-death simulation. The function \code{bd.sim} allows for
#' constant rates, rates varying as a function of time, or time and/or an 
#' environmental variable, as well as age-dependent rates by using a shape 
#' parameter from a Weibull distribution (which can itself also be 
#' time-dependent). Extinction and speciation rates can be supplied 
#' independently, so that one can combine any types of scenarios for birth and 
#' death rates. The function \code{find.lineages} separates birth-death 
#' simulations into monophyletic clades so one can generate fossil records and 
#' phylogenies (see below) for clades with a specific mother species. This is
#' particularly useful for simulations with multiple starting species. See 
#' \code{?bd.sim} and \code{?find.lineages} for more information.
#' 
#' All birth-death simulation functions return a \code{sim} object, which is a 
#' list of vectors containing speciation times, extinction times, status
#' (extant or extinct) and parent identity for each species of the simulation.
#' We supply methods for summarizing and printing \code{sim} objects in a more
#' informative manner (see Visualization below). See \code{?sim} for more 
#' information.
#' 
#' @section Fossil record simulation:
#' The package provides users with a similarly diverse array of scenarios for 
#' preservation rates in generating fossil records from birth-death 
#' simulations. The function \code{sample.clade} accepts constant, 
#' time-varying, and environmentally dependent rates. Users might also supply a
#' model describing the distribution of fossil occurrences over a species 
#' duration to simulate age-dependent sampling. See \code{?sample.clade} for 
#' more information.
#' 
#' @section Phylogeny generation:
#' We believe it is imperative to be able to generate fossil records and
#' phylogenetic trees from the same underlying process, so the package provides
#' \code{make.phylo}, a function that takes a simulation object of the form 
#' returned by \code{bd.sim} and generates a \code{phylo} object from the APE
#' package. One can then use functions such as \code{ape::plot.phylo} and
#' \code{ape::drop.fossil} to plot the phylogeny or analyze the
#' phylogeny of extant species. Since APE is not required for any function in 
#' the package, it is a suggested but not imported package. Note that, as 
#' above, the function \code{find.lineages} allows users to separate clades 
#' with mother species of choice, the results of which can be passed to
#' \code{make.phylo} to generate separate phylogenies for each clade. See 
#' \code{?make.phylo} and \code{?find.lineages} for more information.
#' 
#' Note: If a user wishes to perform the opposite operation - transform a
#' \code{phylo} object into a \code{sim} object, perhaps to use paleobuddy for
#' sampling on phylogenies generated by other packages, see 
#' \code{?phylo.to.sim}.
#' 
#' @section Visualization:
#' \code{paleobuddy} provides the user with a number of options for visualizing
#' a \code{sim} object besides phylogenies. The \code{sim} object returned by 
#' birth-death simulation functions (see above) has summary and plot methods. 
#' \code{summary(sim)} gives quantitative details of a \code{sim} objective, 
#' namely the total and extant number of species, and summaries of species
#' durations and speciation times. \code{plot(sim)} plots births, deaths, and
#' diversity through time for that realization. The function \code{draw.sim}
#' draws longevities of species in the simulation, allowing for customization
#' through the addition of fossil occurrences (which can be time points or
#' ranges), and vertical order of the drawn longevities.
#' 
#' @section Utility functions:
#' The package makes use of a few helper functions for simulation and testing
#' that we make available for the user. \code{rexp.var} aims to emulate the
#' behavior of \code{rexp}, the native R function for drawing an exponentially
#' distributed variate, with the possibility of supplying a time-varying
#' rate. The function also allows for a shape parameter, in which case the 
#' times drawn will be distributed as a Weibull, possibly with time-varying
#' parameters, for age-dependent rates. \code{var.rate.div} calculates the 
#' expected diversity of a birth-death process with varying rates for any time 
#' period, which is useful when testing the birth-death simulation functions. 
#' Finally, \code{binner} takes a vector of fossil occurrence times and a 
#' vector of time boundaries and returns the number of occurrences within each
#' time period. This is mostly for use in the \code{sample.clade} function. See
#' \code{?rexp.var}, \code{?var.rate.div} and \code{?binner} for more 
#' information.
#' 
#' @author 
#' Bruno do Rosario Petrucci, Matheus Januario and Tiago B. Quental
#' 
#' Maintainer: Bruno do Rosario Petrucci <petrucci@iastate.edu>
#' 
#' @examples 
#' 
#' # here we present a quick example of paleobuddy usage
#' # for a more involved introduction, see the \code{overview} vignette
#'
#' # make a vector for time
#' time <- seq(0, 10, 0.1)
#' 
#' # speciation rate
#' lambda <- function(t) {
#'   0.15 + 0.03*t
#' }
#' 
#' # extinction rate
#' mu <- 0.08
#' 
#' # these are pretty simple scenarios, of course
#' # check the examples in ?bd.sim for a more comprehensive review
#' 
#' # diversification
#' d <- function(t) {
#'   lambda(t) - mu
#' }
#' 
#' # calculate how many species we expect over 10 million years
#' div <- var.rate.div(rate = d, n0 = 1, t = time)
#' # note we are starting with 3 species (n0 = 3), but the user
#' # can provide any value - the most common scenario is n0 = 1
#' 
#' # plot it
#' plot(time, rev(div), type = 'l', main = "Expected diversity",
#'      xlab = "Time (My)", ylab = "Species",
#'      xlim = c(max(time), min(time)))
#' 
#' # we then expect around 9 species 
#' # alive by the present, seems pretty good
#' 
#' # set seed
#' set.seed(1)
#' 
#' # run the simulation
#' sim <- bd.sim(n0 = 1, lambda = lambda, mu = mu, 
#'               tMax = 10, nFinal = c(20, Inf))
#' # nFinal controls the final number of species
#' # here we throw away simulations with less than 20 species generated
#' 
#' # draw longevities
#' draw.sim(sim)
#'
#' # from sim, we can create fossil records for each species
#' # rho is the fossil sampling rate, see ?sample.clade
#' samp <- sample.clade(sim = sim, rho = 0.75, tMax = 10,
#'                      bins = seq(10, 0, -1))
#' # note 7 out of the 31 species did not leave a fossil - we can in this way
#' # simulate the incompleteness of the fossil record
#' 
#' # we can draw fossil occurrences as well, and order by extinction time
#' draw.sim(sim, fossils = samp, sortBy = "TE")
#' 
#' # take a look at the phylogeny
#' if (requireNamespace("ape", quietly = TRUE)) {
#'   ape::plot.phylo(make.phylo(sim), root.edge = TRUE)
#'   ape::axisPhylo()
#' }
#' 
#' @docType package
#' @name paleobuddy
NULL