make.bisse.split {diversitree} | R Documentation |
Create a likelihood function for a BiSSE model where the tree is partitioned into regions with different parameters.
make.bisse.split(tree, states, nodes, split.t, unresolved=NULL, sampling.f=NULL, nt.extra=10, strict=TRUE, control=list())
tree |
An ultrametric bifurcating phylogenetic tree, in
|
states |
A vector of character states, each of which must be 0 or
1, or |
nodes |
Vector of nodes that will be split (see Details). |
split.t |
Vector of split times, same length as |
unresolved |
Unresolved clade information: see section below for structure. |
sampling.f |
Vector of length 2 with the estimated proportion of
extant species in state 0 and 1 that are included in the phylogeny.
A value of |
nt.extra |
The number of species modelled in unresolved clades (this is in addition to the largest observed clade). |
strict |
The |
control |
List of control parameters for the ODE solver. See
details in |
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).
TODO: Describe nodes
and split.t
here.
Richard G. FitzJohn
pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01)
set.seed(546)
phy <- tree.bisse(pars, max.taxa = 30, x0 = 0)
Here is the phylogeny:
plot(phy, show.node.label = TRUE, label.offset = 0.1, font = 1, cex = 0.75,
no.margin = TRUE)
Here is a plain BiSSE function for comparison:
lik.b <- make.bisse(phy, phy$tip.state)
lik.b(pars) # -93.62479
## [1] -93.62
Split this phylogeny at three points: nd15, nd18 and nd26
nodes <- c("nd15", "nd18", "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")
## Error: plot.new has not been called yet
To make a split BiSSE function, pass the node locations and times in:
lik.s <- make.bisse.split(phy, phy$tip.state, nodes.i)
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)
pars4
## [1] 0.10 0.20 0.03 0.03 0.01 0.01 0.10 0.20 0.03 0.03 0.01 0.01 0.10 0.20
## [15] 0.03 0.03 0.01 0.01 0.10 0.20 0.03 0.03 0.01 0.01
Run the likelihod calculation:
lik.s(pars4) # -93.62479
## [1] -93.62
These are basically identical (to acceptable tolerance)
lik.s(pars4) - lik.b(pars)
## [1] 3.847e-09
You can use the labelled nodes rather than indices:
lik.s2 <- make.bisse.split(phy, phy$tip.state, nodes)
identical(lik.s(pars4), lik.s2(pars4))
## [1] TRUE
This also works where some tips are unresolved clades. Here are a few:
unresolved <- data.frame(tip.label = c("sp12", "sp32", "sp9", "sp22",
"sp11"), Nc = c(2, 5, 3, 2, 5), n0 = c(1, 4, 3, 2, 4), n1 = c(1, 1, 0, 0,
1))
Plain BiSSE with unresolved clades:
lik.u.b <- make.bisse(phy, phy$tip.state, unresolved = unresolved)
lik.u.b(pars) # -139.3688
## [1] -139.4
Split BiSSE with unresolved clades:
lik.u.s <- make.bisse.split(phy, phy$tip.state, nodes, unresolved = unresolved)
lik.u.s(pars4) # -139.3688
## [1] -139.4
lik.u.b(pars) - lik.u.s(pars4) # numerical error only
## [1] -9.006e-09