In statistical phylogenetics, studying how species evolved helps people to understand evolution better. As many questions are arising from evolutionary biology and ecology, one interesting research question could be: how could traits of a group of related species behave to adapt the changing environment? For example, when studying marine species(Watanabe12052015, ), a scientist may be interested in understanding the moving speed and moving style by comparing fin structures in various kind of swordfish. One useful tool to track down their evolutionary information is incorporating a phylogenetic tree into analysis. A phylogenetic tree is a branching diagram that infers evolutionary relationships among a group of species. Given a tree
and traits (e.g. fin lengths or total lengths of fish in center-meter), we could use statistical approach to study ancestral status for species as well as how one trait could be related to the other trait. From mathematical perspective, changing of trait value or status during evolutionary history can be viewed as a stochastic random variable defined on time/status domain. In the case of continuous trait, letbe a trait of a species observed at time . The dynamic behavior of , when applied for studying trait evolution, can be assumed as a solution to the following stochastic differential equation (SDE)
In the left hand side of Eq. 1, represents the amount of change in an infinitesimal time . In the right hand side of Eq. 1, the deterministic term is referred to a drift coefficient that measures the amount of change in an infinitesimal time while is called the diffusion coefficient that amplifies the trait change according to the random changing environment measured by where is a Wiener process having continuous paths and independent Gaussian increments (i.e. ) and is the model parameters.
In literature, there have been statistical methods developed for traits evolution by applying continuous stochastic processes ranging from Gaussian process (felsenstein1985 ; hansen1996 ) or non-Gaussian processes (Blomberg067363, ; jhwueng2019, ). Currently one of the most popular continuous process for trait evolution can be credited to the Ornstein Uhlenbeck(OU) process (hansen1997, ). An OU stochastic variable solves the SDE in Eq. (1) with and . The OU process provides a suitable interpretation in describing natural selection in evolution and ecology context. The constant parameter is intepreted as the optimum status (evolutionary niche in ecology context) of . The parameter is called a constraining force that pulls trait back to the optimum . The parameter is called the rate of evolution and measures the speed of the random change.
Many works have been developed by expanding the OU model through considering more sophisticated and complex biological phenomenon. Those models used a generalized OU process to describe trait change along the tree. The generalized OU model for trait evolution is built by assuming pertinent processes for model parameters and . Therefore, the trait solves the following SDE:
Currently several works have been focus on the conditions by assuming as a constant, or as either a constant or with a stochastic dynamics during the evolutionary process (see butler2004 , omeara2006 , and beaulieu2012 ). By assuming following a pertinent process, solves the following SDE:
In particular, in the case of and , hansen08
created an OUBM model for optimal regression analysis built under the assumption that the optimumhas a linear relationship with predictors. Jhwueng2014 expanded the OUBM model to OUOU model by allowing an Ornstein-Uhlenbeck process for the dynamic of (i.e. and ). Those models are applied to study the adaptive relationship of traits building upon its optimum with where is a set of predictors, are regression parameters. See application sections in (hansen08, ; Jhwueng2014, ; jhwueng2016adaptive, ).
For the rate evolution in Eq. (2), instead of considering constant value or piecewise constant value (omeara2006, ), it is also reasonable to assume that the rate of evolution follows another pertinent process. Under this assumption, is a solution to another SDE: In literature, jhwueng2016adaptive considered the rate to be a Brownian motion where and is a constant.
In this work, observing that there are needs and possibilities to create models for more sophisticated and realistic biological applications, we expand previous existed models in two folds. First, as the rate is regarded as non-negative for , we intend to incorporate a Cox-Ingersoll-Ross(CIR) process (cox1985, ) for . In this case, solves the following SDE:
where is a constant force, is the optimum of and is the rate of change for . In CIR process, the distribution of future values of conditioned in current value has a distribution of where and and
is a non central chi-squared distribution. Notice that in Eq. (4), the diffusion coefficient involves a term which indicates that the Eq. (4) is neither a linear SDE nor
a normal distributed stochastic variable. Hence statistical inference on the parameter estimation under our new model will be different from the framework in(hansen08, ; Jhwueng2014, ) using the multivariate normal distribution for jointly modeling trait evolution. Secondly, we assume that there exists an interaction relationship between the optimum and predictors as following
where the term is the interaction between the th and the th predictors with regression parameter . Note that this model is different from the phylogenetic ancova model in fuen16 where the optimum is not considered with relationship to the predictors as shown in Eq. (5).
When jointly modeling adaptive trait evolution using Eqs. (2), (3), (4) and (5), the distribution for the trait of a species is constrained by the dynamic assumption of the rate parameter via either a Brownian motion case (jhwueng2016adaptive, ) or the CIR process case in Eq. (4).
However, given as a CIR process is not a Gaussian random variable and has intractable model likelihood. Conceiving this, we propose an algorithm under the approximation Bayesian computation(ABC) framework for statistical inference. We describe our framework into the following sections. Section 2 illustrates the general construction of adaptive model under a various of assumption of pertinent processes for and . We call our new models the OUBMCIR model for following generalized OU process with following a BM and following a CIR process. And we called the OUOUCIR model for following a generalized OU process with following an OU process and following a CIR process. Section 3 contains methods on simulating traits under each model. We make an attempt to derive the solution as explicitly as possible for the purpose of applying tree traversal algorithm (fel2004, ) to simulate trait status on the internal nodes and tips on the tree. We conduct statistical inference for parameter estimation under ABC in Section 4. Currently we mainly use the R package abc for inference after traits are simulated from Section 5. We provide empricial analsis on analyzing data from literature in Section 6. We conclude our study in Section 7. The scripts and their brief description developed in work project can be accessed at Github: https://github.com/djhwueng/ououcir.
2.1 Property of adaptive trait models
We start this section by first introducing some definitions of the SDE property. In Eq. (1), the SDE is a linear SDE if and are linear function of . That is, A linear SDE is autonomous if all coefficients are constants, is homogenous if and and is linear in the additive sense if .
These properties could provide some information on the distribution of . For instance, the SDE for in OUBM model (hansen08, ) with and is a linear additive non-autonomous SDE. In the OUBM model, since both and are BMs, the solution for the SDE in Eq. (1) is represented as a linear combination of two BMs. As dynamics of each BM can be treated as a normal random variable, we can conclude that is normal random variable in OUBM model. In this case, we can implement normal distribution to analyze data. We categorize the properties of SDE of as well as and in Table 1.
|OUBM||(✓, ✓, -)||(n, ✓, -)||(✓, ✓, -)||(✓, ✓, -)||(hansen08, )|
|OUOU||(✓, ✓, -)||(n, ✓, -)||(✓, ✓, -)||(✓, ✓, -)||(Jhwueng2014, )|
|OUBMBM||(✓, ✓, ✓)||(n, ✓, ✓)||(✓, ✓, ✓)||(n, ✓, ✓)||(jhwueng2016adaptive, )|
|OUOUBM||(✓, ✓, ✓)||(n, ✓, ✓)||(✓, ✓, ✓)||(n, ✓, ✓)||(jhwueng2016adaptive, )|
|OUBMCIR||(✓, ✓, n)||(n, ✓, n)||(✓, ✓, n)||(n, ✓, n)||This work|
|OUOUCIR||(✓, ✓, n)||(n, ✓, n)||(✓, ✓, n)||(n, ✓, n)||This work|
2.2 Solution of Model
into a system of SDE for the random vectoras , where is the drift vector, is the diffusion vector, and is the associated independent Brownian process random vector and is a transpose of a vector . By assuming that the force parameters are time invariant(), the model can be represented as
For homogeneous model assuming the rate of evolution as a time invariant constant (i.e. and in OUBM model and in OUOU model), we have and is a constant diagonal matrix. In this case, given the initital condition at , the system of SDE described by Eq. (6) has a unique solution . In this case, the expected value of can be calculated straightforwardly as
while the second moment of the random vector, denoted by
, can uniquely be determined by solving the system of an ordinary differential equationwhere . Once the first and second moment of are identified, because is a normal random vector, its first component is a normal random variable. We can also work from Eq. (2) on the assumption that is a constant. The solution is a linear combination of normal random variable which is again a normal random variable under the assumption of BM for (hansen08, ) or OU for (Jhwueng2014, ).
On the other hand, however, for OUBMBM, OUOUBM, OUBMCIR, and OUOUCIR model, as the rate of evolution follows a certain pertinent process, the distribution of is not as straightforward to work through. We show that fails to be a normal distributed random vector. We first demonstrate this using the new proposed OUOUCIR model with
Due to assumption of using CIR process for the rate parameter and the stationary distribution of a CIR random variable is not a normal random variable, the solution to the system of equation in Eq. (6) is intractible and not likely to be normal distribution.
Moreover, even for following a Brownian motion, we claims that fails to be a normal random variable. For the OUBMBM model in (jhwueng2016adaptive, ), the solution for the SDE in Eq. (2) under OUBMBM model is
where and are standard Wiener processes.
By direct calculation on the stochastic integral, we have which is a normal random variable with mean
and variance (by Itô’s isometry). However, the stochastic integral in 2⃝ fails to be a normal random variable. To show this, for simplicity we assume that and are two identical and independent Wiener processes. Let , by Itô’s lemma we have . Since neither nor is a normal random variable, 2⃝ fails to be a normal distributed. This indicates that to in Eq. (7) can not a normal random variable.
2.3 Multiple optimal regression with interaction
In this section, we describe how to implement the interaction in Eq. (5) into our model. To start, we use an example of two predictors for illustration. The general case can be extended accordingly. Given that the linear relationship between the optimum and predictors with interaction is , by differentiating on both side of the equation with respect to , we have
where is a diffusion process satisfies the SDE as following
By the SDE of in Eq. (3) and assumptions of stochastic calculus with , we have . In the case of assuming either a BM or an OU process, we have which implies . Similarly for in Eq. (9) for either BM or OU process, we have and . The relationship between and given the predictor traits and can be derived with expanding using Eq. (8) and represented as
The general case of optimum regression on the predictors with interaction can be extended from above with assumption with the form
By applying the same technique from above, we have
where is the correlation between two Wiener processes (i.e. ).
Then using the same technique on and compare it with , we have
Eq. (12) suggests that depends on the predictors s which are stochastic variable, in order to quantify , we consider to use expected value of .
When is a Brownian motion, since and , we have
When is an OU process, we have
where and .
3 Simulate trait along tree
Given a tree with known topology and length, we simulate tip as well as ancestral states using tree traversal algorithm (fel2004, ) under a specified model . In particular, when the distribution is known, for instance, under Brownian motion trait value of a species at time conditioned on its ancestor on is a normal random variable with mean and variance . (i.e. ). Under OU process, is a normal random variable with mean , and variance
. Moreover, under either BM or OU process the tip can be simulated directly under the joint distribution (i.e.where is the trait vector at tip of the tree, is the mean vector, and is the variance covariance structure for (jhwueng13, )).
Given the prior information on model parameters, our goal is to simulate th response trait , and predictor traits at the tip. We describe our method for simulating trait under each model using the given parameters values.
3.1 OUBM & OUOU model
For OUBM model, the model parameters are and , and regression parameters are . We first simulate predictor traits s on each node/tip of tree given . The optimal value can then be calculated via given and . Then use to simulate (see (hansen08, ) for the formula of and ).
For OUOU model, model parameters are , and , and regression parameters are We simulate predictor trait s on each node/tip of tree using . The optimal value can be calculated via to obtain on each nodes. We use to simulate , by (see (Jhwueng2014, ) for the formula of and ).
Note that since the OUBM model and OUOU model are both of multivariate normal distributions, trait values at tips can be simulated directly given the specified mean vector and variance structure .
3.2 OUBMBM model
In OUBMBM model, the model parameters are , and regression parameters are . We first simulate predictor traits s on each node/tip of tree given . The optimal value can then be calculated via given and . To simulate s at the nodes/tips, we first look at the solution in Eq. (2) for :
For 1⃝, as we assume the optimum follows Brownian motion (i.e. ), the term is a stochastic integral of Brownian motion with respect to time and equals to .
Since , we have the integral which is a normal random variable with mean and variance
In 2⃝, since the rate is assumed as BM (i.e. ), we have . Hence 2⃝ is a stochastic integral that involves an integral of Brownian motion with respect to another Brownian motion . Note 2⃝ is not a normal distributed random variable (see section 2.1). In order to draw sample from 2⃝, we use function int.st in R package Sim.DiffProc (guid17, ) to simulate the trajectory of this stochastic integral. We assume and are two independent and identical processes. We then use median of the trajectory as a sample for 2⃝. Given the parameter values, we can apply tree traversal algorithm to simulate sample on node/tip conditioned on its ancestor .
3.3 OUOUBM model
In OUOUBM model, model parameters are , and regression parameters . We first simulate predictor traits s on each node/tip of tree using . The optimum on each node and tip can be calculated as . To simulate s, since the solution in Eq. (2) under OUOUBM model is
For 1⃝, because is an OU process with where is optimum of and is the initital condition. The integral becomes
Note that a⃝ and b⃝ are both definite integrals with a⃝ = and b⃝ = . In c⃝, the term is a normal random variable with mean and variance . The integrand in c⃝ defined as is a normal random variable with mean and variance (by Itô Isometry) . So c⃝ = is again a normal random variable. Because is not an invertible function, it is not likely to identify the distribution of directly using change of variable. We alternatively use linear approximation for with at where and to obtain an candidate of distribution of which is a normal random variable with mean and variance .
For 2⃝, as the rate is a BM, we can simulate samples use the method for the 2⃝ described in the OUBMBM model.
3.4 OUBMCIR model
In OUBMCIR model, the model parameters are , and regression parameters are , . We first use to simulate predictor trait s and then use to obtain on each node/tip. To simulate s, since the solution in Eq. (2) is
For 1⃝, since the optimum is a BM (i.e ), we can draw using the expected value and variance as shown in the 1⃝ in the OUBMBM model.
For 2⃝, it is a stochastic integral of a CIR random variable with respect to Brownian motion . Note that follows a scaled non-central chi-squared distribution where and
is a non-central chi-squared distribution with degree of freedomand non-centrality parameter (jhwueng2019, ).
The distribution of the random variable conditioned on can be seen as the sum of three independent random variables (see prop. 4 Eq. 2.10 in chan10 ). Moreover, glasserman11 and chan10 showed that the exact distribution of , conditional on and can be representd by infinite sums and mixtures of gamma random variables (see prop 4. in (chan10, )). For our case, to simulate sample in 2⃝, we first simulate on each node along the tree using tree traversal as in (jhwueng2019, ). We next simulate sample for the random variable . Since the solution to the CIR SDE in Eq. (4) is given by