simulate_phylproc {pcmabc}

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 birth.params, see par0. The time entry of the phenotype vector is at the moment relative to an unspecified (from fbirth's perspective) speciation event on the phylogeny. Hence, it cannot be used as for writing a time-inhomogenous speciation function. The speciation process is assumed to be time homogeneous in the current implementation. The package has support for two inbuilt rate functions. The can be indicated by passing a string in fbirth: either "rate_id" or "rate_const". The string "rate_const" corresponds to a constant rate and has to have the rate's value in field called $rate. However, a switching of rates is allowed. If the value of the first trait exceeds a certain threshold (provided in field $switch of birth parameters), then the rate is changed to the value in $rate2, see body of hidden function .f_rate_const(), in file rates.R. The string "rate_id" corresponds to the .f_rate_id() function, in file rates.R. If the birth parameters are NULL, then the rate equals the value of the first phenotype. However, a number of linear, threshold and power transformations of the rate are possible.The field varnum indicates the index of the variable to take as the one influencing the rate (remember to add 1 for the time entry). Then, if x stands for the trait influencing the branching rate it is transformed into a rate by the following fields in the following order. Set rate<-x and let params correspond to the list containing the branching parameters.

• substractbaserate<-max(0,params$rate-substractbase) rate<-abs(rate)

• p and if is.null(params$raise.at.end) || !params$raise.at.end rate<-rate^params$p

• base and !is.null(params$const)
if (res<params$base){rate<-params$const}

• base and is.null(params$const)
if (res<params$base){rate<-0}

• invbase and !is.null(params$const)
if (res>params$invbase){rate<-params$const}

• invbase and is.null(params$const)
if (res>params$invbase){rate<-0}

• scalerate<-rate/params$scale

• p and if params$raise.at.endrate<-rate^params$p res<-abs(res)

• maxvalif(rate>params$maxval){rate<-params$maxval}

fdeath

A function that returns the birth rate at a given moment. Its structure is the same as fbirth. The current version of the package does not provide any support for any inbuilt function.

fbirth.params

The parameters of the birth rate, to be provided to fbirth. They have to be a named list.

fdeath.params

The parameters of the death rate, to be provided to fdeath. They have to be a named list.

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): time, params, X0 and step. The parameter time is the duration of the simulation, i.e. the length of the branch. The parameter params is a list of the parameters governing the simulation, what will be passed here is the list
phenotype.model.params, without the fields fixed, abstepsd and
positivevars. It is up to the function in phenotype.model to interpret the provided parameters. X0 stands for the value at the start of the branch and step is a control parameter governing the time step size of the simulation. The pcmabc package has inbuilt support for simulating the trait as as stochastic differential equation by the yuima package with the function simulate_sde_on_branch(). However the user needs to write the phenotype.model function that translates the vector of parameters into an object that is understandable by yuima and then call simulate_sde_on_branch(), see the Example.

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.contemporary tips at height tree.height the simulation is stopped. However, there is no conditioning on this value, it is just an upper bound. It may just happen that due to the birth and death rate functions the process stops before reaching the target number of tips.

n.tips.total

The total (contemporary and extinct) number of tips to generate. If equals -1, then ignored. Otherwise when the tree reaches n.tips.total tips the simulation is stopped. However, there is no conditioning on this value, it is just an upper bound. It may just happen that due to the birth and death rate functions the process stops before reaching the target number of tips.

step

The time step size for simulating the phenotype. If not provided, then calculated as min(0.001,tree.height/1000).

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 phylo format.

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 phyltree$edge) of the tree. Each entry is a matrix with first row equalling time (relative to the start of the branch) and the next rows correspond to the trait value(s).

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

PCM_ABC

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)