R: Time-varing Birth-Death Models

make.bd.t {diversitree}

R Documentation

Time-varing Birth-Death Models

Description

Create a likelihood function for the birth-death model, where birth and/or death rates are arbitrary functions of time.

Usage

make.bd.t(tree, functions, sampling.f=NULL, unresolved=NULL,
          control=list())

Arguments

`tree`	An ultrametric bifurcating phylogenetic tree, in `ape` “phylo” format.
`functions`	A named list of functions of time. See details.
`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`	Not yet included: present in the argument list for future compatibility with `make.bd`.
`control`	List of control parameters for the ODE solver. See details in `make.bisse`.

Details

There is significant overhead in using this function, compared with the plain birth-death (make.bd) models. Usually, all of the calculations are done directly, or in C for the ODE based models. However, because of the generality here – parameters may be any function of time – we must call back to R many times during the integration process. The overhead is something like a factor of 10 compared with the ODE-based constant rate birth death model, and a factor of 40 compared with the direct-calculation constant rate birth death model.

Author(s)

Richard G. FitzJohn

Examples

First, show equivalence to the plain Birth-death model. This is not a very interesting use of the functions, but it serves as a useful check.

Here is a simulated 25 species tree for testing.

set.seed(1)
pars <- c(0.1, 0.03)
phy <- trees(pars, "bd", max.taxa = 25)[[1]]

Next, make three different likelihood functions: a "normal" one that uses the direct birth-death calculation, an "ode" based one (that uses numerical integration to compute the likelihood, and is therefore not exact), and one that is time-varying, but that the time-dependent functions are constant.t().

lik.direct <- make.bd(phy)
lik.ode <- make.bd(phy, control = list(method = "ode"))
lik.t <- make.bd.t(phy, list(constant.t, constant.t))

lik.direct(pars)  # -22.50267

## [1] -22.5

ODE-based likelihood calculations are correct to about 1e-6.

lik.direct(pars) - lik.ode(pars)

## [1] -3.391e-06

The ODE calculation agrees exactly with the time-varying (but constant) calculation.

lik.ode(pars) - lik.t(pars)

## [1] 0

Next, make a real case, where speciation is a linear function of time.

lik.t2 <- make.bd.t(phy, list(linear.t, constant.t))

Confirm that this agrees with the previous calculations when the slope is zero

pars2 <- c(pars[1], 0, pars[2])
lik.t2(pars2) - lik.t(pars)

## [1] 0

However, notice that this is much slower than the previous versions:

system.time(lik.direct(pars))  # ~ 0.001

##    user  system elapsed 
##   0.001   0.000   0.000

system.time(lik.ode(pars))  # ~ 0.004

##    user  system elapsed 
##   0.002   0.000   0.003

system.time(lik.t(pars))  # ~ 0.020

##    user  system elapsed 
##   0.010   0.000   0.009

system.time(lik.t2(pars2))  # ~ 0.040

##    user  system elapsed 
##   0.027   0.001   0.027


fit <- find.mle(lik.direct, pars)
fit.t2 <- find.mle(lik.t2, pars2)

No significant improvement in model fit:

anova(fit, time.varying = fit.t2)

##              Df lnLik  AIC ChiSq Pr(>|Chi|)
## minimal       2 -22.1 48.2                 
## time.varying  3 -21.9 49.9 0.341       0.56

[Package diversitree version 0.9-4 ]