A simple hierarchical model for heterogeneity in the evolutionary correlation on a phylogenetic tree

A hierarchical set of models

One obvious limitation of the approach illustrated in Fig. 1 is that it only considers two possible alternative models for the evolutionary variance-covariance matrix among traits: one in which both the evolutionary variances and the evolutionary correlation are constant; and a second in which the two mapped regimes have no similarity in evolutionary process for the two traits (Fig. 1C). In fact, it is possible to identify a number of different alternative models between these two extremes. We refer to these models as hierarchical in a similar way to how the set of sequence evolution models used for phylogeny estimation is hierarchical (Posada & Crandall, 1998): each model of increasing complexity has another simpler model as a special case.

Table 1 lists a total of eight alternative models (our original two models, from above, and six others). In square parentheses after each model, we have also provided the alphanumeric code that has been used to denominate the different models in the phytools (Revell, 2012) R software function evolvcv.lite where these models are implemented.

The eight models of Table 1 are as follows: model (1) common rates, common correlation; model (2) different rates, common correlation; model (2b) different rates for trait 1 only, common correlation; model (2c) different rates for trait 2 only, common correlation; model (3) common rates, different correlation; model (3b) different rates for trait 1 only, different correlation; model (3c) different rates for trait 2 only, different correlation; finally, model (4) no common structure between the two different evolutionary variance-covariance matrices of the multivariate Brownian process.

When we analyze this complete set of models for our centrarchid dataset of Fig. 1, we find that the best fitting model (that is, the model with the highest log-likelihood) is the no common structure model in which the Brownian evolutionary variance-covariance matrix is free to differ in all possible ways depending on the mapped regime. It is, in fact, a logical necessity that model 4 has a log-likelihood that is greater than or equal to the next best model. This is because model 4, our no common structure model, has all of our other seven models as special cases. On the other hand, the best supported model (that is, the model that’s best-justified by our data taking into account model complexity; Burnham & Anderson, 2002) is model 3c (different rates in trait 2, relative buccal length, and different correlations; Table 1), indicated with bold font in the table.

Note that some other software, such as the mvMORPH R package of Clavel, Escarguel & Merceron (2015) and the ratematrix package of Caetano & Harmon (2017), also fits alternative models for multivariate Brownian evolution – such as a model in which the rate of evolution for different traits or for different regimes are constrained to be equal, a model in which the evolutionary correlation between traits, r, is constrained to be 0, or alternative matrix decompositions (such as models in which matrix “shape” or “orientation” are permitted to differ between regimes).

An empirical example: South American rock- and non-rock-dwelling lizards

In addition to the centrarchid data, above, we also applied the method to a morphological dataset of South American iguanian lizards (members of the lizard family Tropiduridae sensu lato; Toyama, 2017). For this example, we mapped habitat use of rock-dwelling vs. non-rock-dwelling (Revell et al., 2007; Goodman, Miles & Schwarzkopf, 2008) on a phylogeny of 76 lizard species. Our phylogeny was obtained from Pyron, Burbrink & Wiens (2013), but pruned to contain only the taxa of the present study, and rescaled to have a total length of 1.0. (We rescaled the tree to unit length merely so that our parameter wouldn’t need to be represented using scientific notation. Relative model fits should be completely unaffected by this rescaling.)

To set our regimes, we used a single Maximum Parsimony reconstruction of the discrete trait (rock- vs. non-rock-dwelling) on our phylogeny, in which we fixed all transitions between regimes to be located at the precise midpoint of each edge containing a state change in our reconstruction (Fig. 2). Just as in the centrarchid case, in an empirical study we would probably recommend using multiple reconstructions from a statistical method such as stochastic character mapping (Huelsenbeck, Nielsen & Bollback, 2003), and then averaging the results across stochastic maps (e.g., O’Meara et al., 2006; but see Revell, 2013a and the Discussion for some limitations associated with this common practice). This general workflow is illustrated in the Appendix.

We next fit all eight of the models listed in Table 1 to a dataset consisting of body size and relative dorsoventral body depth from Toyama (2017), both calculated using geometrical definitions for size and shape (Mosimann, 1970; Klingenberg, 2016). Since rock-dwelling has previously been suggested to favor the evolution of dorsoventral flattening (e.g., Revell et al., 2007; Goodman, Miles & Schwarzkopf, 2008), we hypothesized that the evolutionary correlation between body size and depth, while generally positive across this group, could decrease or become negative in rock-dwelling forms due to ecological selection to decouple body depth from size. In Table 2, we show the parameter estimates, model fits, and Akaike weights (see section below) of all eight models from this analysis, sorted by model weight.

