castor (version 1.3.6)

fit_bm_model: Fit a Brownian motion model for multivariate trait evolution.

Description

Given a rooted phylogenetic tree and states of one or more continuous (numeric) traits on the tree's tips, fit a multivariate Brownian motion model of correlated co-evolution of these traits. This estimates a single diffusivity matrix, which describes the variance-covariance structure of each trait's random walk. The model assumes a fixed diffusivity matrix on the entire tree.

Usage

fit_bm_model(tree, tip_states, check_input=TRUE)

Arguments

tree

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge.

tip_states

A numeric vector of size Ntips, or a 2D numeric matrix of size Ntips x Ntraits, specifying the numeric state of each trait at each tip in the tree.

check_input

Logical, specifying whether to perform some basic checks on the validity of the input data. If you are certain that your input data are valid, you can set this to FALSE to reduce computation.

Value

A list with the following elements:

diffusivity

Either a single non-negative number (if tip_states was a vector) or a 2D quadratic non-negative-definite matrix (if tip_states was a 2D matrix). The fitted diffusivity matrix of the multivariate Brownian motion model.

loglikelihood

The log-likelihood of the fitted model, given the provided tip states data.

Details

The BM model is defined by the stochastic differential equation $$ dX = \sigma \cdot dW $$ where \(W\) is a multidimensional Wiener process with Ndegrees independent components and \(\sigma\) is a matrix of size Ntraits x Ndegrees. Alternatively, the same model can be defined as a Fokker-Planck equation for the probability density \(\rho\): $$ \frac{\partial \rho}{\partial t} = \sum_{i,j}D_{ij}\frac{\partial^2\rho}{\partial x_i\partial x_j}. $$ The matrix \(D\) is referred to as the diffusivity matrix (or diffusion tensor), and \(2D=\sigma\cdot\sigma^T\). Note that \(\sigma\) can be obtained from \(D\) by means of a Cholesky decomposition.

The function uses phylogenetic independent contrasts (Felsenstein, 1985) to retrieve independent increments of the multivariate random walk. The diffusivity matrix \(D\) is then fitted using maximum-likelihood on the intrinsic geometry of positive-definite matrices.

If tree$edge.length is missing, each edge in the tree is assumed to have length 1. The tree may include multifurcations (i.e. nodes with more than 2 children) as well as monofurcations (i.e. nodes with only one child). Note that multifurcations are internally expanded to bifurcations, prior to model fitting. The asymptotic time complexity of this function is O(Nedges*Nsimulations*Ntraits).

References

J. Felsenstein (1985). Phylogenies and the Comparative Method. The American Naturalist. 125:1-15.

See Also

simulate_bm_model, get_independent_contrasts

Examples

# NOT RUN {
# generate a random tree
tree = generate_random_tree(list(birth_rate_intercept=1), 10000)$tree


# Example 1: Scalar case
# - - - - - - - - - - - - - -
# simulate scalar continuous trait on the tree
D = 1
tip_states = simulate_bm_model(tree, diffusivity=D)$tip_states

# estimate original diffusivity from the generated data
fit = fit_bm_model(tree, tip_states)
cat(sprintf("True D=%g, fitted D=%g\n",D,fit$diffusivity))


# Example 2: Multivariate case
# - - - - - - - - - - - - - - -
# simulate vector-valued continuous trait on the tree
D = get_random_diffusivity_matrix(Ntraits=5)
tip_states = simulate_bm_model(tree, diffusivity=D)$tip_states

# estimate original diffusivity matrix from the generated data
fit = fit_bm_model(tree, tip_states)

# compare true and fitted diffusivity matrices
cat("True D:\n"); print(D)
cat("Fitted D:\n"); print(fit$diffusivity)
# }