make.bd.t  R Documentation 
Create a likelihood function for the birthdeath model, where birth and/or death rates are arbitrary functions of time.
make.bd.t(tree, functions, sampling.f=NULL, unresolved=NULL,
control=list(), truncate=FALSE, spline.data=NULL)
tree 
An ultrametric bifurcating phylogenetic tree, in

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 
control 
List of control parameters for the ODE solver. See
details in 
truncate 
Logical, indicating if functions should be truncated to zero when negative (rather than causing an error). May be scalar (applying to all functions) or a vector (of length 2). 
spline.data 
List of data for splinebased time functions. See details 
.
Richard G. FitzJohn
## First, show equivalence to the plain Birthdeath 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(.1, .03)
phy < trees(pars, "bd", max.taxa=25)[[1]]
## Next, make three different likelihood functions: a "normal" one that
## uses the direct birthdeath calculation, an "ode" based one (that
## uses numerical integration to compute the likelihood, and is
## therefore not exact), and one that is timevarying, but that the
## timedependent functions are constant.t().
lik.direct < make.bd(phy)
lik.ode < make.bd(phy, control=list(method="ode"))
lik.t < make.bd.t(phy, c("constant.t", "constant.t"))
lik.direct(pars) # 22.50267
## ODEbased likelihood calculations are correct to about 1e6.
lik.direct(pars)  lik.ode(pars)
## The ODE calculation agrees exactly with the timevarying (but
## constant) calculation.
lik.ode(pars)  lik.t(pars)
## Next, make a real case, where speciation is a linear function of
## time.
lik.t2 < make.bd.t(phy, c("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)
## The time penalty comes from moving to the ODEbased solution, not
## from the time dependence.
system.time(lik.direct(pars)) # ~ 0.000
system.time(lik.ode(pars)) # ~ 0.003
system.time(lik.t(pars)) # ~ 0.003
system.time(lik.t2(pars2)) # ~ 0.003
## Not run:
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)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.