mvmorph: an r package for fitting multivariate evolutionary models to morphometric data

Introduction

Comparative methods for fitting evolutionary models to continuous data are becoming increasingly popular (e.g. Hansen 1997; Pagel 1999; Butler & King 2004; O'Meara et al. 2006; Thomas, Freckleton & Székely 2006; Hansen, Pienaar & Orzack 2008; Beaulieu et al. 2012; Thomas & Freckleton 2012; Adams 2013; Ingram & Mahler 2013). Over the last years, these methods used to study evolutionary patterns of continuous traits within phylogenies have been widely applied on a trait-by-trait basis for estimating shifts in evolutionary rates associated with environmental global changes, the timing of acquisition of key innovations, or the release of biotic interactions (Collar et al. 2009; Price et al. 2011, 2013; Slater 2013), as well as changes in modes of evolution and adaptation towards different optima (Scales, King & Butler 2009; Collar, Schulte & Losos 2011; Frédérich et al. 2013; Knope & Scales 2013; Mahler et al. 2013). However, phenomic evolution is multivariate by nature, as natural selection operates on functionally related organism's traits (Walsh 2007; Walsh & Blows 2009), and because traits may covary due to pleiotropic effects, genetic linkages or different sources of developmental constraints (Lande & Arnold 1983; Felsenstein 1988; Armbruster & Schwaegerle 1996; Armbruster et al. 2014). Thus, multivariate analyses should logically be favoured over univariate ones (Fig. 1) for describing the complex structure of multiple covarying traits and measurements as routinely studied in morphometrics. Moreover, by considering the evolutionary history of phylogenetically related species, we can further assess how and when phenotypic covariations evolved, either in relation or not with biotic/abiotic events that potentially drive evolution (Fig. 2). Here, we present a new r package, mvmorph, dedicated to the phylogenetic comparative analysis of biological data sets with multiple covarying continuous traits.

Phenotypic multivariate evolution has been thoroughly studied over the last decades, particularly in quantitative genetics in relation to the evolution of genetic covariances (Lande 1979; Lande & Arnold 1983; Cheverud 1984; Felsenstein 1988; Zeng 1988; Endler 1995; Arnold, Pfrender & Jones 2001; Steppan, Phillips & Houle 2002) and has led to particular expectations over macroevolutionary times (e.g. Schluter 1996; Steppan 1997; Arnold, Pfrender & Jones 2001; Baker & Wilkinson 2003; Polly 2004; Hohenlohe & Arnold 2008; Revell & Harmon 2008; Shoval et al. 2012; Goswami et al. 2014). Importantly, modelling multivariate phenotypic evolution through trait variances and pairwise covariances allows describing several aspects (e.g. correlated evolution, correlated selection) that separate univariate analyses cannot provide (Fig. 2). Indeed, univariate analyses are merely limited for describing complex structures and generally show less statistical power than multivariate counterparts (Zheng et al. 2009). Moreover, while complex multivariate data sets are often first reduced through ordination methods such as PCA and a subset of the first principal components is then treated as independent univariate traits, evolutionary model inferences from these transformed spaces can be severely biased (Revell 2009; Monteiro 2013; Polly et al. 2013; Uyeda, Caetano & Pennell 2015). Finally, it is a rather common approach in comparative analysis to regressing traits on each other by assuming, often arbitrarily, that one trait is the dependent variable and the other the independent variable. Nevertheless, such an approach implicitly assumes that one trait is evolutionarily responding to the other (Hansen, Pienaar & Orzack 2008; Revell 2010; Hansen & Bartoszek 2012), a critical assumption of biological causality which is actually relaxed in a multivariate setting focused on mere correlations among response variables.

The multivariate models provided by the mvmorph package will be also of particular interest for evolutionary biologists integrating fossil specimens and taxa in their studies, as the package allows estimating models of changes in evolutionary modes and/or rates at a given point in time (e.g. Slater 2013). Indeed, integration of extinct clades in comparative studies is an important step forward in the study of morphological diversification (Slater, Harmon & Alfaro 2012; Slater 2013): while some studies have pointed towards the need to consider phylogenies in among-group comparisons of morphological rates of evolution (O'Meara et al. 2006; Ricklefs 2006; Adams et al. 2009; Revell & Collar 2009), it appears also necessary to consider underlying processes which may affect the realized variances (Hunt 2012). Actually, some processes such as those leading to shifts from the ancestral state to the optima can be estimated only when including fossil taxa in the analysis (Hansen 1997; Hunt, Bell & Travis 2008). Moreover, although the use of non-ultrametric trees (whether including fossil taxa or fast evolving organisms such as viruses) enhances the accuracy and identifiability of parameter estimates (Slater, Harmon & Alfaro 2012; Ho & Ané 2013), it remains computationally limited to some of the current software implementations (Slater 2014).

Description

Log-Likelihood of Multivariate Models

Most phylogenetic comparative analyses are done by fitting a generalized least square (GLS) model of the form:

(eqn 1)

For a multivariate model with m traits and n species, Y is a vector of length mn and ε is a vector of mn residuals normally distributed in , with V a mn × mn phylogenetic multivariate variance–covariance matrix. If we fit an intercept-only model (i.e. a correlational model where we are estimating the phylogenetic signal in Y), then X is a mn × m design matrix, and β – the (notional) ancestral state expectation – is a vector of parameters to estimate of length m (assuming a unique ancestral state for each trait).

As for the univariate case, the fit of multivariate models of continuous trait evolution is assessed by maximization of the log-likelihood eqn 2 (e.g. Martins & Hansen 1997; Blomberg, Garland & Ives 2003; O'Meara et al. 2006).

(eqn 2)

Evolutionary parameters for the various models of phenotypic evolution [Brownian motion (BM), Ornstein–Uhlenbeck (OU), etc.] enter this equation through the phylogenetic variance–covariance matrix V, while β is given by the GLS solution (eqn 3).

(eqn 3)

Brownian motion

The BM model is the most widespread evolutionary model used in comparative methods (e.g. Felsenstein 1985). Its multivariate extension is straightforward and has been already thoroughly described (Felsenstein 2004; Revell & Harmon 2008; Revell & Collar 2009; Freckleton 2012). The variance–covariance matrix related to a multivariate BM process (eqn 4) is the Kronecker product of the phylogenetic variance–covariance matrix C (for which each element C_ij is the height above the root for the common ancestor of species i and j; see Rohlf 2001; Revell & Collar 2009) with the maximum-likelihood estimate of the evolutionary rate matrix R (Revell & Harmon 2008; Revell & Collar 2009; see also Appendix S1).

(eqn 4)

Ornstein–Uhlenbeck

The OU model is an appealing alternative to the BM model as it describes processes where the traits variance is constrained around one or several optima often referred as ‘selective regime' or ‘adaptive zone' optima at the macroevolutionary scale (Hansen 1997, 2012; Butler & King 2004). In comparative method studies, this process has often been misleadingly interpreted in terms of Lande's (1976) model of population undergoing stabilizing selection in a static landscape (see also discussion in Harmon et al. 2010 and Pennell & Harmon 2013).

The OU model for a multivariate trait matrix Y along a lineage is described by the differential equation:

(eqn 5)

where A and R are square matrices of the strength of selection and the random drift, respectively, β is a vector of mean (optimum) values, and W is the diffusion parameter (Butler & King 2009; Bartoszek et al. 2012; Monteiro 2013). The eigenvalues and eigenvectors of the matrix A describe the path towards the optimum, while the matrix R describes the covariations of the random perturbations. The variance–covariance matrix computed between the m traits at tips i and j can be obtained in two different ways, depending on whether the value at the root is estimated or assumed at the stationary point (see Ho & Ané 2013 for the univariate counterpart). Particularly, when A is symmetric positive definite (i.e. A has only real positive eigenvalues), a stationary solution exists (Gardiner 2004) for which the variance–covariance matrix can be computed as follows:

(eqn 6)

This is the implementation found in the r package ouch (Butler & King 2009). Alternatively, if the covariance matrix is conditioned on the root, a more general form which also allows non-symmetric matrices A may be used (Eq. B.3 in Bartoszek et al. 2012; see also Reitan, Schweder & Henderiks 2012 and Appendix S1 for details):

(eqn 7)

This second approach is particularly useful for fitting multivariate models to non-ultrametric trees (e.g. including fossil species) for which the estimation of the root state is feasible. On ultrametric trees, both approaches should converge to the same estimates if the process is stationary.

As for the univariate case, the expectation of the OU process (eqn 8) can be computed as the weighted average of the optima and the common ancestral state at the root of the tree (Hansen 1997).

(eqn 8)

This expectation is used to estimate the root and the optima by entering the coefficients of their relative contribution in the design matrix X of the GLS equation (eqn 3). The various parameterizations implemented for the root state are described and discussed in the Appendix S1.

Early Burst

The mvmorph package allows computing a joint estimate for a decrease or increase in evolutionary rates in a multivariate data set, following Blomberg, Garland & Ives (2003) exponentially accelerating or decelerating (ACDC) model or Harmon et al.'s (2010) Early Burst (EB) model. These models have been used for studying adaptive radiations where evolutionary theory predicts a burst of evolutionary rates and an increase in morphological disparity, followed by a decrease in rates of morphological evolution as ecological niches are filled (Cooper & Purvis 2010; Harmon et al. 2010; Slater et al. 2010; Ingram, Harmon & Shurin 2012; Slater 2013). Here, the phylogenetic matrix C is transformed by the exponential transformation (g parameter) of the ACDC model (Blomberg, Garland & Ives 2003), which is jointly estimated during likelihood maximization of the evolutionary rate matrix (eqn 9).

(eqn 9)

Shifts models

Slater (2013) recently introduced a heterogeneous model allowing for the shift from one evolutionary process to another in the same tree. Dealing with heterogeneous models was also previously proposed by using tree-stretching approaches (O'Meara 2012; but see Slater 2014). Such approaches are particularly appealing when mixing fossil and extant data to assess changes in evolutionary modes in addition to temporal changes in evolutionary tempo (Hunt, Bell & Travis 2008; Slater & Harmon 2013; see Appendix S1). A combination of the various models implemented in mvmorph can be fitted to two successive time slices or mapped era using the mvSHIFT function (Table 1). For instance, the variance–covariance for the ‘ecological release' model proposed by Slater (2013), combining an OU process then followed by a BM process, is computed as follows:

(eqn 10)

where Cbm_ij represents the history between the taxa i and j shared during the time period where the BM process took place; Cou_ij represents the shared history between both taxa accrued during the OU process; and Cou_jj and Cou_ii represent the evolutionary histories (variances) accrued in each lineage during the OU process.

Constraining the parameter spaces

Estimating parameters of multivariate models is a complex task because matrices of parameters must always fulfil specific conditions, for example positive definiteness for variance–covariance matrices. Generally, one of the efficient ways to solve this issue is to use an unconstrained optimization of the likelihood with a parameterization that enforces the matrices to a particular structure (Pinheiro & Bates 1996). Various parameterizations are used in mvmorph (see Appendix S1) for ensuring positive definiteness of rate matrices (Pinheiro & Bates 1996), conditioning matrices on their eigenvalues (e.g. Sy, Taylor & Cumberland 1997; Jaffrézic, Thompson & Pletcher 2004; Bartoszek et al. 2012), or testing significance of interactions between traits by constraining the parameter space and taking advantage of joint likelihood optimization (e.g. Revell & Collar 2009; Adams 2013).

Alternatives

Other r implementations of multivariate comparative approaches already exist that calculate parts of the same or alternative models to those provided in mvmorph. Multivariate extensions of Pagel (1999) branch-length transformations are proposed in motmot (Thomas & Freckleton 2012; no more on the cran repository but accessible from the archive); BM with multiple evolutionary rate matrices can be fitted with evol.vcv as well as a phylogenetic canonical correlation analysis with phyl.cca in phytools (Revell 2012), and geomorph allows computing BM rates for high-dimensional data sets (Adams & Otarola-Castillo 2013). Multivariate OU is available from the ouch package (Butler & King 2009) as well as mvslouch that allows fitting multivariate OU models with a randomly evolving predictor variable (Bartoszek et al. 2012).

Drawbacks

Worked examples

To illustrate how mvmorph works, we describe and comment two simulated examples in the Appendix S2 (also included as vignettes in the mvmorph package). In the first example, we show how multivariate models can out-compete separate univariate analyses for a phylogeny involving two distinct selective regimes corresponding to two separate species groups. In the second worked example, we illustrate how to assess significances of differences in evolutionary covariations by testing a hierarchy of nested models. In addition, the help pages embedded within the mvmorph package also provide several simulated examples illustrating how the various functions work.

Conclusion

Elaboration of the mvmorph package was initially motivated by the need for a unified r package for fitting multivariate models of phenotypic evolution incorporating measurement errors and using flexible ways to map evolutionary scenarios for morphometric data. Obviously, such objectives can be extended to any biological data set made of covarying continuous traits, as well as univariate traits. mvmorph also implements models that were not previously designed for the multivariate case, as well as new ones. Last, but not least, the fast likelihood implementations provided by the mvmorph package may be particularly appealing when dealing with large data sets and assessing the statistical power of comparative methods through modelling approaches (Boettiger, Coop & Ralph 2012). The mvmorph package is in continuous development and can be followed and contributed from the gitHub website (https://github.com/JClavel/mvMORPH). mvmorph (current version: 1.0.5) is an open-source r package available for download on the cran repository (http://cran.r-project.org/web/packages/mvMORPH/).