Create a likelihood function for models of simple Brownian Motion (BM), Ornstein-Uhlenbeck (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.
1 2 3 4
A bifurcating phylogenetic tree, in
A vector of continuous valued character states. This
vector must be named with the tip labels of
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
Should the optimum (often "theta") be considered
a free parameter? The default,
A list of control parameters. See details below.
control argument is a named list of options.
The main option is
control=list(method="vcv") uses a variance-covariance 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.
uses Freckleton (2012)'s contrasts based approach, which is also much
faster on large trees.
method="pruning" is specified,
backend="C" may also be provided, which switch between a slow
(and stable) R calculator and a fast (but less extensively tested) C
backend="R" is currently the default.
The VCV-based functions are heavily based on
For non-ultrametric 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/ou-nonultrametric/ for a discussion about calculations on non-ultrametric trees.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## 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.