Cladogenetic State change Speciation and Extinction Model

Prepare to run ClaSSE (Cladogenetic State change Speciation and Extinction) on a phylogenetic tree and character distribution. This function creates a likelihood function that can be used in maximum likelihood or Bayesian inference.

make.classe(tree, states, k, sampling.f=NULL, strict=TRUE,
              control=list())
  starting.point.classe(tree, k, eps=0.5)

Arguments

tree

An ultrametric bifurcating phylogenetic tree, in ape “phylo” format.

states

A vector of character states, each of which must be an integer between 1 and k. This vector must have names that correspond to the tip labels in the phylogenetic tree (tree$tip.label).

k

The number of states. (The maximum now is 31, but that can easily be increased if necessary.)

sampling.f

Vector of length k where sampling.f[i] is the proportion of species in state i that are present in the phylogeny. A value of c(0.5, 0.75, 1) means that half of species in state 1, three quarters of species in state 2, and all species in state 3 are included in the phylogeny. By default all species are assumed to be known

strict

The states vector is always checked to make sure that the values are integers on 1:k. If strict is TRUE (the default), then the additional check is made that every state is present. The likelihood models tend to be poorly behaved where states are missing, but there are cases (missing intermediate states for meristic characters) where allowing such models may be useful.

control

List of control parameters for the ODE solver. See details in make.bisse.

eps

Ratio of extinction to speciation rates to be used when choosing a starting set of parameters. The procedure used is based on Magallon & Sanderson (2001).

Details

The ClaSSE model with k = 2 is equivalent to but a different parameterization than the BiSSE-ness model. The GeoSSE model can be constructed from ClaSSE with k = 3; see the example below.

make.classe returns a function of class classe. The arguments and default values for this function are:

    f(pars, condition.surv=TRUE, root=ROOT.OBS, root.p=NULL,
      intermediates=FALSE)
  

The arguments of this function are explained in make.bisse. The speciation rate parameters are lambda_ijk, ordered with k changing fastest and insisting on j < k.

With more than 9 states, lambda_ijk and q_ij can be ambiguous (e.g. is q113 1->13 or 11->3?). To avoid this, the numbers are zero padded (so that the above would be q0113 or q1103 for 1->13 and 11->3 respectively). It might be easier to rename the arguments in practice though. More human-friendly handling of large speciation rate arrays is in the works.

starting.point.classe produces a first-guess set of parameters, ignoring character states.

Unresolved clade methods are not available for ClaSSE.

Tree simulation methods are not yet available for ClaSSE.

See also

constrain for making submodels, find.mle for ML parameter estimation, and mcmc for MCMC integration. The help page for find.mle has further examples of ML searches on full and constrained BiSSE models. Things work similarly for ClaSSE, just with different speciation parameters.

make.bisse, make.bisseness, make.geosse, make.musse for similar models and further relevant examples.

References

FitzJohn R.G., Maddison W.P., and Otto S.P. 2009. Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies. Syst. Biol. 58:595-611.

Goldberg E.E. and Igic B. Tempo and mode in plant breeding system evolution. In review.

Maddison W.P., Midford P.E., and Otto S.P. 2007. Estimating a binary character's effect on speciation and extinction. Syst. Biol. 56:701-710.

Magallon S. and Sanderson M.J. 2001. Absolute diversification rates in angiospem clades. Evol. 55:1762-1780.

Magnuson-Ford, K., and Otto, S.P. 2012. Linking the investigations of character evolution and species diversification. American Naturalist, in press.

Author

Emma E. Goldberg

Examples

## 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")
}
#> Warning: non-uniform 'Rounding' sampler used

## GeoSSE equivalence
## Same tree simulated in ?make.geosse
pars <- c(1.5, 0.5, 1.0, 0.7, 0.7, 2.5, 0.5)
names(pars) <- diversitree:::default.argnames.geosse()
set.seed(5)
phy <- tree.geosse(pars, max.t=4, x0=0)

lik.g <- make.geosse(phy, phy$tip.state)
pars.g <- c(1.5, 0.5, 1.0, 0.7, 0.7, 1.4, 1.3)
names(pars.g) <- argnames(lik.g)

lik.c <- make.classe(phy, phy$tip.state+1, 3)
pars.c <- 0 * starting.point.classe(phy, 3)
pars.c['lambda222'] <- pars.c['lambda112'] <- pars.g['sA']
pars.c['lambda333'] <- pars.c['lambda113'] <- pars.g['sB']
pars.c['lambda123'] <-  pars.g['sAB']
pars.c['mu2'] <- pars.c['q13'] <- pars.g['xA']
pars.c['mu3'] <- pars.c['q12'] <- pars.g['xB']
pars.c['q21'] <- pars.g['dA']
pars.c['q31'] <- pars.g['dB']

lik.g(pars.g)   # -175.7685
#> [1] -175.7685
lik.c(pars.c)   # -175.7685
#> [1] -175.7685