simulate_bm_model {castor} | R Documentation |

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.

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

`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 (" |

`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 |

`include_tips` |
Include random states for the tips. If |

`include_nodes` |
Include random states for the nodes. If |

`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` |
Number of random independent simulations to perform. For each node and/or tip, there will be |

`drop_dims` |
Logical, specifying whether singleton dimensions should be dropped from |

The BM model for Ntraits co-evolving traits is defined by the stochastic differential equation

*
dX = σ \cdot dW
*

where *W* is a multidimensional Wiener process with Ndegrees independent components and *σ* is a matrix of size Ntraits x Ndegrees. Alternatively, the same model can be defined as a Fokker-Planck equation for the probability density *ρ*:

*
\frac{\partial ρ}{\partial t} = ∑_{i,j}D_{ij}\frac{\partial^2ρ}{\partial x_i\partial x_j}.
*

The matrix *D* is referred to as the diffusivity matrix (or diffusion tensor), and *2D=σ\cdotσ^T*. Either `diffusivity`

(*D*) or `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).

A list with the following elements:

`tip_states` |
Either |

`node_states` |
Either |

Stilianos Louca

`simulate_ou_model`

,
`simulate_rou_model`

,
`simulate_mk_model`

,
`fit_bm_model`

# 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)

[Package *castor* version 1.6.7 Index]