Create a likelihood function for a birth-death model where the tree is partitioned into regions with different parameters.
An ultrametric bifurcating phylogenetic tree, in
Vector of nodes that will be split (see Details).
Vector of split times, same length as
Probability of an extant species being included in the phylogeny (sampling fraction). By default, all extant species are assumed to be included.
Unresolved clade information. This is a named
vector, with the number of species as the value and names
corresponding to tip labels. Tips that represent a single species
should not be included in this vector. For example
Branching times can be controlled with the
argument. If this is
Inf, split at the base of the branch (as in
0, split at the top (closest to the present, as in
the new option for MEDUSA). If
0 < split.t < Inf then we split
at that time on the tree (zero is the present, with time growing
This function is related to MEDUSA (Alfaro et al. 2009), but does not
include any of the code for efficiently moving between different splits
(split creation here is fairly slow). The primary use for this model is
for generating starting points for state dependent split models (e.g.,
make.bisse.split) and testing a priori splits.
Richard G. FitzJohn
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
set.seed(1) pars <- c(.1, .03) phy <- trees(pars, "bd", max.taxa=30)[] ## Here is the phylogeny: plot(phy, show.node.label=TRUE, label.offset=.1, font=1, cex=.75, no.margin=TRUE) ## Construct the plain likelihood function as a benchmark: lik <- make.bd(phy) lik(pars) # -21.74554 ## Split this phylogeny at three points: nd11, nd13 and nd26 nodes <- c("nd11", "nd13", "nd26") ## This is the index in ape's node indexing: nodes.i <- match(nodes, phy$node.label) + length(phy$tip.label) nodelabels(node=nodes.i, pch=19, cex=2, col="#FF000099") ## To make a split likelihood function, pass the node locations and times in: lik.s <- make.bd.split(phy, nodes) ## The parameters must be a list of the same length as the number of ## partitions. Partition '1' is the root partition, and partition i is ## the partition rooted at the node[i-1] pars4 <- rep(pars, 4) names(pars4) <- argnames(lik.s) ## Run the likelihod calculation: lik.s(pars4) # -21.74554 ## These are basically identical (to acceptable tolerance) lik.s(pars4) - lik(pars) ## You can use the labelled nodes rather than indices: lik.s2 <- make.bd.split(phy, nodes) identical(lik.s(pars4), lik.s2(pars4)) ## All the usual ML/MCMC functions work as before: fit <- find.mle(lik.s, pars4)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.