R: Brownian Motion and Ornstein-Uhlenbeck Character Evolution

make.bm {diversitree}

R Documentation

Brownian Motion and Ornstein-Uhlenbeck Character Evolution

Description

Create a likelihood function for models of simple Brownian Motion or Ornstein-Uhlenbeck (OU) character evolution. This function creates a likelihood function that can be used in maximum likelihood or Bayesian inference.

Usage

make.bm(tree, states, states.sd=0, control=list())
make.ou(tree, states, states.sd=0, control=list())

Arguments

`tree`	A bifurcating phylogenetic tree, in `ape` “phylo” format.
`states`	A vector of continuous valued character states. This vector must be named with the tip labels of `tree`.
`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 `states`.
`control`	A list of control parameters. See details below.

Details

The control argument is a named list of options.

The main option is method. Specifying 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 for make.bm (but is not available for make.ou). The alternative is to specify control=list(method="pruning") which 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.

When method="vcv" 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 VCV-based make.bm function is heavily based on fitContinuous in the geiger package.

For non-ultrametric trees with OU models, computed likelihoods will differ because of the different root treatments.

Author(s)

Richard G. FitzJohn

Examples

Random data (following APE)

data(bird.orders)
set.seed(1)
x <- structure(rnorm(length(bird.orders$tip.label)), names = bird.orders$tip.label)

With the VCV approach

fit1 <- find.mle(make.bm(bird.orders, x), 0.1)

With the pruning calculations

lik.pruning <- make.bm(bird.orders, x, control = list(method = "pruning"))
fit2 <- find.mle(lik.pruning, 0.1)

All the same (need to drop the function from this though)

all.equal(fit1[-7], fit2[-7])

## [1] TRUE

If this is the same as the estimates from Geiger, to within the tolerances expected for the calculation and optimisation:

[This section is not run by default]

fit3 <- fitContinuous(bird.orders, x)
all.equal(fit3$Trait1$lnl, fit1$lnLik)
all.equal(fit3$Trait1$beta, fit1$par, check.attributes = FALSE)

[Ends not run section]

[Package diversitree version 0.9-4 ]