simulate_phylproc: Simultaneously simulate a phylogeny and a trait evolving on...
In pcmabc: Approximate Bayesian Computations for Phylogenetic Comparative Methods

View source: R/simtree.R

simulate_phylproc

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 `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. `substractbase`rate<-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}` `scale`rate<-rate/params$scale `p` and if `params$raise.at.end`rate<-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.

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)

pcmabc documentation built on May 9, 2022, 9:09 a.m.