On the Application of the Autoregressive–Moving-Average Model for Studying Phylogenetic Rate of Trait Evolution

1. Introduction

In evolution, studying the heritable characteristics of the biological population is helpful to understand the diversity between species on our planet Earth [1]. Macroevolution is expected to occur when selection acts on a trait that has a heritable basis of phenotypic variation. During the evolutionary process, speciation results in new species, and the comparison of traits (e.g. height, weight, size, ⋯ etc.) among a group of related species can be made by studying the speed of changes in their characteristics over successive generations [2].

Although one subgroup of species evolved at a faster rate and resulted in a larger variation in the trait, the other subgroup of species evolved at a relatively lower rate and produced a moderate variation of the trait. When evolutionary processes such as natural selection (including sexual selection) and genetic drift act on this variation, certain characteristics become more common or rare within a population [3]. For example, Darwin finches are a group of about 18 species of dull-colored passerine birds on the Galápagos islands [4]. They are well known for their remarkable diversity in form, size, and function of the beak, which is highly adapted to different food sources [5]. Another example is angiosperms, which survive and thrive successfully on our planet. The evolution of fruits is one of the most important characteristics, as fruits not only provide a food source for other species, but also protect seeds and contribute to seed dispersal [6]. The survival and success of fruits require adaptation to their environment. For example, while dragon fruit (Selenicereus) endures temperatures up to 40°C (104°F) for survival, watermelon (Citrullus lanatus) needs temperatures higher than about 25°C (77°F) to thrive. Studying the reproducible properties by the rate of evolution to the resistance of temperature would help us shed light on the evolution of angiosperms themselves and understand their ecological implications.

The rate of evolution is a measurement of the change in an evolutionary lineage over time and can be defined as the ratio of the character displacement over a certain time interval. [7] defined the rate of change between two samples using three quantities: the proportional difference between the sample means, the pooled standard deviation of the samples, and the time interval between the samples. For example, suppose that a character has been measured twice, $t_1$ and $t_2$, where $t_1$ and $t_2$ are expressed as the time before the present in millions of years. The time interval between the two samples can be written as $dt=t_1-t_2$, which is 1 million years if $t_1=6.27$ and $t_2=5.27$. The average value of the character is defined as $x_1$ in the previous sample and $x_2$ in the later sample. Let $y_{t_i}=log\left(x_{t_i}\right),i=1,2$ be the natural logarithm taking into account $x_{t_1}$ and $x_{t_2}$. Then the evolutionary rate ($\sigma$) can be defined as in Eq. (1)

$$\sigma =(y_{t_2}-y_{t_1})/(t_2-t_1).$$

(1)

Next, multiplying the time difference $t_2-t_1$ on both sides of Eq. (1) produces the character difference within a time unit $y_{t_2}-y_{t_1}=\sigma (t_2-t_1).$ Conceptually consider that the character change occurred in infinitesimal time and denote the character displacement by $d{y}_t={lim}_{t_1\to t_2}(y_{t_2}-y_{t_1})={lim}_{t_1\to t_2}\sigma (t_2-t_1)$, then we have the differential equation $d{y}_t=\sigma dt$. Given $y\left(0\right)=y_0$, the solution is $y_t=y_0+\sigma t$ which shows that $y_t$ increased with time t. Unfortunately, this may not be an appropriate model for describing character change in the evolutionary perspective. Instead, one may consider that the variation of the character change adopts a certain dynamic. For example, if one considers that the variation of the character change is proportional to time, then a stochastic variable $W_t$ can be introduced, where $W_t$ is a Wiener process with independent Gaussian distributed increment and $d{W}_t=W_{t+dt}-W_t$ is a normal distributed random variable with mean 0 and variance $dt$ (i.e. $d{W}_t\sim \mathcal{N}(0,dt)$). Thus, the displacement of the continuous trait variable $y_t$ ($=log{x}_t$) solves stochastic differential equation shown in Eq. (2)

$$d{y}_t=\sigma d{W}_t.$$

(2)

Given $y\left(0\right)=y_0$, integrate both sides of Eq. (2), one has

$$y_t=y_0+\sigma W_t,$$

(3)

where $W_0=0$, and $W_t\sim \mathcal{N}(0,t)$.

Given the initial value $y_0$, Eq. (3) describes the dynamic of the trait variable $y_t$ at time t pending the rate parameter $\sigma$. Figure 1 presents one hundred trajectories generated using rates $\sigma =1$ (bottom left) and $\sigma =4$ (bottom right), respectively. It can be seen straightforwardly that while a smaller rate ($\sigma =1$) yields a narrower range of character value, the larger rate ($\sigma =4$) yields a wider range of character value.

For a group of n related species, denote $y_{i,t},i=1,2,\cdots ,n$ as the trait variable for the ith species. Then we can apply the Brownian motion with the rate parameter to explore the dynamic of the trait along the evolutionary history using a phylgoenetic tree $\mathcal{T}$ that represents the relatedness between species. In particular, the estimation of the evolution rate $\sigma$ can be performed using comparative phylogenetic methods (PCM) [8,9,10].

There are models created on the basis of Eq. (3) via considering the constant-rates BM model may not be well addressed for the evolution in many scenarios [11,12,13,14,15]. Those models come with the assumption that the evolution of a species changes over time, the rates of evolution can be modeled as either constants or stochastic variables along times (branch lengths). The trait variable $y_t$ adopts the following dynamics.

$$y_t=\mu +{\sigma }_t{W}_t,$$

(4)

where ${\sigma }_t$ can be constant (i.e. ${\sigma }_t=\sigma ,t>0$) [16], piecewise constant (i.e. ${\sigma }_t={\sigma }_i,{\gamma }_{i-1}\le t\le {\gamma }_i$ where ${\gamma }_i$ and ${\gamma }_{i-1}$ are the successive time regimes) [11] or a random variable modeled by another pertinent process (i.e. ${\sigma }_t\sim \mathcal{D}\left(t\right)$ where $\mathcal{D}\left(t\right)$ is a distribution function for the stochastic variable) [17,18].

A hypothetical tree of four species $A,B,C,D$ and the corresponding simulation of the trajectories $y_{i,t}$s using Brownian motion along a rooted phylogenetic tree under two different rates is shown in Figure 2.

Those models have been broadly applied in many studies. For example, in the evolution of the morphology of the world’s largest flowers (Raffkesianeae: up to 1 meters in diameter), [19] found that the enormous flowers evolved from ancestors with tiny flowers. In the study of the evolution of the size of the plant genome, [20] found that the woody lineages had a stochastic motion rate that was nearly five times slower than the rate of the herbaceous lineages. Although the existing framework has produced rate models, none consider a scenario in which the rate is treated as a time-correlated stochastic variable, which could potentially enhance the study of rate evolution. This highlights the need for our work to apply a time series model [21,22]. Specifically, this approach aims to answer key questions: Are the evolutionary rates of biological traits statistically independent, or are they believed to be phylogenetically serially autocorrelated ?

Note that in modeling the rate of evolution, it is essential to consider its dynamic nature, as the rate can fluctuate, increase, or decrease over time [18,23], rather than a constant. Consider that the implementation of the time-correlated rate evolution ${\sigma }_t$ could possibly provide an alternative to reveal embedded information about species evolution, in this work we intend to expand the model $y_t=\mu +{\sigma }_t{W}_t$ in Eq. (4) within the framework of correlated rate evolution (${\sigma }_t=g\left({\sigma }_s\right|\mathsf{\Theta })$ for $s<t$ where $\mathsf{\Theta }$ is a parameter vector) to model the trait evolution for phylogenetic comparative analysis. In particular, we use the autoregressive moving average (ARMA) time series model that has been widely applied in econometrics to model the rate parameter [24,25]. The description of the methods can be found in Section 2. The simulations are detailed in Section 3. Empirical analyzes are presented in Section 4. The discussions and conclusions are covered in Section 5.

2. Methods

2.1. Trait Evolution with Time Correlated Rate

We describe our procedure as follows. In Section 2.1.1, we adopt phylogenetic ridge regression [26] to perform the analysis using the phylogenetic tree and trait data as input to obtain rate estimates. In Section 2.1.2, we perform a phylogenetic rate ARMA(p,q) regression, w here the rate estimates are incorporated along the tree using a traversal algorithm, treating the data as time series.

2.1.1. Step 1: Phylogenetic Ridge Regression

Given the tree with the topology and branch lengths, one assumed that the observed trait values are a linear combination of past values, the time elapsed, and the rate of evolutionary change. To illustrate, we use a case of 3 taxa tree as shown in Figure 3 where the extant species A, B, and C, and the fossil species D, E, and O have trait values of $(y_A,y_B,y_C)$ and $(y_D,y_E,y_0)$. These traits values can be expressed as a linear combination of the time elapsed and the rate of evolution shown in Eq. (5).

$$\begin{array}{cc}\hfill y_A& =y_0+l_A{\sigma }_A+l_D{\sigma }_D+l_E{\sigma }_E,\hfill \\ \hfill y_B& =y_0+l_B{\sigma }_B+l_D{\sigma }_D+l_E{\sigma }_E,\hfill \\ \hfill y_C& =y_0+l_C{\sigma }_C+l_E{\sigma }_E,\hfill \\ \hfill y_D& =y_0+l_D{\sigma }_D+l_E{\sigma }_E,\hfill \\ \hfill y_E& =y_0+l_E{\sigma }_E.\hfill \end{array}$$

