simulOUCHProcPhylTree: Simulate data on a phylogeny under a (multivariate) OU model

In mvSLOUCH: Multivariate Stochastic Linear Ornstein-Uhlenbeck Models for Phylogenetic Comparative Hypotheses

Simulate data on a phylogeny under a (multivariate) OU model

Description

Simulate data on a phylogeny under a (multivariate) OU model

Usage

simulOUCHProcPhylTree(phyltree, modelParams, regimes = NULL, 
regimes.times = NULL, dropInternal = TRUE, M.error=NULL, fullTrajectory=FALSE,
jumpsetup=NULL,keep_tree=FALSE,step=NULL)

Arguments

phyltree	The phylogeny in phylo format. The tree can be obtained from e.g. a nexus file by the read.nexus() function from the ape package. The "standard" ape node indexing is assumed: for a tree with n tips, the tips should have indices 1:n and the root index n+1. The root.edge field is ignored.
modelParams	List of model parameters of OUOU model as ParamsInModel part of output of ouchModel.
regimes	A vector or list of regimes. If vector then each entry corresponds to each of phyltree's branches, i.e. to each row of phyltree$edge.If list then each list entry corresponds to a tip node and is a vector for regimes on that lineage. If NULL, then a constant regime is assumed on the whole tree.
regimes.times	A list of vectors for each tree node, it starts with 0 and ends with the current time of the species. In between are the times where the regimes (niches) changed. If NULL then each branch is considered to be a regime.
dropInternal	Logical whether the simulated values at the internal nodes should be dropped.
M.error	An optional measurement error covariance structure. The measurement errors between species are assumed independent. The program tries to recognize the structure of the passed matrix and accepts the following possibilities : a single number that is a common measurement error for all tips and species, a m element vector with each value corresponding to a variable, measurement errors are independent between variables and each species is assumed to have the same measurement errors, a m x m ((number of variables) x (number of variables)) matrix, all species will have the same measurement error, a list of length n (number of species), each list element is the covariance structure for the appropriate (numbering according to tree) species, either a single number (each variable has same variance), vector (of length m for each variable), or m x m matrix, the order of the list has to correspond to the order of the nodes in the phyltree object, NULL no measurement error. From version 2.0.0 of mvSLOUCH it is impossible to pass a single joint measurement error matrix for all the species and traits.
fullTrajectory	should the full realization of the process or only node and tip values be returned
jumpsetup	Either NULL or list describing the jump at speciation. In the second case: jumptype In what way does the jump take place. Possible values are "ForBoth" the jump occurs at speciation and is common to both daughter lineages, "RandomLineage" the jump occurs just after speciation affecting exactly one daughter lineage, both descending branches have the same chance of being affected, "JumpWithProb" the jump occurs with probability jumpprob just after speciation independently on each daughter lineage independently. jumpprob A value in [0,1] indicating the probability of a jump taking place, only matters if jumptype is "JumpWithProb" or "JumpWithProb". jumpdistrib The distribution of the jump, currently only can take value "Normal". vMean The expected value of the jump, a vector of appropriate length if the trait is multivariate. mCov The variance of the jump, a matrix of appropriate dimensions if the trait is multivariate.
keep_tree	Logical whether the used tree should be saved inside the output object. Useful for any future reference, but as the tree is enhanced for mvSLOUCH's needs the resulting output object may be very large (it the number of tips is large).
step	The step size of the simulation.

Value

If fullTrajectory is FALSE then returns a matrix with each row corresponding to a tree node and each column to a trait. Otherwise returns a more complex object describing the full realization of the process on the tree. If dropInternal is TRUE, then the entries for the internal nodes are changed to NAs. The ordering of the rows corresponds to the order of the nodes (their indices) in the phylo object. Hence, the first n rows will be the tip rows (by common phylo convention).

Author(s)

Krzysztof Bartoszek

References

Bartoszek, K. (2014) Quantifying the effects of anagenetic and cladogenetic evolution. Mathematical Biosciences 254:42-57.

Bartoszek, K. (2016) A Central Limit Theorem for punctuated equilibrium. arXiv:1602.05189.

Bartoszek, K. and Pienaar, J. and Mostad. P. and Andersson, S. and Hansen, T. F. (2012) A phylogenetic comparative method for studying multivariate adaptation. Journal of Theoretical Biology 314:204-215.

Butler, M.A. and A.A. King (2004) Phylogenetic comparative analysis: a modeling approach for adaptive evolution. American Naturalist 164:683-695.

Hansen, T.F. (1997) Stabilizing selection and the comparative analysis of adaptation. Evolution 51:1341-1351.

Hansen, T.F. and Bartoszek, K. (2012) Interpreting the evolutionary regression: the interplay between observational and biological errors in phylogenetic comparative studies. Systematic Biology 61(3):413-425.

Pienaar et al (in prep) An overview of comparative methods for testing adaptation to external environments.

See Also

hansen, ouchModel, simulOUCHProcPhylTree

Examples

RNGversion(min(as.character(getRversion()),"3.6.1"))
set.seed(12345, kind = "Mersenne-Twister", normal.kind = "Inversion")
### We will first simulate a small phylogenetic tree using functions from ape.
### For simulating the tree one could also use alternative functions, e.g. sim.bd.taxa 
### from the TreeSim package
phyltree<-ape::rtree(5)

## The line below is not necessary but advisable for speed
phyltree<-phyltree_paths(phyltree)

### Define a vector of regimes.
regimes<-c("small","small","large","small","small","large","large","large")

### Define SDE parameters to be able to simulate data under the OUOU model.
OUOUparameters<-list(vY0=matrix(c(1,-1,0.5),nrow=3,ncol=1),
A=rbind(c(9,0,0),c(0,5,0),c(0,0,1)),mPsi=cbind("small"=c(1,-1,0.5),
"large"=c(-1,1,0.5)),Syy=rbind(c(1,0.25,0.3),c(0,1,0.2),c(0,0,1)))

### Now simulate the data.
jumpobj<-list(jumptype="RandomLineage",jumpprob=0.5,jumpdistrib="Normal",
vMean=rep(0,3),mCov=diag(1,3,3))
OUOUdata<-simulOUCHProcPhylTree(phyltree,OUOUparameters,regimes,NULL,jumpsetup=jumpobj)
RNGversion(as.character(getRversion()))

mvSLOUCH documentation built on Nov. 21, 2023, 1:08 a.m.