simulate: Evolve Birth-Death Trees

simulateR Documentation

Evolve Birth-Death Trees


Evolves one or more trees under the BiSSE (Binary State Speciation and Extinction), MuSSE (Multi-State Speciation and Extinction), BiSSE-ness (BiSSE-node enhanced state shift), ClaSSE (Cladogenetic State change Speciation and Extinction), or GeoSSE (Geographic State Speciation and Extinction) model, or a simple character independent birth-death model. For the SSE models, it simultaneously evolves a character that affects speciation and/or extinction, and the tree itself.


trees(pars, type=c("bisse", "bisseness", "bd", "classe", "geosse",
      "musse", "quasse", "yule"), n=1, max.taxa=Inf, max.t=Inf,
      include.extinct=FALSE, ...)

tree.bisse(pars, max.taxa=Inf, max.t=Inf, include.extinct=FALSE,
tree.musse(pars, max.taxa=Inf, max.t=Inf, include.extinct=FALSE,
tree.musse.multitrait(pars, n.trait, depth, max.taxa=Inf, max.t=Inf,
                      include.extinct=FALSE, x0=NA)

tree.quasse(pars, max.taxa=Inf, max.t=Inf, include.extinct=FALSE, x0=NA,
            single.lineage=TRUE, verbose=FALSE)

tree.bisseness(pars, max.taxa=Inf, max.t=Inf, include.extinct=FALSE,

tree.classe(pars, max.taxa=Inf, max.t=Inf, include.extinct=FALSE,

tree.geosse(pars, max.taxa=Inf, max.t=Inf, include.extinct=FALSE,
           x0=NA), max.taxa=Inf, max.t=Inf, include.extinct=FALSE)
tree.yule(pars, max.taxa=Inf, max.t=Inf, include.extinct=FALSE)

prune(phy, to.drop=NULL)



Vector of parameters. The parameters must be in the same order as an unconstrained likelihood function returned by make.x, for tree type x. The MuSSE simulator automatically detects the appropriate number of states, given a parameter vector.


Type of tree to generate: May be "bisse" or "bd".


How many trees to generate?


Maximum number of taxa to include in the tree. If Inf, then the tree will be evolved until max.t time has passed.


Maximum length to evolve the phylogeny over. If Inf (the default), then the tree will evolve until max.taxa extant taxa are present.


Logical: should extinct taxa be included in the final phylogeny? And should extinct trees be returned by trees?


Initial character state at the root (state 0 or 1). A value of NA will randomly choose a state from the model's equilibrium distribution for a BiSSE, ClaSSE, or GeoSSE model, but a non-NA value must be specified for MuSSE and QuaSSE.

n.trait, depth

For tree.musse.multitrait only, these specify the number of binary traits and the style of parameters (with the same meaning as in make.musse.multitrait). The pars argument then needs to be in the same order as a likelihood function created by make.musse.multitrait with these arguments (this interface may be improved in future – email me if you find this annoying).


(tree.quasse only): Start simulation with a single lineage? If FALSE, then the simulation starts with two lineages in state x0 (i.e., immediately following a speciation event).


(tree.quasse only): print verbose details about tree simuations. This can be reassuring for really large trees.


Additional arguments


A phylogeny, possibly with extinct species, produced by one of the tree evolving functions.


Optional vector with the species names to drop.


The phylogeny will begin from a single lineage in state x0, but the final phylogeny will include only branches above the first split.

tree.bisse may return an extinct phylogeny, and trees might return extinct phylogenies if include.extinct is TRUE.


A phylo phylogenetic tree (ape format), or for bisse.trees, a list of phylo trees.

The trees will have an element tip.state that contains the binary state information.


There are some logic problems around the creation of zero and one species trees; this will cause occasional errors when running the above functions. Things will change to fix this soon. All these functions may change in the near future.


Richard G. FitzJohn


## Due to a change in sample() behaviour in newer R it is necessary to
## use an older algorithm to replicate the previous examples
if (getRversion() >= "3.6.0") {
  RNGkind(sample.kind = "Rounding")

pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01)
phy <- tree.bisse(pars, max.taxa=30, x0=0)

h <- history.from.sim.discrete(phy, 0:1)
plot(h, phy)

## Retain extinct species:
phy2 <- tree.bisse(pars, max.taxa=30, x0=0, include.extinct=TRUE)
h2 <- history.from.sim.discrete(phy2, 0:1)
plot(h2, phy2)

#### MuSSE:
## Two states
pars <- c(.1, .2, .03, .04, 0.05, 0.1)
phy <- tree.musse(pars, 20, x0=1)

h <- history.from.sim.discrete(phy, 1:2)
plot(h, phy)

## A 3-state example where movement is only allowed between neighbouring
## states (1 <-> 2 <-> 3), and where speciation and extinction rates
## increase moving from 1 -> 2 -> 3:
pars <- c(.1,  .15,  .2,  # lambda 1, 2, 3
          .03, .045, .06, # mu 1, 2, 3
          .05, 0,         # q12, q13
          .05, .05,       # q21, q23
          0,   .05)       # q31, q32

phy <- tree.musse(pars, 30, x0=1, include.extinct=TRUE)

h <- history.from.sim.discrete(phy, 1:3)
plot(h, phy, cex=.7)

## And with extinct taxa pruned:
phy2 <- prune(phy)
h2 <- history.from.sim.discrete(phy2, 1:3)
plot(h2, phy2, cex=.7)

## This can all be done in one step (and is by default):
phy <- tree.musse(pars, 30, x0=1)
h <- history.from.sim.discrete(phy, 1:3)
plot(h, phy, cex=.7)

diversitree documentation built on May 29, 2024, 4:38 a.m.