(5)

One can formulate into the system of linear equation in Eq. (6)

$$\left[\begin{array}{c}{y}_A\\ y_B\\ y_C\\ y_D\\ y_E\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_0\\ y_0\\ y_0\\ y_0\end{array}\right]+\left[\begin{array}{ccccc}{l}_A& 0& 0& l_D& l_E\\ 0& l_B& 0& l_D& l_E\\ 0& 0& l_C& l_D& l_E\\ 0& 0& 0& l_D& l_E\\ 0& 0& 0& l_D& l_E\end{array}\right]\left[\begin{array}{c}{\sigma }_A\\ {\sigma }_B\\ {\sigma }_C\\ {\sigma }_D\\ {\sigma }_E\end{array}\right].$$

(6)

Given empirical data where the tip values are $y_A$, $y_B$, and $y_C$, and the ancestral values $y_4$ and $y_5$ may be known through fossil records or remain unknown. Assuming the currently provided data $\mathbf{y}=(y_A,y_B,y_C)$, and the phylogenetic tree $\mathcal{T}$ with a given set of known branch lengths $L=\{l_A,l_B,l_C,l_D,l_E\}$, the tip states $(y_A,y_B,y_C)$ can be written into a matrix form in Eq. (7)

$$\left[\begin{array}{c}{y}_A\\ y_B\\ y_C\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_0\\ y_0\end{array}\right]+\left[\begin{array}{ccccc}{l}_A& 0& 0& l_D& l_E\\ 0& l_B& 0& l_D& l_E\\ 0& 0& l_C& 0& l_E\end{array}\right]\left[\begin{array}{c}{\sigma }_A\\ {\sigma }_B\\ {\sigma }_C\\ {\sigma }_D\\ {\sigma }_E\end{array}\right].$$

(7)

Use the can written into Equation (8)

$$\widehat{Y}=L\mathsf{\Sigma },$$

(8)

where

$$\widehat{\mathsf{\Sigma }}=\underset{\mathsf{\Sigma }}{argmin}\sum _i{(y_i-{\left(L\mathsf{\Sigma }\right)}_i)}^2+\lambda \sum _j{\mathsf{\Sigma }}_j^2,$$

(9)

which can be written as

$$\widehat{\mathsf{\Sigma }}={(L^{T}L+\lambda I)}^{-1}{L}^{T}y.$$

(10)

In general, given tip data $\mathbf{y}=\{y_1,y_2,\cdots ,y_n\}\in {\mathcal{R}}^n$, and tree topology $\mathcal{T}$ of n taxa with branch length set $\mathcal{B}=\{l_{j_1},l_{j_2},\cdots ,l_{j_m}\}$ where $(j_1,j_2,\cdots ,j_m):=\mathcal{J}$ is internal node index set with the corresponding with the ancestral-descedant relationship delineated by the tree topology, we use the R package rrphylo [27] to perform the analysis for estimation of phenotypic evolutionary rates ${\left\{{\sigma }_{j_s}\right\}}_{j_s\in \mathcal{J}}$ using Eq. (10) with adopting the leave-one-out cross-validation (LOOCV) to search for the optimal $\lambda$. In this method, each observation is removed one at a time, and the remaining data is used to make predictions. The error for the excluded observation is calculated, and this process is repeated for all observations. The total LOOCV error is computed for various $\lambda$ values, and the $\lambda$ with the smallest error is selected as the best regularization parameter, improving model performance [26].

2.1.2. Step 2: Phylogenetic Rate ARMA(p,q) for Trait Evolution

Next, we implements the ARMA model for studying rate ${\left\{{\sigma }_{j_s}\right\}}_{j_s\in \mathcal{J}}$. Given tree topology $\mathcal{T}$ of n taxa and with tip node $\mathbf{y}=\{y_1,y_2,\cdots ,y_n\}\in {\mathcal{R}}^n$ and internal node index set $\mathcal{J}=(j_1,j_2,\cdots ,j_m)$ and the corresponding branch length set $\mathcal{B}=\{l_{j_1},l_{j_2},\cdots ,l_{j_m}\}$ with the ancestral-descedant relationship delineated by the tree topology, by utilizing phylogenetic ridge regression in Eq. (10) one can get a set of estimated rates ${\mathsf{\Sigma }}_{\mathcal{J}}={\left\{{\sigma }_{j_s}\right\}}_{s=1}^m$ where index $j_s$ corresponds to the internal index, then the phylogenetic ARMA($p,q$) model of rates defined by a mixture of phylogenetic autoregressive (AR) and phylogenetic moving average (MA) models given by the Eq. (11):

$${\sigma }_{j_s}=\sum _{i=1}^p{\varphi }_i{\sigma }_{j_{s-i}}+\sum _{j=1}^q{\theta }_j{ϵ}_{j_{s-j}}+{ϵ}_{j_s},$$

(11)

where ${\sigma }_{j_s}$ represents the time series at time $j_s$, ${\varphi }_i$ are coefficients for the AR terms, ${\theta }_j$ are coefficients for the MA terms, and ${ϵ}_{j_s}$ denotes the white noise error terms following a normal distribution with mean zero and variable proportional to the branch length $j_s$ as shown in Eq. (12)

$${ϵ}_{j_s}\sim \mathcal{N}(0,{\tau }^2{l}_{j_s}).$$

(12)

For a phylogenetic ARMA(1,1) rate model, denote ${\sigma }_{\mathrm{des}}={\sigma }_{j_s}$ as the descendent (${\sigma }_{j_s}$) and ${\sigma }_{\mathrm{anc}}={\sigma }_{j_{s-1}}$ as the ancestor, the model in Eq. (11) can be represented as in Eq. (13).

$${\sigma }_{\mathrm{des}}=\varphi {\sigma }_{\mathrm{anc}}+\theta {ϵ}_{\mathrm{anc}}+{ϵ}_{\mathrm{des}},\mathrm{where}{ϵ}_{\mathrm{des}}\sim \mathcal{N}(0,{\tau }^2{l}_{\mathrm{des}}).$$

(13)

An exemplary example of the 3 taxa case is shown in Figure 4.

We are interested in the joint distribution $L\left(\mathsf{\Theta }\right|\mathbf{y},\mathbf{\sigma },\mathbf{ϵ},\mathcal{T})$ where $\mathsf{\Theta }$ is the parameter vector (i.e. $\mathsf{\Theta }=({\left\{{\varphi }_i\right\}}_{i=1}^p,{\left\{{\theta }_j\right\}}_{j=1}^q,\tau )$ for the ARMA $(p,q)$. On a specific branch $s\in \mathcal{J}$ with branch length $l_s$, we assume the error ${ϵ}_{j_s}$ follow a normal distribution with mean conditioning on the ancestor and variance with branch length $l_{j_s}$. The logarithmic likelihood on a branch is a one- dimensional ARMA $(p,q)$ that has likelihood in Eq. (14).

$$logL\left(\mathsf{\Theta }\right|\mathbf{y},{\sigma }_{j_s},{ϵ}_{j_s},\mathcal{T})={\displaystyle \frac{1}{2}}\left(-log\left(2\pi \right)-log\left(\tau \right)-{\displaystyle \frac{{ϵ}_{j_s}}{\tau }}\right).$$

(14)

Assuming branches are independent, the full log likelihood given the tree $\mathcal{T}$ and trait $\mathbf{Y}$, $\mathbf{\sigma }={\left\{{\sigma }_s\right\}}_{s\in \mathcal{J}},\mathbf{ϵ}={\left\{{ϵ}_s\right\}}_{s\in \mathcal{J}}$ can be written as

$$logL\left(\mathsf{\Theta }\right|\mathbf{y},\mathbf{\sigma },\mathbf{ϵ},\mathcal{T})=logL\left(\mathsf{\Theta }\right|{ϵ}_0,{\sigma }_0,\mathcal{T})+\sum _{s\in \mathcal{J}}logL\left(\mathsf{\Theta }\right|\mathbf{y},{\sigma }_s,{ϵ}_s,\mathcal{T}),$$

(15)

where ${ϵ}_0$ is the residual at the root and the product operator $s\in \mathcal{J}$ follows the corresponding tree topology $\mathcal{T}$.

To perform model inference on parameter estimation, given trait vector $\mathbf{y}={(y_1,y_2,\cdots ,y_n)}^t$ on the tip, and a phylogenetic tree $\mathcal{T}$ with known branch length set $\mathcal{J}$ where $l_0=0$ is the branch length of the root node, and each branch $l_{j_s}$. Use y and the tree $\mathcal{T}$ to obtain the rate estimates at the root ${\sigma }_0$ from the Brownian motion model [12]

$${\sigma }_0=\sqrt{{(\mathbf{y}-\widehat{\mu }\mathbf{1})}^t{\mathbf{C}}^{-1}(\mathbf{y}-\widehat{\mu }\mathbf{1})/n},$$

(16)

where $\widehat{\mu }=1{\mathbf{C}}^{-1}y)/(1^t{\mathbf{C}}^{-1}1)$ is the MLE estimate for the ancestral at the root under the Brownian motion model, and $C$ is the variance covariance matrix transformed from the tree $\mathcal{T}$ and 1 = (1, 1, ⋯, 1)^t is vectors of 1s. Then we use the maximum likelihood approach for parameter estimation and inference. Refer to Section 6.2 Appendix for expression and parameter estimation a first few order phylorateARMA((p, q)) model.

