simulate_phylproc {pcmabc} | R Documentation |
Simultaneously simulate a phylogeny and a trait evolving on it.
Description
The function does simulates a phylogeny and phenotypic dataset under user defined models of trait and phylogeny evolution. In particular the phenotype may influence the branching dynamics.
Usage
simulate_phylproc(tree.height, simul.params, X0, fbirth, fdeath=NULL,
fbirth.params=NULL, fdeath.params=NULL, fsimulphenotype="sde.yuima",
n.contemporary=-1, n.tips.total=1000, step=NULL)
Arguments
tree.height |
The height of the desired output tree. The simulated tree is conditioned to be of a certain height. |
simul.params |
The parameters of the stochastic model to simulate the phenotype. They should be passed as a named list. |
X0 |
The value of the ancestral state at the start of the tree. |
fbirth |
A function that returns the birth rate at a given moment.
The fist parameter of the function has to correspond to the value of the
phenotype (it is a vector first element is time and the others the values of the trait(s))
and the second to the list of birth parameters
|
fdeath |
A function that returns the birth rate at a given moment. Its structure
is the same as |
fbirth.params |
The parameters of the birth rate, to be provided to |
fdeath.params |
The parameters of the death rate, to be provided to |
fsimulphenotype |
The name of a function to simulate the phenotype over a period of time.
The function has to have four parameters (in the following order not by name):
The phenotype is simulated prior to the simulation of the speciation/extinction events on a lineage (see Details). Hence, it is not possible (at the moment) to include some special event (e.g. a cladogenetic jump) at branching. Such dynamics are only possible at the start of the lineage, i.e. when the new lineage separated from the main lineage. Hence, cladogenetic change can only be included as an event connected with a lineage (subpopulation) breaking off. |
n.contemporary |
The number of contemporary species to generate. If equals -1, then ignored.
Otherwise when the tree reaches |
n.tips.total |
The total (contemporary and extinct) number of tips to generate.
If equals -1, then ignored. Otherwise when the tree reaches
|
step |
The time step size for simulating the phenotype. If not provided,
then calculated as |
Details
The tree is simulate by means of a Cox process (i.e. Poisson process
with random rate). First the trait is simulated along the spine of a
tree, i.e. a lineage of duration tree.height
. Then, along this
spine the birth and death rates are calculated (they may depend
on the value of the phenotype). The maximum for each rate is calculated
and a homogeneous Poisson process with the maximum rate is simulated.
Then, these events are thinned. Each event is retained with probability
equalling true rate divided by maximum of rate (p. 32, Sheldon 2006).
All speciation events are retained until the first death event.
Value
tree |
The simulated tree in |
phenotype |
A list of the trajectory of the simulated phenotype.
The i-th entry of the list corresponds to the trait's evolution on
the i-th edge (as in i-th row of |
root.branch.phenotype |
The simulation on the root branch. A matrix with first row being time and next rows the simulated trait(s). |
Author(s)
Krzysztof Bartoszek
References
Bartoszek, K. and Lio', P (2019). Modelling trait dependent speciation with Approximate Bayesian Computation. Acta Physica Polonica B Proceedings Supplement 12(1):25-47.
Sheldon R. M. (2006). Simulation. Elsevier Academic Press.
See Also
Examples
## simulate 3d OUBM model with id branching rate
set.seed(12345)
simulate_mvsl_sde<-function(time,params,X0,step){
A <- c(paste("(-",params$a11,")*(x1-(",params$psi1,"))
-(",params$a12,")*(x2-(",params$psi2,"))-(",params$b11,")*x3",sep=""),
paste("(-",params$a21,")*(x1-(",params$psi1,"))
-(",params$a22,")*(x2-(",params$psi2,"))-(",params$b21,")*x3",sep=""),0)
S <- matrix( c( params$s11, params$s12, 0, 0, params$s22
, 0, 0, 0, params$s33), 3, 3,byrow=TRUE)
yuima.3d <- yuima::setModel(drift = A, diffusion = S,
state.variable=c("x1","x2","x3"),solve.variable=c("x1","x2","x3") )
simulate_sde_on_branch(time,yuima.3d,X0,step)
}
birth.params<-list(scale=1,maxval=2)
sde.params<-list(a11=2.569531,a12=0,a21=0,a22=28.2608,b11=-5.482939,
b21=-34.806936,s11=0.5513215,s12=1.059831,s22=1.247302,s33=1.181376,
psi1=-2.4590422,psi2=-0.6197838)
X0<-c(5.723548,4.103157,3.834698)
step<-0.25 ## for keeping example's running time short <5s as CRAN policy,
## in reality should be much smaller e.g. step<-0.001
simres<-simulate_phylproc(3.5, sde.params, X0, fbirth="rate_id", fdeath=NULL,
fbirth.params=NULL, fdeath.params=NULL, fsimulphenotype=simulate_mvsl_sde,
n.contemporary=5, n.tips.total=-1, step=step)