make.bd {diversitree}  R Documentation 
Prepare to run a constant rate birthdeath model on a
phylogenetic tree. This fits the Nee et al. 1994 equation,
duplicating the birthdeath
function in ape. Differences with
that function include (1) the function is not constrained to positive
diversification rates (mu can exceed lambda), (2) [eventual] support
for both random taxon sampling and unresolved terminal clades (but see
bd.ext
), and (3) run both MCMC and MLE fits to birth death
trees.
make.bd(tree, sampling.f=NULL, unresolved=NULL, times=NULL, control=list())
make.yule(tree, sampling.f=NULL, unresolved=NULL, times=NULL, control=list())
starting.point.bd(tree, yule=FALSE)
tree 
An ultrametric bifurcating phylogenetic tree, in

times 
Vector of branching times, as returned by

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

yule 
Should the starting point function return a Yule model (zero extinction rate)? 
control 
List of control parameters. The element 
make.bd
returns a function of class bd
.
This function has argument list (and default values)
f(pars, prior=NULL, condition.surv=TRUE)
The arguments are interpreted as
pars
A vector of two parameters, in the order
lambda
, mu
.
prior
: a valid prior. See make.prior
for
more information.
condition.surv
(logical): should the likelihood
calculation condition on survival of two lineages and the speciation
event subtending them? This is done by default, following Nee et
al. 1994.
The function "ode" method is included for completeness, but should not be taken too seriously. It uses an alternative ODEbased approach, more similar to most diversitree models, to compute the likelihood. It exists so that other models that extend the birthdeath models may be tested.
Richard G. FitzJohn
Nee S., May R.M., and Harvey P.H. 1994. The reconstructed evolutionary process. Philos. Trans. R. Soc. Lond. B Biol. Sci. 344:305311.
constrain
for making submodels, find.mle
for ML parameter estimation, mcmc
for MCMC integration,
and make.bisse
for statedependent birthdeath models.
## Simulate a tree under a constant rates birthdeath model and look at
## the maximum likelihood speciation/extinction parameters:
set.seed(1)
phy < trees(c(.1, .03), "bd", max.taxa=25)[[1]]
lik < make.bd(phy)
## By default, optimisation gives a lambda close to 0.1 and extremely
## small mu:
fit < find.mle(lik, c(.1, .03))
coef(fit)
## The above optimisation uses the algorithm \link{nlm} for
## compatibility with ape's \link{birthdeath}. This can be slightly
## improved by using \link{optim} for the optimisation, which allows
## bounds to be specified:
fit.o < find.mle(lik, c(.1, .03), method="optim", lower=0)
coef(fit.o)
logLik(fit.o)  logLik(fit) # slight improvement
## Special case methods are worked out for the Yule model, for which
## analytic solutions are available. Compare a direct fit of the Yule
## model with one where mu is constrained to be zero:
lik.yule < make.yule(phy)
lik.mu0 < constrain(lik, mu ~ 0)
## The same to a reasonable tolerance:
fit.yule < find.mle(lik.yule, .1)
fit.mu0 < find.mle(lik.mu0, .1)
all.equal(fit.yule[1:2], fit.mu0[1:2], tolerance=1e6)
## There is no significant improvement in the fit by including the mu
## parameter (unsurprising as the ML value was zero)
anova(fit.o, yule=fit.yule)
## Optimisation can be done without conditioning on survival:
fit.nosurv < find.mle(lik, c(.1, .03), method="optim", lower=0,
condition.surv=FALSE)
coef(fit.nosurv) # higher lambda than before
## Look at the marginal likelihoods, computed through MCMC (see
## \link{mcmc} for details, and increase nsteps for smoother
## plots [takes longer]).
samples < mcmc(lik, fit$par, nsteps=500,
lower=c(Inf, Inf), upper=c(Inf, Inf), w=c(.1, .1),
fail.value=Inf, print.every=100)
samples$r < with(samples, lambda  mu)
## Plot the profiles (see \link{profiles.plot}).
## The vertical lines are the simulated parameters, which match fairly
## well with the estimated ones.
col < c("red", "blue", "green3")
profiles.plot(samples[c("lambda", "mu", "r")], col.line=col, las=1,
legend="topright")
abline(v=0, lty=2)
abline(v=c(.1, .03, .07), col=col)
## Sample the phylogeny to include 20 of the species, and run the
## likelihood search assuming random sampling:
set.seed(1)
phy2 < drop.tip(phy, sample(25, 5))
lik2 < make.bd(phy2, sampling.f=20/25)
fit2 < find.mle(lik2, c(.1, .03))
## The ODE based version gives comparable results. However, it is
## about 55x slower.
lik.ode < make.bd(phy, control=list(method="ode"))
all.equal(lik.ode(coef(fit)), lik(coef(fit)), tolerance=2e7)