For empirical analysis, we analyze the data using the model up to order $p=2,q=2$. For model selection, the Akaike information criterion $AIC=2k-logL$ and its sample size correction version $AICc=AIC+2k(k+1)/(n-k-1)$ [28] and the weighted parameter are calculated using equation $w_k=exp\left(-0.5\times \Delta {\mathrm{AIC}}_k\right)/{\sum }_{k=1}^{K}exp\left(-0.5\times \Delta {\mathrm{AIC}}_k\right)$, where $w_i$ is the weight for the models i and $\Delta {\mathrm{AIC}}_i={\mathrm{AIC}}_i-{min}_i\left\{{\mathrm{AIC}}_i\right\}$. The weighted parameter estimate $\tilde{\xi }$ is calculated as $\tilde{\xi }={\sum }_{k=1}^K{w}_i\times {\widehat{\xi }}_k$. Here, ${\widehat{\xi }}_i$ represents the MLE estimate from model i. This method incorporates uncertainty in the selection of the model by averaging parameters based on AIC weights [29].

2.2. Testing of AR(1) Effect

We conduct inference to test for the existence of the autoregressive effect. This is equivalent to testing whether $\varphi =0$, as shown in Section 2.2.1. Additionally, we investigate whether the autoregressive effect varies across different clades of the tree. For this, we pose the null hypothesis $H_0:{\varphi }_1={\varphi }_2$, which is discussed in Section 2.2.2.

2.2.1. Test of Correlated Rate Evolution

Given the AR(1) process: ${\sigma }_{j_s}=\varphi {\sigma }_{j_s-1}+{ϵ}_{j_s}$ where ${ϵ}_{j_s}\sim \mathcal{N}(0,{\tau }^2{l}_{j_s})$ is white noise, the asymptotic variance of the MLE $\widehat{\varphi }$ is: $\mathrm{Var}\left(\widehat{\varphi }\right)={\tau }^2/{\sum }_{j_s\in \mathcal{J}}{\sigma }_{j_s}^2$. The null hypothesis is given by: $H_0:\varphi ={\varphi }_0$ where ${\varphi }_0$ is a hypothesized value. The alternative hypothesis for a two-sided test is: $H_1:\varphi \ne {\varphi }_0$. We set ${\varphi }_0=0$ to test for the AR(1) effect. To perform the test, first fit an AR(1) model to the data and obtain the estimate $\widehat{\varphi }$ of $\varphi$ and its standard error $\mathrm{SE}\left(\widehat{\varphi }\right)$. The test statistic is calculated as:

$$t={\displaystyle \frac{\widehat{\varphi }-{\varphi }_0}{\mathrm{SE}\left(\widehat{\varphi }\right)}}.$$

(17)

This test statistic is compared with the critical value of the t distribution with $n-1$ degrees of freedom, where n is the sample size. The standard deviation estimate can be derived from the Hessian matrix obtained during the optimization process, which is then used to conduct hypothesis testing.

In the context of hypothesis testing for phylogenetic rate model in an AR(1) model, the process begins by establishing the null hypothesis, $H_0:\varphi ={\varphi }_0$, where ${\varphi }_0$ represents a specific hypothesized value. The alternative hypothesis for a two-sided test is $H_1:\varphi \ne {\varphi }_0$. The methodology for testing these hypotheses involves several steps: (i) Firstly, an AR(1) model is fitted to the data. From this model, the estimate $\widehat{\varphi }$ of $\varphi$ and its standard error $SE\left(\widehat{\varphi }\right)$ are obtained. (ii) The next step involves calculating the test statistic using the formula $t={\displaystyle \frac{\widehat{\varphi }-{\varphi }_0}{SE\left(\widehat{\varphi }\right)}}$. (iii) Finally, this test statistic is compared against the critical value from the t-distribution, which is determined by the sample size $J=\left|\mathcal{J}\right|$ where $\left|\mathcal{J}\right|$ is the number of the branches of the tree and corresponds to $J-1$ degrees of freedom. This comprehensive method is essential to determine whether the hypothesized value ${\varphi }_0$ is a plausible value for $\varphi$ in the AR(1) model.

2.2.2. Testing Heterogeneity Rates on Sub Clade

We consider whether heterogeneity autoregressive effect exists on the two subclade, as shown in Figure 5.

The null Hypothesis ($H_0$) poses that a single tree with n taxa and use one parameter $\varphi$ while the alternative Hypothesis ($H_a$) states that 2 autoregressive parameters ${\varphi }_1$ and ${\varphi }_2$ would more approporate when tree the is divided into two subtrees from the root, and parameters are estimated independently.

The null hypothesis in the likelihood ratio test states that the simpler model (in this case, the single combined tree) is sufficient to explain the variability in the data. The alternative hypothesis suggests that dividing the tree into two subtrees provides a better fit. The likelihood ratio test (LRT) statistic is computed as:

$$\mathrm{LRT}=-2\left(log{L}_{H_0}-log{L}_{H_a}\right),$$

(18)

where $log{L}_{H_0}$ is the log-likelihood under the null hypothesis and $log{L}_{H_a}$ is the log-likelihood under the alternative hypothesis. The test statistic follows a chi-square distribution with the degrees of freedom $df=1$.

The procedure in our framework is summarized in Algorithm 1.

Algorithm 1: Procedure of Phylogenetic ARMA(p,q) Rate Data Analysis

1: Input: a phylogenetic tree $\mathcal{T}$ with m branch lengthes $\mathcal{B}=\{b_{j_1},b_{j_2},\cdots ,b_{j_m}\}$, and quantitative trait values $\mathbf{y}={(y_1,y_2,\cdots ,y_n)}^t$ of n taxa. 2: Output: estimated rate of evolution ${\mathsf{\Sigma }}_{\mathcal{J}}={\left\{{\tilde{\sigma }}_{j_s}\right\}}_{s=1}^m$, the parameter estimates ( $\widehat{\varphi }s$, $\widehat{\theta }$s, ${\widehat{\tau }}^2$), the AICc values, the test statistics t-value for $H_0:\varphi =0$, and ${\chi }_2$-value for $H_0:{\varphi }_1={\varphi }_2$ and their corresponding p values. 3: Apply the phylogenetic ridge regression using the tree $\mathcal{T}$, branch lengths $\mathcal{B}$, and trait value $\mathbf{Y}$ using Eq. (10) to obtain the rate estimates ${\left\{{\tilde{\sigma }}_{j_s}\right\}}_{s=1}^m$. 4: Use ${\left\{{\tilde{\sigma }}_{j_s}\right\}}_{s=1}^m$ as input and apply the tree traveral algorithm on the tree $\mathcal{T}$ to fit the phylonetic rate ARMA(p,q) model. 5: Evaluate the AICc to compare the phyrateARMA($p,q$) model, and compute the Akaike weight and the weighted estimate. 6: Test of the significance of the AR(1) relationship in the two scenarios: Test of correlated rate evolution $H_0:\varphi =0$. Testing heterogeneity of AR (1) on two subclades $H_0:{\varphi }_1={\varphi }_2$. 7: Return the results.

3. Simulation

We assess the performance of the model using simulation. Four types of trees are used for the assessment: star tree, balanced tree, left tree, coalescent tree, and birth-and-death tree, each of size $16,32,64,128$. The initial samples are drawn from an independent normal distribution. The initial estimate for ${\sigma }_0$ by applying the formula in Eq. (16) at the root is estimated by the Brownian motion model using the R package geiger. For each type of tree, for each data set and for each prior set, a million replicates of the sample will be generated to assess the posterior sample distribution.

Let the tree $\mathcal{T}$ have m branches, with each branch representing a transition between an ancestral node $j_s-1$ and a descendant node $j_s$. The tree is traversed in postorder (from root to tips), and the trait values for each node are computed based on the following AR(1) process.

For each branch, the evolutionary rate at the descendant node $j_s$ is determined by the rate at its ancestor $j_s-1$, with some noise (innovation) added. The autoregressive equation is shown in Eq. (19)

$${\sigma }_{j_s}=\varphi {\sigma }_{j_s-1}+{ϵ}_{j_s}$$

(19)

Initially, the rate value at the root node is set to zero: ${\beta }_0=0$.

After simulating the rates using Eq. (19) for all nodes in the tree, the final trait values at the tree tips (taxa) are computed by accumulating the rates along the path from the root to each tip.

Let $\mathcal{J}$ denote the set of nodes from the root to tip j, and ${\ell }_{j_s}$ denote the branch length leading to node $j_s$. The trait value at each tip j is:

$$y_j=\sum _{j_s\in \mathcal{J}}{\sigma }_{j_s}{\ell }_{j_s}$$

(20)

where $y_j$ is the final accumulated trait value at tip j.

The final output consists of two arrays: (1) the array of simulated evolutionary rates for each node, $\mathbf{\sigma }$, and (2) the trait values at the tree tips, $\mathbf{y}$ summarized by the following equations Eq. (21)

$$\begin{array}{cc}\hfill {\sigma }_{j_s}=\varphi {\sigma }_{j_s-1}+{ϵ}_{j_s},& ϵ\sim \mathcal{N}(0,{\tau }^2{l}_{j_s}),\hfill \\ y_j=\sum _{j_s\in \mathcal{J}}{\sigma }_{j_s}{\ell }_{j_s}.& \end{array}$$

(21)

This simulation captures the evolutionary rate changes along a phylogeny under an AR(1) model, accounting for both ancestral states and random innovation at each branch. We create a function simulate_trait_data_given_phi which simulates the trait evolution process along a phylogenetic tree using an autoregressive process (AR(1)) for the rates of evolution. The simulation setup is shown in Algorithm 2

Algorithm 2: Estimating Type I Error and Power with Phylogenetic AR(1) rate Model

1: Input: tree type {balanced, coalescent, pure birth, birth death }, taxa size $\{16,32,64,128\}$, simulations 1000, $\alpha =0.05$, ${\varphi }_0=0$, ${\sigma }_0=1$, ${\varphi }_a=\{0.05,0.10,0.15,\cdots .0.95,1\}$. 2: Output: Type I error and power for each taxon size. 3: for each tree type do 4:     Generate tree $\mathcal{T}$ with branch length set $\mathcal{B}={\left\{l_{j_s}\right\}}_{j_s\in \mathcal{J}}$ 5:     for each taxa do 6:         Initialize Type I Error $T_1=0$ ▹Estimate Type I Error 7:         for $i=1\mathrm{to}\mathrm{sims}$ do 8:            Use Eq. (21) with $\varphi ={\varphi }_0$ to simulate trait data y 9:            Use Eq. (10) to obtain rate estimates ${\left\{{\tilde{\sigma }}_{j_s}\right\}}_{j_s\in \mathcal{J}}$ 10:            Optimize the neg log likelihood in Eq. (15) to obtain MLE parameters $\widehat{\varphi },\widehat{\theta },\widehat{\tau }$. 11:            Use Eq. (17) to test $\varphi =0$ using MLE parameter, and its standard error and p value. 12:            if $p<\alpha$ then 13:                $T_1=T_1+1$ 14:            end if 15:         end for 16:         Store $T_1$ for this tree and taxa size ▹Estimate Power 17:         for each ${\varphi }_a$ do 18:            Initialize power pow $=0$ 19:            for $i=1\mathrm{to}\mathrm{sims}$ do 20:                Use Eq. (21) with $\varphi ={\varphi }_a$ to simulate trait data y 21:                Use Eq. (10) to obtain rate estimates ${\left\{{\tilde{\sigma }}_{j_s}\right\}}_{j_s\in \mathcal{J}}$ 22:                Optimize the neg log likelihood in Eq. (15) to obtain the MLE parameters $\widehat{\varphi },\widehat{\theta },\widehat{\tau }$. 23:                Use Eq. (17) to test $\varphi =0$ using the MLE parameter and its standard error and p value. 24:                if $p<\alpha$ then 25:                    pow= pow $+1$ 26:                end if 27:            end for 28:            Store pow for this tree and taxa size 29:         end for 30:     end for 31: end for 32: Return$T_1,pow$

Figure 6 and Table 1 show how the statistical power to detect an effect changes with different values of ${\varphi }_a$ and different taxa levels.

From Figure 6 and Table 1, one can envision that the power of the test changes depending on the number of taxa. Generally, as the number of taxa increases, the power also increases. Next, for all taxa levels, as ${\varphi }_a$ increases, the power of the test also increases. This could mean that as the parameter ${\varphi }_a$ increases, the ability to detect an effect or a difference becomes stronger. Third, all curves converge or come closer together at higher values of ${\varphi }_a$, indicating that the power could be more consistent across different levels of taxa at these higher values. This simulation provide the evidence that our methods fit the statistical properties as expect in the performance of evaluating the statistical power and type I error.

4. Empirical Analysis

4.1. Genome Size Rate Evolution in Herbaticus and Woody Plant

This study collected and analyzed data sets on the genome size in pg of herbaticus and woody plant. The primary data source is [30], we includes 35 herbaticus and 33 woody species in this study. The main objective is to test whether there is a time-series relationship in the evolutionary rate of these plant’ genome size.

Figure 7 shows the phylgoenetic tree for herbaticus, woody and their joint tree.

4.1.1. Genome Size: Combined Herbaticus and Woody Species

Table 2 presents the parameter estimates for various autoregressive (AR) and autoregressive moving average (ARMA) models. ARMA(1,1) has a ${\varphi }_1$ value of $-0.600$ and a ${\theta }_1$ value of $0.738$, indicating a moderate correlation between consecutive rates of evolution, with a residual variance of ${\tau }^2=4.663$. The weighted estimates combine values across the models, with ${\varphi }_1=-0.050$ and ${\theta }_1=0.202$, and ${\tau }^2=4.722$, providing a comprehensive view of the underlying evolutionary process.

Using the phylogenetic ARMA(1,1) rate model described in Eq. (13), the equation for ${\sigma }_{\mathrm{des}}$ demonstrates how the descendant’s rate of evolution is influenced by multiple factors. The term $\varphi {\sigma }_{\mathrm{anc}}$ indicates that -5.0% of the ancestor’s evolutionary rate (${\sigma }_{\mathrm{anc}}$) is carried over to the descendant, reflecting a weak and inverse correlation between the two rates. The second term, $\theta {ϵ}_{\mathrm{anc}}$, with a coefficient of $0.202$, adjusts for the random variation (or noise) in the ancestor’s rate. This positive value suggests that random fluctuations in the ancestor’s rate have a small reinforcing effect on the descendant. Finally, ${ϵ}_{\mathrm{des}}$, with variance proportional to ${\tau }^2=4.722$, represents the new random variation specific to the descendant. The weighted estimates offer a balanced perspective on how both inherited traits and random fluctuations shape evolutionary changes across generations.

In Table 3, AR(2) ranks first with the lowest AICc value of $537.870$ and the highest Akaike weight ($w=0.530$), indicating the best fit among the models. AR(1) follows with an AICc of $539.850$ and a weight of $0.2$, suggesting a close fit but slightly less optimal than AR(2). ARMA(1,1) ranks third with an AICc of $540.070$ and a weight of $0.180$. Models ARMA(2,1) and ARMA(2,2) show considerably lower Akaike weights, indicating that increasing model complexity beyond AR(2) does not significantly improve the fit.

4.1.2. Genome Size: Herbaticus Species

The parameter estimates in Table 4 for the herbaticus models show that AR(1) has a ${\varphi }_1$ value of $0.470$ with residual variance ${\tau }^2=0.009$, indicating a moderate influence from the previous rate of evolution. AR(2) introduces a second autoregressive term (${\varphi }_2=-0.019$) and maintains a similar variance of ${\tau }^2=0.009$. The ARMA(1,1) model has ${\varphi }_1=0.491$ and ${\theta }_1=-0.024$, with the same variance of ${\tau }^2=0.009$. More complex models, such as ARMA(2,1) and ARMA(2,2), show mixed ${\varphi }_1$ and ${\theta }_1$ values, but they maintain the same low variance of ${\tau }^2=0.009$.

The weighted estimate balances these models, showing ${\varphi }_1=0.427$, ${\theta }_1=0.021$, and ${\tau }^2=0.009$. This suggests that the descendant’s rate of evolution is influenced by $42.7$% of the ancestor’s rate, with a small adjustment of $2.1$% due to noise from the ancestor, all within a low variance process.

In Table 5, AR(1) ranks first with the lowest AICc value of $-110.550$ and the highest Akaike weight ($w=0.610$), making it the best-fitting model for the herbaticus data. ARMA(1,1) follows in second place with an AICc of $-108.550$ and a weight of $0.220$, while AR(2) ranks fourth with an AICc of $-105.690$ and a weight of $0.050$. More complex models like ARMA(2,1) and ARMA(2,2) have even lower Akaike weights, indicating that they provide little additional explanatory power compared to simpler models.

4.1.3. Genome Size: Woody Species

The parameter estimates in Table 6 for the woody models indicate that AR(1) has a ${\varphi }_1$ value of $-0.013$ and a residual variance of ${\tau }^2=0.002$, suggesting a very weak negative influence from the previous rate of evolution. AR(2) introduces a second autoregressive term (${\varphi }_2=-0.302$) and slightly increases the variance to ${\tau }^2=0.003$. ARMA(1,1) shows stronger dependence with ${\varphi }_1=0.359$ and ${\theta }_1=-0.446$, while maintaining a variance of ${\tau }^2=0.002$.