Although the weight of evidence is distributed among our top four models in the table, the most notable aspect of all of the best-supported models for these data is that they each allow the evolutionary correlation (r) to differ between the two different mapped regimes on the tree. Models that do not allow the evolutionary correlation to differ by regime (models 1, 2, 2b, and 2c from Table 1) each received less than 10% support.

We found that the evolutionary correlation between body size and size-adjusted body depth was positive in non-rock-dwelling lizards, indicating that larger lizards tended to evolve proportionally greater body depth (Table 2). By contrast, rock-dwelling forms actually exhibited a negative evolutionary correlation between body size and size-adjusted body depth. This is because larger rock-dwelling animals do not tend to evolve proportionally greater body depths. To the contrary, their size-adjusted body depth actually decreases. This is largely consistent with what’s expected given behavioral and biomechanical considerations (Revell et al., 2007; Goodman, Miles & Schwarzkopf, 2008).

A small simulation test of the method

In addition to the empirical applications given above, we tested the method using a small simulation experiment as follows. We first generated forty 100-taxon pure-birth random phylogenetic trees. On each of these trees, we simulated the history of a three-state discrete character. We rejected and repeated any simulation in which any of the three states of the trait was not observed in at least 20 tips. An example simulated tree with evolutionary regimes is given in Fig. 3A.

For all of the forty random trees, we simulated data under each of the eight models of Table 1. To begin each simulation, we first drew values for $\log ({\sigma }_1^2)$ and $\log ({\sigma }_2^2)$ for the two traits from a standard normal distribution (that is to say, ${\sigma }_1^2$ and ${\sigma }_2^2$ were randomly sampled from a log-normal distribution); and we drew a random value or values of the correlation coefficient (r) from a uniform distribution on the interval −1 to 1. Naturally, we sampled different numbers of values for $\log ({\sigma }_1^2)$, $\log ({\sigma }_2^2)$, and r depending on the model that was being used for simulation. For instance, a model with three mapped regimes (e.g., Fig. 3A) and different rates for trait 1, equal rates for trait 2, and different correlations between traits 1 and 2, would involve randomly sampling three values for $\log ({\sigma }_1^2)$, one value for $\log ({\sigma }_2^2)$, and three values for r from their respective distributions. Our simulation procedure does not fix any specific difference in the rates or evolutionary correlations between regimes. Nonetheless, it will on average result in a geometric mean ratio of the highest evolutionary rate over the lowest (for any variable σ² simulation) of around 5.4; and a mean difference between the highest evolutionary correlation and the lowest (for any variable r simulation) of about 1.0.

An example simulated dataset generated using our procedure for different rates (trait 1), and different correlations (model 3b) is shown in Fig. 3B. In this example, we simulated the data using an evolutionary correlation between traits x₁ and x₂ that was positive for regimes 1 and 3, but strongly negative for regime 2 (Fig. 3B).

After completing the numerical simulations, we then proceeded to fit each of the same eight models to each simulated dataset. For each fitted model, we computed AIC and Akaike weights as follows (Akaike, 1974; Burnham & Anderson, 2002).

(4) $$AI{C}_i=2k-2\ln (l_i)w_i={\displaystyle \frac{e^{-\mathrm{\Delta }AI{C}_i/2}}{\mathrm{\Sigma }{e}^{-\mathrm{\Delta }AI{C}_j/2}}}$$

Here, AIC_i is the value of AIC for the ith model; k is the number of parameters in the model; ln(l_i) is the log-likelihood of the ith model; and $\mathrm{\Delta }$AIC_i is the difference in AIC between the ith model and the model with the minimum AIC score in the set. In general, we should prefer the model with the lowest overall value for AIC, and can interpret the Akaike model weights (w), from Eq. (4), as a measure of the strength of evidence supporting each of the models in our set (Akaike, 1974; Wagenmakers & Farrell, 2004).

After fitting all eight models to each of our 40 × 8 = 320 simulated datasets, we next simply calculated the fraction of times in which the generating model was selected as the ‘best’ model (as well as second best, third best, and so on). These results are summarized in Table 3. In general, we found that the generating model tended to be selected as the best or second best model over 86% of the time in simulation, under the simulation conditions described above (Table 3).

In addition, we also calculated the average weight ( $\bar{w}$) of each of the forty datasets for each model. These results are summarized in Fig. 4. This analysis shows that the generating model (in rows) also tended to have the highest average Akaike model weight (in columns; Fig. 4).

