simulate_bm_model {castor}

Simulate a Brownian motion model for multivariate trait co-evolution.

Description

Given a rooted phylogenetic tree and a Brownian motion (BM) model for the co-evolution of one or more continuous (numeric) unbounded traits, simulate random outcomes of the model on all nodes and/or tips of the tree. The function traverses nodes from root to tips and randomly assigns a multivariate state to each node or tip based on its parent's previously assigned state and the specified model parameters. The generated states have joint distributions consistent with the multivariate BM model. Optionally, multiple independent simulations can be performed using the same model.

Usage

simulate_bm_model(tree, diffusivity=NULL, sigma=NULL,
                  include_tips=TRUE, include_nodes=TRUE, 
                  root_states=NULL, Nsimulations=1, drop_dims=TRUE)

Arguments

tree

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

diffusivity

Either NULL, or a single number, or a 2D quadratic positive definite symmetric matrix of size Ntraits x Ntraits. Diffusivity matrix ("D") of the multivariate Brownian motion model (in units trait^2/edge_length). The convention is that if the root's state is fixed, then the covariance matrix of a node's state at distance L from the root will be 2LD (see mathematical details below).

sigma

Either NULL, or a single number, or a 2D matrix of size Ntraits x Ndegrees, where Ndegrees refers to the degrees of freedom of the model. Noise-amplitude coefficients of the multivariate Brownian motion model (in units trait/sqrt(edge_length)). This can be used as an alternative way to specify the Brownian motion model instead of through the diffusivity D. Note that sigma\cdot\sigma^T=2D (see mathematical details below).

include_tips

Include random states for the tips. If FALSE, no states will be returned for tips.

include_nodes

Include random states for the nodes. If FALSE, no states will be returned for nodes.

root_states

Numeric matrix of size NR x Ntraits (where NR can be arbitrary), specifying the state of the root for each simulation. If NR is smaller than Nsimulations, values in root_states are recycled in rotation. If root_states is NULL or empty, then the root state is set to 0 for all traits in all simulations.

Nsimulations

Number of random independent simulations to perform. For each node and/or tip, there will be Nsimulations random states generated.

drop_dims

Logical, specifying whether singleton dimensions should be dropped from tip_states and node_states, if Nsimulations==1 and/or Ntraits==1. If drop_dims==FALSE, then tip_states and tip_nodes will always be 3D matrices.

Details

The BM model for Ntraits co-evolving traits 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. Either diffusivity (D) or sigma (\sigma) may be used to specify the BM model, but not both.

If tree$edge.length is missing, each edge in the tree is assumed to have length 1. The tree may include multi-furcations (i.e. nodes with more than 2 children) as well as mono-furcations (i.e. nodes with only one child). The asymptotic time complexity of this function is O(Nedges*Nsimulations*Ntraits).

Value

A list with the following elements:

tip_states

Either NULL (if include_tips==FALSE), or a 3D numeric matrix of size Nsimulations x Ntips x Ntraits. The [r,c,i]-th entry of this matrix will be the state of trait i at tip c generated by the r-th simulation. If drop_dims==TRUE and Nsimulations==1 and Ntraits==1, then tip_states will be a vector.

node_states

Either NULL (if include_nodes==FALSE), or a 3D numeric matrix of size Nsimulations x Nnodes x Ntraits. The [r,c,i]-th entry of this matrix will be the state of trait i at node c generated by the r-th simulation. If drop_dims==TRUE and Nsimulations==1 and Ntraits==1, then node_states will be a vector.

Author(s)

Stilianos Louca

See Also

simulate_ou_model, simulate_rou_model, simulate_mk_model, fit_bm_model

Examples

# generate a random tree
tree = generate_random_tree(list(birth_rate_intercept=1),max_tips=10000)$tree

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

# plot histogram of simulated tip states
hist(tip_states, breaks=20, xlab="state", main="Trait probability distribution", prob=TRUE)

# Example 2: Multivariate case
# - - - - - - - - - - - - - - -
# simulate co-evolution of 2 traits with 3 degrees of freedom
Ntraits  = 2
Ndegrees = 3
sigma    = matrix(stats::rnorm(n=Ntraits*Ndegrees, mean=0, sd=1), ncol=Ndegrees)
tip_states = simulate_bm_model(tree, sigma=sigma, drop_dims=TRUE)$tip_states

# generate scatterplot of traits across tips
plot(tip_states[,1],tip_states[,2],xlab="trait 1",ylab="trait 2",cex=0.5)