The weighted estimate balances the models, with ${\varphi }_1=0.023$, ${\theta }_1=-0.044$, and ${\tau }^2=0.002$. This suggests that the descendant’s rate is slightly positively influenced by the ancestor’s rate ($2.3$%) but slightly negatively impacted by the noise ($-4.4$%), with very low overall variance ($0.002$), indicating minimal change between generations.

In Table 7, AR(1) ranks first with the lowest AICc value of $-192.770$ and the highest Akaike weight ($w=0.860$), indicating it is the best-fitting model for the woody data. ARMA(1,1) ranks second with an AICc of $-188.710$ and a much smaller weight ($w=0.110$), while AR(2) ranks last with an AICc of $-180.880$ and no Akaike weight ($w=0.000$). More complex models, such as ARMA(2,1) and ARMA(2,2), also have minimal Akaike weights ($w=0.010$), suggesting they provide little additional explanatory power and overcomplicate the fit for the woody data.

From above one may foreseen that the woody primate and the herbaticus genome size possess different autoregressive effect for the rates of the evolution, we provide the test in the following sections.

4.1.4. Testing Autocorrelation Rate

Interpretation of AR(1) Test Results

Table 8 show the AR1 test result

Table 8 presents the results of the AR(1) test under the null hypothesis $H_0:\varphi =0$ (no autocorrelation) and the alternative hypothesis $H_a:\varphi \ne 0$. The herbaticus tree shows a significant positive autocorrelation with a large z-value of 5.100 and a p-value of 0.000, strongly rejecting the null hypothesis. This suggests a significant effect of the previous rate of evolution on the current rate in the herbaticus dataset.

In contrast, the woody tree show no significant autocorrelation, with z-values of $-0.121$, and p-value $0.903$ indicates that the null hypothesis cannot be rejected for these datasets, suggesting no significant relationship between the past and current rates of evolution in the woody trees. This may due to woody tree with slower rates (5 times slower than the herbaticus tree [20]).

The combined tree shows no significant autocorrelation, with z-values of $1.119$, and p-values of $0.263$. The result suggests that there is no significant relationship between the past and current rates of evolution in the combined trees.

Testing Heterogeneity Rates on subclades

Table 9 presents the results of the heterogeneity test in evolutionary rates across the herbaticus, woody, and combined (dino) subclades using various autoregressive models (AR(1), AR(2), ARMA(1,1), ARMA(2,1), ARMA(2,2)). The log-likelihood values for the herbaticus ($log{L}_{\mathrm{di}}$), woody ($log{L}_{\mathrm{no}}$), and combined ($log{L}_{\mathrm{dino}}$) datasets are provided, along with the chi-squared statistic and p-values.

All models yield highly significant results (p-values = 0.00), indicating strong evidence for rate heterogeneity between the subclades. The highest chi-squared statistic is observed for AR(1) (847.16), suggesting that it explains the greatest amount of variance in the data compared to the other models.

4.2. Body Mass Rate Evolution in Diurnal and Nocturnal Primate

This study collected and analyzed data sets on the body weights of male and female primates. The primary data source are [33] and Animal Diversity Web [34], which includes diurnal and nocturnal male and female primates. The main objective is to test whether there is a time-series relationship in the evolutionary rate of these primates’ body weights. The original data set contains 88 primate species, 34 diurnal species and 54 nocturnal species. Body weight for males, females, or combined is collected. The analysis is as follows, combined analysis of male and female primate data. For the male-only and female-only analysis case, see Section 6.1.2 for the corresponding analysis.

4.2.1. Body Mass: Combined Diurnal and Nocturnal Species

Table 10 presents the parameter estimates for various autoregressive moving average (ARMA) and autoregressive (AR) models. ARMA(1,1) shows the highest ${\varphi }_1$ value of $0.875$ and ${\theta }_1$ of $-0.652$, indicating a significant correlation between consecutive rates of evolution, with residual variance ${\tau }^2=1.103$. The weighted estimates combine the values across the models, with ${\varphi }_1=0.779$ and ${\theta }_1=-0.566$, and $\tau =1.109$ providing a balanced view of the underlying process.

By the phylgoenetic ARMA (1,1) rate model in Eq. (13), as the equation for ${\sigma }_{\mathrm{des}}$ shows how the descendant’s rate of evolution is determined by a combination of factors. The term $\varphi {\sigma }_{\mathrm{anc}}$ suggests that 77.9% of the ancestor’s evolutionary rate (${\sigma }_{\mathrm{anc}}$) carries over to the descendant, indicating a strong correlation between the two rates. The second term, $\theta {ϵ}_{\mathrm{anc}}$, which has a coefficient of $-0.566$, adjusts for the random variation (or noise) in the ancestor’s rate. This negative value means that random fluctuations in the ancestor’s rate are counterbalanced, effectively reducing the influence of ancestor-specific randomness on the descendant. Finally, ${ϵ}_{\mathrm{des}}$ has variance proportional to ${\tau }^2=1.109$ representing the new random variation specific to the descendant. The weighted estimates provide a balanced view of the process by combining these terms, giving insight into how both inherited traits and random fluctuations influence evolutionary changes across generations.

In Table 11, ARMA(1,1) ranks first with the lowest AICc value of $417.080$ and the highest Akaike weight ($w=0.610$), indicating the best fit among the models. ARMA(2,1) and ARMA(2,2) follow, but with considerably lower Akaike weights, suggesting that increasing model complexity beyond ARMA (1,1) does not produce substantial fit improvements. AR(1) is the lowest, highlighting its reduced explanatory power compared to more complex models.

4.2.2. Body Mass: Diurnal Species

The parameter estimates in Table 10 for the diurnal models show that AR(1) has a ${\varphi }_1$ value of $0.422$ with residual variance ${\tau }^2=0.129$, indicating a moderate level of influence from the previous rate of evolution. AR(2) introduces a second autoregressive term (${\varphi }_2=-0.206$) and slightly increases the variance to ${\tau }^2=0.131$. ARMA(1,1) and more complex models, such as ARMA(2,1) and ARMA(2,2), demonstrate negative ${\varphi }_1$ and ${\theta }_1$ values, but these more complex models have similar or slightly lower variances (${\tau }^2=0.127$). The weighted estimate balances between models, showing ${\varphi }_1=0.260$, ${\theta }_1=0.222$, and ${\tau }^2=0.127$. The weighted estimate shows that the descendant’s rate is influenced by $26.0$% of the ancestor’s rate, with a positive adjustment of $22.2$% from the ancestor’s noise, while the process has a relatively low variance of $0.127$.

In Table 13, AR (1) is the first with the lowest AICc value of $45.680$ and the highest Akaike weight ($w=0.480$), making it the best-fitting model for the diurnal data. ARMA(1,1) ranks second with an AICc of $47.170$ and a lower weight ($w=0.230$), followed by AR(2) with an AICc of $47.710$ ($w=0.170$). More complex models such as ARMA(2,1) and ARMA(2,2) rank lowest, with minimal Akaike weights, indicating that they provide little additional explanatory power.

4.2.3. Body Mass: Nocturnal Species

The parameter estimates in Table 14 for the nocturnal models indicate that AR(1) has a ${\varphi }_1$ value of $-0.213$ and a residual variance of ${\tau }^2=0.003$, suggesting a weak negative influence from the previous rate of evolution. AR(2) introduces a second autoregressive term (${\varphi }_2=-0.080$) and slightly increases the variance to ${\tau }^2=0.004$. ARMA (1,1) has larger ${\varphi }_1=0.593$ and ${\theta }_1=-0.901$, with a similar variance of ${\tau }^2=0.004$, indicating a stronger dependency on the autoregressive and moving average terms. The weighted estimate shows a balance between models with ${\varphi }_1=-0.181$, ${\theta }_1=-0.037$, and ${\tau }^2=0.003$. The weighted estimate indicates that the descendant’s rate is slightly negatively influenced by the ancestor’s rate ($-18.1$%) and its noise ($-3.7$%), with very low overall variance ($0.003$), suggesting minimal change between generations.

In Table 15, AR (1) is the first with the lowest AICc value of $-255.450$ and the highest Akaike weight ($w=0.930$), indicating it is the best-fitting model for the nocturnal data. AR(2) and ARMA(1,1) rank much lower, with AICc values of $-248.710$ and $-248.690$, respectively, both having a very small Akaike weight ($w=0.030$). More complex models such as ARMA(2,1) and ARMA(2,2) rank lowest, with minimal Akaike weights, indicating they add little explanatory power and overcomplicate the fit for the nocturnal data.

From above one may foreseen that the nocturnal primate and the diurnal primate body mass possess different autoregressive effect for the rates of the evolution, we provide the test in the following sections.

4.2.4. Testing Autocorrelation Rate

Interpretation of AR(1) Test Results

Table 16 show the AR1 test result

Table 16 presents the results of the AR(1) test under the null hypothesis $H_0:\varphi =0$ (no autocorrelation) and the alternative hypothesis $H_a:\varphi \ne 0$. The combined tree and the diurnal tree show significant positive autocorrelation with large z-values ($4.229$ and $4.237$, respectively) and p values of $0.000$, indicating a strong rejection of the null hypothesis. This suggests a significant effect of the previous evolution rate on the current rate in these datasets.