For each generating model, the next highest average Akaike model weights tended to be observed in models of similar complexity. For instance, when the generating model was model 4 (no common structure), we found the highest average model weight for model 4 (0.65); and then the next highest average model weights for models 3b (different rates for trait 1, different correlations; 0.14) and 3c (different rates for trait 2, different correlations; 0.13). Likewise, when the generating model was model 1 (common rates, common correlation), we found the highest average model weight for model 1 (0.37), and the next highest average model weight for model 3 (common rates, different correlation; 0.14).

Lastly, having established that the method tends to select the generating model used for simulation, we also proceeded to assess parameter estimation. Here, we computed both the correlation between the generating and estimated values for each model (to measure precision); and the mean difference (to measure accuracy or bias). Correlations for all parameters and all models are given in Table 4. Table 5 shows the mean difference between the generating and estimated parameter values for each model.

Notes on implementation

The model and methods of this study have been implemented for the R statistical computing environment (R Core Team, 2021), and all simulations and analyses for this study were done using R.

The method that we describe in the article is implemented as the function evolvcv.lite of the phytools R package (Revell, 2012). phytools in turn depends on the core phylogenetics R packages ape (Paradis & Schliep, 2019) and phangorn (Schliep, 2011), as well as on a number of other R libraries (Venables & Ripley, 2002; Ligges & Mächler, 2003; Lemon, 2006; Plummer et al., 2006; Chasalow, 2012; Becker et al., 2018; Gilbert & Varadhan, 2019; Azzalini & Genz, 2020; Qiu & Joe, 2020; Warnes, Bolker & Lumley, 2020; Goulet et al., 2021; Pinheiro et al., 2021).

Though the simulations and empirical analyses presented herein use only two or three mapped regimes, the current implementation of this method in phytools, evolvcv.lite, permits an unlimited number of mapped regimes. The reader should keep in mind, however, that the number of parameters to be estimated will rise in direct proportion to the number of regimes.

Limitations

Although we believe that the method of this article extends existing methodology in an interesting way, it comes with its own limitations.

The most evident limitation of the approach described herein and implemented in the phytools evolvcv.lite function is that it is restricted to analyzing no more than two quantitative traits at a time. This practical constraint arises from the fact that we have chosen to parameterize our set of nested models using a common (or different) correlation and common (or different) evolutionary rates paradigm. This paradigm does not extend as easily to more than two traits. The reasons for this are entirely pragmatic. With only two quantitative traits, the likelihood function is defined (and can thus be evaluated) for any values of ${\sigma }_1^2$ and ${\sigma }_2^2>0$ in each regime, and for any −1 < r < 1. For three or more quantitative traits, on the other hand, many combinations of the correlation coefficients r_1,2, r_1,3, and so on, can lead to an evolutionary covariance matrix that is uninvertible (non-positive definite), and thus a likelihood that cannot be computed. This is a solvable problem and one of us [LJR] is presently working to extend evolvcv.lite to an arbitrary number of continuous characters; however, this has required implementing numerical optimization with dynamic bounds and has thus far proved challenging. Ultimately, we hope that this functionality will be added to the phytools R package at some time in the not too distant future.

Likewise, evolvcv.lite only includes multivariate Brownian motion as a model of trait evolution. Expanding beyond Brownian motion, for instance to include the Ornstein-Uhlenbeck (OU) model (e.g., Butler & King, 2004; Beaulieu et al., 2012), is a highly interesting future direction. Even more so than extending to more than two traits, this addition would present a number of additional challenges and is thus beyond the scope of our present work. Not least among these is the substantial increase in model complexity that a multi-regime multivariate OU model would entail. (In particular, each of σ and α, the stochastic and stabilizing selection or ‘rubber band’ parameter in the OU model, respectively, would require a separate m × m matrix for each regime.)

Finally, our approach requires that different regimes on the tree are hypothesized a priori by the investigator. In this way, it follows a paradigm established in Butler & King (2004) and O’Meara et al. (2006), and re-described in a variety of other places since (e.g., Revell, 2008; Beaulieu et al., 2012; Revell et al., 2021; and others). An alternative approach would be to, for example, use reversible-jump Markov Chain Monte Carlo (rjMCMC) to allow our data to inform the location of regime shifts on the tree (e.g., Eastman et al., 2011; Uyeda & Harmon, 2014; also see Revell, 2021 for a different approach using penalized likelihood). In fact, you could go one step further still and sample both the number and position of regime shifts, and model complexity (as given in Table 1) from their joint posterior probability distribution. We feel that this is certainly an intriguing future direction, but out of scope for this relatively modest contribution.

Other considerations

