This function simulates stochastic birth-death trees. Simulation can be performed conditioning on n, on t, or on both simultaneously. If the both, then (for optional argument method="rejection") rejection sampling is performed whereby trees are simulated given b and t until a tree containing n taxa is found. The giving-up point can be set using the optional argument max.count. Simulations can also be performed in continuous time (the default) or discrete time; the difference being that wait times in the continous-time simulation come from the exponential distribution; whereas waiting times in discrete-time simulations come from the geometric distribution. In addition, discrete-time simulations allow for the possibility that multiple speciation events can occur at (exactly) the same time, so long as they are on separate branches. Finally, sometimes for stopping criterion n in discrete-time there will be a number of tips different from n. This indicates that the last event contained more than one speciation event, and a warning is printed.
method="direct" is presently experimental. It does not really perform direct sampling; however waiting times & birth or death events are sampled first - with only wait-times consistent with n and t being retained. This rejection sampling occurs one layer earlier than for method="rejection". This results in a significant (several-fold) speed-up of the code and enables sampling conditioned on n and t simultaneously for much higher b and d. At the present time, extant.only=TRUE does not work for this mode, nor does type="discrete".
Note that if ape=FALSE, then the function will run faster, and the tree is theoretically compatible with the ape "phylo" standard; however some downstream errors with functions such as bind.tree have been observed.
Lastly, under the taxon number stopping criterion (n) for a non-zero extinction rate (d>0) sometimes a tree containing fewer than n extant tips is returned because it has gone completely extinct before the end of the simulation.
pbtree(b=1, d=0, n=NULL, t=NULL, scale=NULL, nsim=1, type=c("continuous",
"discrete"), ...)A tree or list of trees as an object of class "phylo" or "multiPhylo", respectively.
birth rate or speciation rate for type="continuous"; the probability of speciating per time-step for type="discrete".
death rate or extinction rate for type="continuous"; the probability of going extinct per time-step for type="discrete".
desired number of species (i.e., taxa-stop criterion).
total time for simulation (i.e., time-stop criterion).
if set, rescales tree to have total length scale.
number of simulated trees to return.
string to indicate whether to simulate trees in continuous or discrete time. If the former, then wait times between speciation events are drawn from an exponential distribution; whereas if the latter then wait times comes from a geometric distribution.
optional arguments including ape, a logical value indicating whether to return nodes in a 'ape' compatible ordering (default is TRUE); extant.only a logical value indicating whether or not to return only extant species (defaults to FALSE); max.count a numeric value indicating the maximum number of iterations to run is sampling conditioned on both n and t (defaults to 1e5); method gives the method used for simultaneously conditioning on n and t - options are "rejection" and "direct"; tip.label, a vector of tip labels (only works for n!=NULL); and, finally, quiet, a logical value indicating whether or not to suppress certain message (defaults to FALSE).
Liam Revell liam.revell@umb.edu
Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol., 3, 217-223.
# simulate a pure-birth tree with 400 tips, scaled to a length of 1.0
tree<-pbtree(n=400,scale=1)
# simulate a pure-birth tree conditioning on n & t
tt<-log(50)-log(2)
tree<-pbtree(n=50,t=tt)
Run the code above in your browser using DataLab