make.bm  R Documentation 
Create a likelihood function for models of simple Brownian Motion (BM), OrnsteinUhlenbeck (OU), or Early Burst (EB) character evolution, or BM on a “lambda” rescaled tree. This function creates a likelihood function that can be used in maximum likelihood or Bayesian inference.
make.bm(tree, states, states.sd=0, control=list())
make.ou(tree, states, states.sd=0, with.optimum=FALSE, control=list())
make.eb(tree, states, states.sd=0, control=list())
make.lambda(tree, states, states.sd=0, control=list())
tree 
A bifurcating phylogenetic tree, in 
states 
A vector of continuous valued character states. This
vector must be named with the tip labels of 
states.sd 
An optional vector of measurement errors, as standard
deviation around the mean. If a single value is given it is used
for all tips, otherwise the vector must be named as for

with.optimum 
Should the optimum (often "theta") be considered
a free parameter? The default, 
control 
A list of control parameters. See details below. 
The control
argument is a named list of options.
The main option is method
. Specifying
control=list(method="vcv")
uses a variancecovariance matrix
based approach to compute the likelihood. This is similar to the
approach used by geiger, and is the default.
Two alternative approaches are available.
control=list(method="pruning")
uses the transition density
function for brownian motion along each branch, similar to how most
methods in diversitree are computed. This second approach is much
faster for very large trees. control=list(method="contrasts")
uses Freckleton (2012)'s contrasts based approach, which is also much
faster on large trees.
When method="pruning"
is specified, backend="R"
or
backend="C"
may also be provided, which switch between a slow
(and stable) R calculator and a fast (but less extensively tested) C
calculator. backend="R"
is currently the default.
The VCVbased functions are heavily based on fitContinuous
in
the geiger
package.
For nonultrametric trees with OU models, computed likelihoods may differ because of the different root treatments. This is particularly the case for models where the optimum is estimated.
For the EB model, the parameter intepretation follows geiger; the 'a' parameter is equivalent to log(g) in Bloomberg et al. 2003; when negative it indicates a decelerating rate of trait evolution over time. When zero, it reduces to Brownian motion.
Richard G. FitzJohn
See https://www.zoology.ubc.ca/prog/diversitree/examples/ounonultrametric/ for a discussion about calculations on nonultrametric trees.
## Random data (following APE)
data(bird.orders)
set.seed(1)
x < structure(rnorm(length(bird.orders$tip.label)),
names=bird.orders$tip.label)
## Not run:
## With the VCV approach
fit1 < find.mle(make.bm(bird.orders, x), .1)
## With the pruning calculations
lik.pruning < make.bm(bird.orders, x, control=list(method="pruning"))
fit2 < find.mle(lik.pruning, .1)
## All the same (need to drop the function from this though)
all.equal(fit1[7], fit2[7])
## If this is the same as the estimates from Geiger, to within the
## tolerances expected for the calculation and optimisation:
fit3 < fitContinuous(bird.orders, x)
all.equal(fit3$Trait1$lnl, fit1$lnLik)
all.equal(fit3$Trait1$beta, fit1$par, check.attributes=FALSE)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.