Interestingly, the nocturnal tree exhibits a significant negative autocorrelation with a z value of $-2.763$ and a p value of $0.006$. Despite the negative z-value, the small p-value suggests a significant relationship. This indicates that while each tree shows a significant effect, the nocturnal tree behaves differently from the combined and diurnal trees, potentially indicating unique evolutionary dynamics in nocturnal species.

Testing Heterogeity Rates on subclades

Table 17 presents the results of the test of heterogeneity in evolutionary rates in the diurnal, nocturnal and combined (dino) subclades using various autoregressive models (AR(1), AR(2), ARMA(1,1), ARMA(2,1), ARMA(2,2)). The log-likelihood values for the diurnal ($logL$_di), nocturnal ($logL$_no), and combined ($logL$_dino) data are shown, along with the chi-squared statistic and p-values. All models show highly significant results (p-values = 0.00), indicating strong evidence for rate heterogeneity between subclades. The chi-squared statistic is highest for AR1 ($637.06$), suggesting it explains the greatest amount of variance in the data compared to the other models.

5. Discussion and Conclusions

In this study, we first assumed that species evolve along a phylogenetic tree, with time length, evolutionary rate, and trait values expressed as a linear combination. Using the given trait values and branch lengths of the phylogenetic tree, we applied ridge regression and used the R package RRphylo to estimate the evolutionary rate. Next, we modeled the evolutionary rate using the ARMA(p,q) models, fitting these models to ancestral trait values, and then applied these models to estimate present-day trait values.

In the empirical analysis, we tested whether the trait evolution of diurnal and nocturnal populations exhibited time-series correlation. In the dataset, which combined male and female primates body mass, the results showed an autoregressive effect in the separate analysis of diurnal and nocturnal populations, indicating a time-series relationship. To observe the evolutionary dynamics, the branches might have different autoregressive coefficients ($\varphi$). These branches can be analyzed to understand how differences between them affect the overall model. Further analysis could deepen the exploration of the effects of this statistical model in the biological context, providing a more comprehensive scientific explanation to more precisely understand and interpret evolutionary processes and phylogenetic patterns in biological systems.

Several methods exist for estimating autoregressive (AR) and moving average (MA) model parameters. Maximum Likelihood Estimation (MLE) maximizes a likelihood function, assuming normally distributed errors, while Least Squares Estimation minimizes squared residuals for AR parameters but is less suited for MA parameters. The Yule-Walker Equations use sample autocorrelations to estimate AR parameters but not MA ones, whereas the Innovations Algorithm estimates both AR and MA parameters using residuals. Burg’s Method improves AR estimates by minimizing forecast errors, and Bayesian Estimation employs prior distributions and observed data, often via MCMC. The Bootstrapping Method starts by estimating the AR(1) parameters ${\varphi }_1$ and ${\varphi }_2$ for both time series. Then, bootstrapping is used to resample with replacement from each series multiple times. For each bootstrapped sample, calculate $\Delta \varphi ={\widehat{\varphi }}_1-{\widehat{\varphi }}_2$. Next, obtain the empirical distribution of $\Delta \varphi$. If 0 falls outside the 95% confidence interval, there is evidence that the two series have different AR(1) parameters. Each method has advantages depending on the time series characteristics, analysis goals, and available resources.

The proposed research on Covid evolution could be extended by investigating how regional factors such as vaccination rates, public health interventions, and population density contribute to the mutation rate of SARS-CoV-2, or by exploring correlations between mutations and viral traits like transmissibility or immune escape. For land mammals, an extension could examine how adaptability in diet breadth, climate niche, and range size has changed over evolutionary periods or incorporate phylogenetic data to study how related species respond to environmental pressures. Lastly, the study on latitude’s influence on elevation could explore how this relationship differs across biomes or is affected by human activities like deforestation and urbanization, providing a more comprehensive understanding of climate-related changes. These extensions would enhance the versatility and application of the phylogenetic autoregressive moving average regression model across various domains.

6. Appendix

6.1. Code and Scripts

Script and data can be accessed through the following link https://tonyjhwueng.info/phyarmarate.

6.1.1. Model Script

• Main empirical analysis: https://tonyjhwueng.info/phyarmarate/mainPhyarima.R
• Model phylorates: https://tonyjhwueng.info/phyarmarate/testphyloratesV3.r

6.1.2. Figures and Table

• Figure 1: https://tonyjhwueng.info/phyarmarate/twosigbmv3.html
• Figure 2: https://tonyjhwueng.info/phyarmarate/4spetreemap.html
• Figure 3: https://tonyjhwueng.info/phyarmarate/yLs.pptx
• Figure 4: https://tonyjhwueng.info/phyarmarate/armarate2.pptx
• Figure 5: https://tonyjhwueng.info/phyarmarate/tworatetree.html
• Figure 6, Table 1: https://tonyjhwueng.info/phyarmarate/onerateT1powerSummarizer.html
• Section 4.1Figure 7Table 2, Table 3, Table 4, Table 5, Table 6, Table 7, Table 8, Table 9: https://tonyjhwueng.info/phyarmarate/plantWoHeGsize.html
• Section 4.2Figure 8, Table 10, Table 11, Table 12, Table 13, Table 14, Table 15, Table 16, Table 17: https://tonyjhwueng.info/phyarmarate/primateDiNoBodyMass.html
• Additional:

  −

  Male primate: https://tonyjhwueng.info/phyarmarate/primateDiNoBodyMassMale.html

  −

  Female primate: https://tonyjhwueng.info/phyarmarate/primateDiNoBodyMassFemale.html

• Section 6.3: https://tonyjhwueng.info/phyarmarate/TestPhyArimaExample.html

6.2. Phylogenetic ARMA(p,q) Model for Rate Evolution

Given tree topology $\mathcal{T}$ of n taxa and with tip node $\mathbf{y}=(y_1,y_2,\cdots ,y_n)$ and internal node $\mathcal{J}=(j_1,j_2,\cdots ,j_m)$ and the corresponding branch length set $\mathcal{B}=\{l_{j_1},l_{j_2},\cdots ,l_{j_m}\}$ with the ancestral relationship, with fitting the RRphylo [27], one can get a set of estimated rates ${\mathsf{\Sigma }}_{\mathcal{J}}={\left\{{\sigma }_{j_s}\right\}}_{s=1}^m$ for the tree, the phylogenetic ARMA($p,q$) model of rates is a mixture of phylogenetic autoregressive (AR) and phylogenetic moving average (MA) models given by the equation:

$${\sigma }_{j_s}=c+\sum _{i=1}^p{\varphi }_i{\sigma }_{j_{s-i}}+\sum _{j=1}^q{\theta }_j{ϵ}_{j_{s-j}}+{ϵ}_{j_s},\mathrm{where}{ϵ}_{j_s}\sim \mathcal{N}(0,{\tau }^2{l}_{j_s})$$

(22)

where ${\sigma }_{j_s}$ represents the time series at time $j_s$, c is a constant, ${\varphi }_i$ are coefficients for the AR terms, ${\theta }_j$ are coefficients for the MA terms, and ${ϵ}_{j_s}$ denotes the white noise error terms with mean 0 and variance ${\tau }^2{l}_{j_s}$.

6.2.1. PhyRateARMA(1,0) Model

The phylogenetic AR(1) process for rate evolution along the is written as

$${\sigma }_{j_s}=\varphi {\sigma }_{j_s^{\prime }}+{ϵ}_{j_s},$$

(23)

where nodes $j_s,j_s^{\prime }\in \mathcal{J}$ and ${ϵ}_{j_s}\sim \mathcal{N}(0,{\tau }^2{l}_{j_s}).$

The negative log-likelihood is:

$$L(\varphi ,{\tau }^2|\mathcal{T},{\mathsf{\Sigma }}_{\mathcal{J}})={\displaystyle \frac{m}{2}}log\left(2\pi \right)+{\displaystyle \frac{m}{2}}log\left({\tau }^2\right)+{\displaystyle \frac{1}{2{\tau }^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }}\in {\mathsf{\Sigma }}_{\mathcal{J}}}{({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }})}^2.$$

(24)

where $m=\left|\mathcal{J}\right|$ is the number of branches, and ${\sigma }_{j_s}$ and ${\sigma }_{j_s^{\prime }}$ are the two successive rates along the branch length $l_{j_s}$ and $l_{j_s^{\prime }}$, respectively.

To find the maximum likelihood estimates of $\varphi$ and ${\tau }^2$ for the phyAR(1) rate process, we aim to minimize the negative log-likelihood. This requires setting the derivatives of the negative log-likelihood with respect to $\varphi$ and ${\tau }^2$ to zero and solving these equations for the parameters.

The first order conditions for the derivatives are: $-L_{\varphi }=0$, and $-L_{{\tau }^2}=0.$ For $\varphi$, the derivative simplifies to:

$$-L_{\varphi }=-{\displaystyle \frac{1}{{\tau }^2}}\sum _{({\sigma }_{j_s},{\sigma }_{j_s^{\prime }})\in {\mathsf{\Sigma }}_{\mathcal{J}}}({\sigma }_{j_s^{\prime }}{\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }}^2)=0,$$

