View source: R/modelbdsplit.R
make.bd.split  R Documentation 
Create a likelihood function for a birthdeath model where the tree is partitioned into regions with different parameters.
make.bd.split(tree, nodes, split.t, sampling.f=NULL, unresolved=NULL)
tree 
An ultrametric bifurcating phylogenetic tree, in

nodes 
Vector of nodes that will be split (see Details). 
split.t 
Vector of split times, same length as 
sampling.f 
Probability of an extant species being included in the phylogeny (sampling fraction). By default, all extant species are assumed to be included. 
unresolved 
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 split.t
argument. If this is Inf
, split at the base of the branch (as in
MEDUSA). If 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
backwards).
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
set.seed(1)
pars < c(.1, .03)
phy < trees(pars, "bd", max.taxa=30)[[1]]
## 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[i1]
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.