In the current article we have focused on our model for quantitative trait evolution and have left consideration about how regimes are obtained or set on the tree mostly aside. It is entirely appropriate to use this method for circumstances in which there is no ambiguity about the mapping of regimes onto the edges and clades of the tree, such as, for instance, when fitting a model with heterogeneous evolutionary correlations between monophyletic groups or across difference time periods. Often, however, investigators may be inclined to use our approach to study the influence of a mapped discrete character on the evolutionary process for their quantitative traits. In that case, we would recommend considering an approach that allows the explicit incorporation of uncertainty, such as integrating over a set of character histories sampled in proportion to their probability under a model (Huelsenbeck, Nielsen & Bollback, 2003; e.g., Price, Friedman & Wainwright, 2015), as illustrated in the Appendix of this article for our rock- and non-rock-dwelling lizard empirical example. As shown in Revell (2013a), however, this common practice can also result in certain biases, such as causing estimated evolutionary rates to resemble each other more closely than the generating rates (Revell, 2013a; also see May & Moore, 2019). An interesting approach that was suggested by Caetano & Harmon (2019) involves first generating a posterior sample of discrete character histories using stochastic mapping, and then sampling from this set during Bayesian MCMC. This seems worthy of further examination.

In addition, herein we have focused on model selection under simulation conditions in which the true, generating model is invariably included among the set of models we fit to our data. Under these conditions, we tended quite strongly to select the correct model, or if not the correct model, one very similar to it (Table 3; Fig. 4). We do not explicitly consider the inevitable condition in which the true model is not among those of our model set. For other types of comparative methods, model inadequacy can result in inference errors that are quite grave (e.g., Rabosky & Goldberg, 2015; Beaulieu & O’Meara, 2016). Indeed, the general problem of model adequacy is underappreciated in phylogenetic comparative methods (Boettiger, Coop & Ralph, 2012). Some promising approaches to the problem have already been developed (e.g., Pennell et al., 2015; Uyeda, Zenil-Ferguson & Pennell, 2018; Duchene et al., 2018) and we recommend that it be the subject of continued research.

Lastly, we have chosen to concentrate on using maximum likelihood for estimation. In fact, it would be entirely conceivable to take a Bayesian MCMC approach (e.g., following Caetano & Harmon, 2017, 2019), or perhaps even use rjMCMC to integrate over the different models of Table 1 in proportion to their posterior probability, as described in the previous subsection. One advantage of this approach is that it provides a natural framework within which to incorporate uncertainty about the fitted model and its parameters. One disadvantage is that it requires us to decide on appropriate prior distributions for our model parameters, evaluate convergence of the MCMC to the posterior, and summarize the posterior sample (Caetano & Harmon, 2019). It does not obviate any of the other limitations or considerations detailed above.

Relationship to existing methods

The method of this paper is an extension of an existing model described in Revell & Collar (2009). In that article, the authors describe an approach to analyzing heterogeneity in the evolutionary correlation, Brownian rate, and Brownian covariance between different pre-specified branches or clades of a phylogenetic tree. In Revell & Collar (2009) the authors compare two alternative models: one in which evolutionary rates and the evolutionary correlation between characters are free to differ between regimes; and a second in which they are not. Here, we add a set of six additional intermediate models (for two quantitative traits) in which regimes are permitted to share various aspects of their evolutionary process in common (rates for one character or another, evolutionary correlation), while differing in others. This method was also extended to a Bayesian framework by Caetano & Harmon (2017, 2019), and as such the work we present here is also closely related to this research.

Readers of this article who are familiar with phylogenetic comparative methods might observe that it’s also possible to model multivariate trait evolution in which the relationship between traits changes as a function of a discrete factor using a phylogenetic generalized analysis of covariance (Grafen, 1989; Rohlf, 2001; Revell, 2010; Mundry, 2014; Fuentes-G et al., 2016). In this case, we would simply fit a linear model in which a single dependent variable (y) varied as a function of a discrete factor (the tip regime), a continuous variable (x), and their interaction (to permit differences in slope between regimes), while assuming that the residual error in y has a multivariate normal distribution given by the structure of the tree (Rohlf, 2001; Revell, 2010; Fuentes-G et al., 2016). Indeed, this is a valid approach for asking how the relationship between traits changes among lineages of a reconstructed phylogeny. We nonetheless feel that our method adds value for many investigators because it permits an arbitrary (not just tip) mapping of discrete regimes, because it doesn’t require the user to specify dependent and independent variables in the model, because it easily allows us to take into account sampling error in the estimation of species’ means (following Ives, Midford & Garland, 2007), because it’s readily extensible to more than two traits whose correlations might also be expected to change as a function of the mapped regimes, and, finally, because it’s more directly connected to a hypothesized evolutionary process for the traits on our phylogeny (Hohenlohe & Arnold, 2008; Revell & Harmon, 2008).