(25)

which leads to:

$$\widehat{\varphi }={\displaystyle \frac{{\sum }_{({\sigma }_{j_s},{\sigma }_{j_s^{\prime }})\in {\mathsf{\Sigma }}_{\mathcal{J}}}{\sigma }_{j_s^{\prime }}{\sigma }_{j_s}}{{\sum }_{({\sigma }_{j_s},{\sigma }_{j_s^{\prime }})\in {\mathsf{\Sigma }}_{\mathcal{J}}}{\sigma }_{j_s^{\prime }}^2}}.$$

(26)

For ${\tau }^2$, the derivative is:

$$-L_{{\tau }^2}={\displaystyle \frac{m}{2{\tau }^2}}-{\displaystyle \frac{1}{2{\left({\tau }^2\right)}^2}}\sum _{({\sigma }_{j_s},{\sigma }_{j_s^{\prime }})\in {\mathsf{\Sigma }}_{\mathcal{J}}}{({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }})}^2=0,$$

(27)

leading to:

$${\widehat{\tau }}^2={\displaystyle \frac{1}{m}}\sum _{({\sigma }_{j_s},{\sigma }_{j_s^{\prime }})\in {\mathsf{\Sigma }}_{\mathcal{J}}}{({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }})}^2.$$

(28)

These formulations allow us to estimate the parameters that best fit our phylogenetic AR(1) model based on the observed evolutionary rates along a phylogenetic tree.

6.2.2. PhyRateARMA(2,0) Model

A phylogenetic autoregressive model of order 2 (Phylo AR(2)) can be expressed for evolutionary rate changes along a phylogenetic tree. The rate at node t is modeled as:

$${\sigma }_{j_s}={\varphi }_1{\sigma }_{j_s^{\prime }}+{\varphi }_2{\sigma }_{j_s^{″}}+{ϵ}_{j_s},$$

(29)

where ${\sigma }_{j_s^{″}}$ is the rate on the next down successive branch.

Given a phylogenetic lineage with rates ${\sigma }_1,{\sigma }_2,\dots ,{\sigma }_T$, the likelihood function for the Phylo AR(2) model parameters ${\varphi }_1,{\varphi }_2$ and ${\tau }^2$ is:

$$L({\varphi }_1,{\varphi }_2,{\tau }^2|\mathcal{T},{\mathsf{\Sigma }}_{\mathcal{J}})=\prod _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }},{\sigma }_{j_s^{″}}\in {\mathsf{\Sigma }}_{\mathcal{J}}}{\displaystyle \frac{1}{\sqrt{2\pi {\tau }^2}}}exp\left(-{\displaystyle \frac{{({\sigma }_{j_s}-{\varphi }_1{\sigma }_{j_s^{\prime }}-{\varphi }_2{\sigma }_{j_s^{″}})}^2}{2{\tau }^2}}\right)$$

(30)

The log-likelihood function simplifies to:

$$logL({\varphi }_1,{\varphi }_2,{\tau }^2|\mathcal{T},{\mathsf{\Sigma }}_{\mathcal{J}})=-{\displaystyle \frac{m-2}{2}}log\left(2\pi {\tau }^2\right)-{\displaystyle \frac{1}{2{\tau }^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }},{\sigma }_{j_s^{″}}\in {\mathsf{\Sigma }}_{\mathcal{J}}}{({\sigma }_{j_s}-{\varphi }_1{\sigma }_{j_s^{\prime }}-{\varphi }_2{\sigma }_{j_s^{″}})}^2$$

(31)

To find the maximum likelihood estimates (MLE) for ${\varphi }_1$, ${\varphi }_2$, and ${\tau }^2$, the log-likelihood function’s partial derivatives are taken with respect to each parameter, set to zero, and solved for these parameters.

The equations for partial derivatives are:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial {\varphi }_1}}L({\varphi }_1,{\varphi }_2,{\tau }^2|\mathcal{T},{\mathsf{\Sigma }}_{\mathcal{J}})& ={\displaystyle \frac{-1}{{\tau }^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }},{\sigma }_{j_s^{″}}\in {\mathsf{\Sigma }}_{\mathcal{J}}}({\sigma }_{j_s}-{\varphi }_1{\sigma }_{j_s^{\prime }}-{\varphi }_2{\sigma }_{j_s^{″}}){\sigma }_{j_s^{\prime }},\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial {\varphi }_2}}logL({\varphi }_1,{\varphi }_2,{\tau }^2|\mathcal{T},{\mathsf{\Sigma }}_{\mathcal{J}})& ={\displaystyle \frac{-1}{{\tau }^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }},{\sigma }_{j_s^{″}}\in {\mathsf{\Sigma }}_{\mathcal{J}}}({\sigma }_{j_s}-{\varphi }_1{\sigma }_{j_s^{\prime }}-{\varphi }_2{\sigma }_{j_s^{″}}){\sigma }_{j_s^{″}},\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial {\tau }^2}}logL({\varphi }_1,{\varphi }_2,{\tau }^2|\mathcal{T},{\mathsf{\Sigma }}_{\mathcal{J}})& =-{\displaystyle \frac{m-2}{2{\tau }^2}}+{\displaystyle \frac{1}{2{\left({\tau }^2\right)}^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }},{\sigma }_{j_s^{″}}\in {\mathsf{\Sigma }}_{\mathcal{J}}}{({\sigma }_{j_s}-{\varphi }_1{\sigma }_{j_s^{\prime }}-{\varphi }_2{\sigma }_{j_s^{″}})}^2.\hfill \end{array}$$

(34)

Solving these equations provides the MLEs for ${\varphi }_1$, ${\varphi }_2$, and ${\tau }^2$ in a phylogenetic context.

6.2.3. PhyRateARMA(1,1) Model

Consider an ARMA(1, 1) model for a time series ${\left\{{\sigma }_{j_s}\right\}}_{j_s}\in \mathcal{J}$, which can be written as:

$${\sigma }_{j_s}={\varphi }_1{\sigma }_{j_s^{\prime }}+{ϵ}_{j_s^{\prime }}+{\theta }_1{ϵ}_{j_s^{″}},$$

(35)

where ${ϵ}_{j_s^{″}}$ represents the corresponding white noise error terms with variance ${\tau }^2{l}_{j_s^{″}}$ associated with the noise of the evolutionary process.

The likelihood function for the ARMA(1, 1) process, assuming that the errors are normally distributed, is given by the product of the probability densities of the innovations ${\eta }_t$, which are the one-step forecast errors:

$$\begin{array}{cc}\hfill L(\varphi ,\theta ,{\tau }^2|\mathcal{T},{\mathsf{\Sigma }}_{\mathcal{J}})& =\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }}\in {\mathsf{\Sigma }}_{\mathcal{J}}}{\displaystyle \frac{1}{\sqrt{2\pi {\tau }^2}}}exp\left(-{\displaystyle \frac{{({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }}-\theta {\epsilon }_{j_s^{\prime }})}^2}{2{\tau }^2}}\right),\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill logL(\varphi ,\theta ,{\tau }^2|\mathcal{T},{\mathsf{\Sigma }}_{\mathcal{J}})& =-{\displaystyle \frac{m-1}{2}}log\left(2\pi \right)-{\displaystyle \frac{m-1}{2}}log\left({\tau }^2\right)\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & -{\displaystyle \frac{1}{2{\tau }^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }}\in {\mathsf{\Sigma }}_{\mathcal{J}}}{({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }}-\theta {\epsilon }_{j_s^{\prime }})}^2,\hfill \end{array}$$

(38)

The MLE estimates for $\varphi$, $\theta$, and ${\tau }^2$ are found by taking the partial derivatives of the log-likelihood function with respect to these parameters and setting them to zero:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial logL}{\partial \varphi }}& =-{\displaystyle \frac{1}{{\tau }^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }}\in {\mathsf{\Sigma }}_{\mathcal{J}}}({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }}-\theta {\epsilon }_{j_s^{\prime }}){\sigma }_{j_s^{\prime }}=0,\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial logL}{\partial \theta }}& =-{\displaystyle \frac{1}{{\tau }^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }}\in {\mathsf{\Sigma }}_{\mathcal{J}}}({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }}-\theta {\epsilon }_{j_s^{\prime }}){\epsilon }_{j_s^{\prime }}=0,\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial logL}{\partial {\tau }^2}}& =-{\displaystyle \frac{m-1}{2{\tau }^2}}+{\displaystyle \frac{1}{2{\left({\tau }^2\right)}^2}}\sum _{{\sigma }_{j_s},{\sigma }_{j_s^{\prime }}\in {\mathsf{\Sigma }}_{\mathcal{J}}}{({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }}-\theta {\epsilon }_{j_s^{\prime }})}^2=0.\hfill \end{array}$$

(41)

These equations generally require numerical optimization to solve because the innovations ${({\sigma }_{j_s}-\varphi {\sigma }_{j_s^{\prime }}-\theta {\epsilon }_{j_s^{\prime }})}^2$ depend on both $\varphi$ and $\theta$ in a non-linear way. To estimate the parameters $\varphi$ and $\theta$, we use iterative numerical optimization methods that handle nonlinearity. Specifically, we apply L-BFGS-B method [36] by initialized with a starting point and boundary constraints, to maximize the likelihood function. This is done by the R function optim.

6.2.4. PhyRateARMA(2,2) Model

The ARMA(2,2) model is a mixture of autoregressive (AR) and moving average (MA) models, used to model and forecast time series data. It is specified as ARMA(p,q) where p is the order of the AR part and q is the order of the MA part.

An phylogenetic ARMA(2,2) rate model is given by the equation:

$${\sigma }_{j_s}={\varphi }_1{\sigma }_{j_s^{\prime }}+{\varphi }_2{\sigma }_{j_{s^{\prime \prime }}}+{\theta }_1{ϵ}_{j_s^{\prime }}+{\theta }_2{ϵ}_{j_{s^{\prime \prime }}}+{ϵ}_{j_s}$$

(42)

where ${\sigma }_{j_s}$ represents the time series at time $j_s$, c is a constant, ${\varphi }_1,{\varphi }_2$ are coefficients for the AR terms, ${\theta }_1,{\theta }_2$ are coefficients for the MA terms, and ${ϵ}_{j_s}$ denotes the white noise error terms.

Similar to ARMA(1,1), to estimate the parameters $({\varphi }_1,{\varphi }_2)$, $({\theta }_1,{\theta }_2)$ and $\tau$, we again apply L-BFGS-B method [36] to maximize the likelihood function done by the R function optim.

6.3. Test Example

Figure 9:

Figure 9. A rooted phylogenetic tree $\mathcal{T}$ of 8 taxa.

Figure 9. A rooted phylogenetic tree $\mathcal{T}$ of 8 taxa.

Trait values in Table 18.

Table 18. Trait value.

${\mathrm{y}}_1$	${\mathrm{y}}_2$	${\mathrm{y}}_3$	${\mathrm{y}}_4$	${\mathrm{y}}_5$	${\mathrm{y}}_6$	${\mathrm{y}}_7$	${\mathrm{y}}_8$
$-0.03$	$-0.09$	$-0.45$	$-0.85$	$-0.51$	$-0.05$	$-1.31$	$-0.57$

Table 18. Trait value.

${\mathrm{y}}_1$	${\mathrm{y}}_2$	${\mathrm{y}}_3$	${\mathrm{y}}_4$	${\mathrm{y}}_5$	${\mathrm{y}}_6$	${\mathrm{y}}_7$	${\mathrm{y}}_8$
$-0.03$	$-0.09$	$-0.45$	$-0.85$	$-0.51$	$-0.05$	$-1.31$	$-0.57$

Table 19. Edge lengths obtained by tree $\mathcal{T}$ with branch lengths using ape package [37].

${\mathrm{eg}}_1$	${\mathrm{eg}}_2$	${\mathrm{eg}}_3$	${\mathrm{eg}}_4$	${\mathrm{eg}}_5$	${\mathrm{eg}}_6$	${\mathrm{eg}}_7$
$0.143$	$0.143$	$0.143$	$0.143$	$0.286$	$0.286$	$0.143$
${\mathrm{eg}}_8$	${\mathrm{eg}}_9$	${\mathrm{eg}}_{10}$	${\mathrm{eg}}_{11}$	${\mathrm{eg}}_{12}$	${\mathrm{eg}}_{13}$	${\mathrm{eg}}_{14}$
$0.143$	$0.143$	$0.143$	$0.143$	$0.286$	$0.286$	$0.571$

Table 19. Edge lengths obtained by tree $\mathcal{T}$ with branch lengths using ape package [37].

${\mathrm{eg}}_1$	${\mathrm{eg}}_2$	${\mathrm{eg}}_3$	${\mathrm{eg}}_4$	${\mathrm{eg}}_5$	${\mathrm{eg}}_6$	${\mathrm{eg}}_7$
$0.143$	$0.143$	$0.143$	$0.143$	$0.286$	$0.286$	$0.143$
${\mathrm{eg}}_8$	${\mathrm{eg}}_9$	${\mathrm{eg}}_{10}$	${\mathrm{eg}}_{11}$	${\mathrm{eg}}_{12}$	${\mathrm{eg}}_{13}$	${\mathrm{eg}}_{14}$
$0.143$	$0.143$	$0.143$	$0.143$	$0.286$	$0.286$	$0.571$

Use (10) (and package RRphylo) to obtain rate estimate $\left\{{\sigma }_s\right\}$ in Table 20.

Table 20. Rate estimates $\left\{{\widehat{\sigma }}_{j_s}\right\}$ from rrphylo.

${\mathrm{nd}}_1$	${\mathrm{nd}}_2$	${\mathrm{nd}}_3$	${\mathrm{nd}}_4$	${\mathrm{nd}}_5$	${\mathrm{nd}}_6$	${\mathrm{nd}}_7$	${\mathrm{nd}}_8$
$-0.320$	$-0.034$	$0.053$	$-0.070$	$-1.247$	$-0.178$	$-0.445$	$0.003$
${\mathrm{nd}}_9$	${\mathrm{nd}}_{10}$	${\mathrm{nd}}_{11}$	${\mathrm{nd}}_{12}$	${\mathrm{nd}}_{13}$	${\mathrm{nd}}_{14}$	${\mathrm{nd}}_{15}$
$0.023$	$-0.041$	$0.006$	$-0.062$	$-0.027$	$-0.054$	$-0.168$

Table 20. Rate estimates $\left\{{\widehat{\sigma }}_{j_s}\right\}$ from rrphylo.

${\mathrm{nd}}_1$	${\mathrm{nd}}_2$	${\mathrm{nd}}_3$	${\mathrm{nd}}_4$	${\mathrm{nd}}_5$	${\mathrm{nd}}_6$	${\mathrm{nd}}_7$	${\mathrm{nd}}_8$
$-0.320$	$-0.034$	$0.053$	$-0.070$	$-1.247$	$-0.178$	$-0.445$	$0.003$
${\mathrm{nd}}_9$	${\mathrm{nd}}_{10}$	${\mathrm{nd}}_{11}$	${\mathrm{nd}}_{12}$	${\mathrm{nd}}_{13}$	${\mathrm{nd}}_{14}$	${\mathrm{nd}}_{15}$
$0.023$	$-0.041$	$0.006$	$-0.062$	$-0.027$	$-0.054$	$-0.168$

Table 21. Maximum likelihood for parameter estimates under each model.

Model	${\mathit{\varphi }}_1$	${\mathit{\varphi }}_2$	${\mathit{\theta }}_1$	${\mathit{\theta }}_2$	${\mathit{\tau }}^2$
AR(1)	$0.296$				$0.130$
AR(2)	$0.303$	$-0.017$			$0.141$
ARMA(1,1)	$0.582$		$-0.396$		$0.125$
ARMA(2,1)	$-0.087$	$0.669$	$-0.396$		$0.125$
ARMA(2,2)	$-0.046$	$0.629$	$-0.120$	$-0.276$	$0.125$
$\tilde{\xi }$	$0.322$	$0.070$	$-0.126$	$-0.008$	$0.130$

Table 21. Maximum likelihood for parameter estimates under each model.

Model	${\mathit{\varphi }}_1$	${\mathit{\varphi }}_2$	${\mathit{\theta }}_1$	${\mathit{\theta }}_2$	${\mathit{\tau }}^2$
AR(1)	$0.296$				$0.130$
AR(2)	$0.303$	$-0.017$			$0.141$
ARMA(1,1)	$0.582$		$-0.396$		$0.125$
ARMA(2,1)	$-0.087$	$0.669$	$-0.396$		$0.125$
ARMA(2,2)	$-0.046$	$0.629$	$-0.120$	$-0.276$	$0.125$
$\tilde{\xi }$	$0.322$	$0.070$	$-0.126$	$-0.008$	$0.130$

Table 22. Model output. k number of parameter, AICc weight w.

Model	k	$log\mathit{L}$	$\mathit{A}\mathit{I}\mathit{C}\mathit{c}$	w	Rank
AR1	2	$8.66$	$-13.33$	$0.60$	1st
ARMA1	3	$8.67$	$-11.33$	$0.22$	2nd
ARMA21	4	$8.67$	$-9.33$	$0.08$	3th
AR2	3	$7.52$	$-9.04$	$0.07$	4th
ARMA22	5	$8.67$	$-7.33$	$0.03$	5th

Table 22. Model output. k number of parameter, AICc weight w.

Model	k	$log\mathit{L}$	$\mathit{A}\mathit{I}\mathit{C}\mathit{c}$	w	Rank
AR1	2	$8.66$	$-13.33$	$0.60$	1st
ARMA1	3	$8.67$	$-11.33$	$0.22$	2nd
ARMA21	4	$8.67$	$-9.33$	$0.08$	3th
AR2	3	$7.52$	$-9.04$	$0.07$	4th
ARMA22	5	$8.67$	$-7.33$	$0.03$	5th

6.4. Trait Dataset for Empirical Analysis

For plant data and primate data, the tree files in newick format as well as the trait dataset in csv files can be accessed at https://tonyjhwueng.info/phyarmarate/dataset