Quantifying and Comparing Phylogenetic Evolutionary Rates for Shape and Other High-Dimensional Phenotypic Data

.—Many questions in evolutionary biology require the quantiﬁcation and comparison of rates of phenotypic evolution. Recently, phylogenetic comparative methods have been developed for comparing evolutionary rates on a phylogeny for single, univariate traits ( (cid:2) 2 ), and evolutionary rate matrices ( R ) for sets of traits treated simultaneously. However, high-dimensional traits like shape remain under-examined with this framework, because methods suited for such data have not been fully developed. In this article, I describe a method to quantify phylogenetic evolutionary rates for high-dimensional multivariate data ( (cid:2) 2mult ), found from the equivalency between statistical methods based on covariance matrices and those based on distance matrices (R-mode and Q-mode methods). I then use simulations to evaluate the statisticalperformanceofhypothesis-testingproceduresthatcompare (cid:2) 2mult fortwoormoregroupsofspeciesonaphylogeny. Under both isotropic and non-isotropic conditions, and for differing numbers of trait dimensions, the proposed method displays appropriate Type I error and high statistical power for detecting known differences in (cid:2) 2mult among groups. In contrast, the Type I error rate of likelihood tests based on the evolutionary rate matrix ( R ) increases as the number of trait dimensions ( p ) increases, and becomes unacceptably large when only a few trait dimensions are considered. Further, likelihood tests based on R cannot be computed when the number of trait dimensions equals or exceeds the number of taxa in the phylogeny (i.e., when p ≥ N ). These results demonstrate that tests based on (cid:2) 2mult provide a useful means of comparing evolutionary rates for high-dimensional data that are otherwise not analytically accessible to methods based on the evolutionary rate matrix. This advance thus expands the phylogenetic comparative toolkit for high-dimensional phenotypictraitslikeshape.Finally,Iillustratetheutilityofthenewapproachbyevaluatingratesofheadshapeevolutioninalineageof Plethodon salamanders.[Evolutionaryrates;geometricmorphometrics;macroevolution;morphologicalevolution; phylogenetic comparative method.].

What is the tempo of evolutionary change, and how rapidly do phenotypic traits diversify? For decades evolutionary biologists have examined the tempo of evolution as one means to understand the accumulation of phenotypic diversity and the processes responsible for evolutionary diversification (Simpson 1944;Gingerich 1993;Foote 1997;Sidlauskas 2008;Gingerich 2009). Numerous hypotheses have been proposed to explain the observation that lineages differ in their phenotypic evolutionary rates (e.g., Sepkoski 1978;Blomberg et al. 2003;Butler and King 2004;Collar et al. 2009) and predict how such differences may result in distinct macroevolutionary trends across groups of organisms (e.g., Gingerich 1993;Harmon et al. 2003). Indeed, the exploration of evolutionary rates has been a perennial component of macroevolution since the field's inception.
In recent years a renaissance has occurred in the study of evolutionary rates which has coincided with rapid analytical developments in phylogenetic comparative biology. These analytical approaches estimate the net rate of phenotype change over time along a phylogeny ( 2 ), based on some underlying model that characterizes the evolutionary process (described below). With these methods, rates of phenotypic evolution can be quantified and used to identify differences in evolutionary rates among clades or shifts in evolutionary rates within clades (O'Meara et al. 2006; Thomas et al. 2006;Revell 2008;Revell and Harmon 2008;Eastman et al. 2011;Beaulieu et al. 2012;Revell 2012). Methods for comparing evolutionary rates between traits on a phylogeny have also been developed (e.g., Adams 2013). Together, these approaches are part of the ongoing mathematical unification of evolutionary rates, models of phenotypic evolution, and evolutionary hypothesis testing. This union represents the current forefront of phylogenetic comparative biology.
Most empirical studies examine rates of evolution for single, univariate traits such as body size (e.g., Thomas et al. 2009;Harmon et al. 2010), or some composite measure derived from multiple traits, such as principal component scores (e.g., Collar et al. 2009;Mahler et al. 2010;Betancur-R et al. 2012). However, organismal phenotypes can also be characterized multivariately, either by a set of traits treated simultaneously or by complex, multi-dimensional traits such as shape. In these instances, evolutionary change corresponds to a shift in the position of a species in a multivariate trait space whose axes correspond to trait dimensions. Figure 1 shows the changes over time for a bivariate trait evolving along a phylogeny (examples along ancestor-descendant sequences are found in : Polly 2004;Bookstein 2013). For multivariate traits, evolutionary rates have been approximated using Mahalanobis distances between pairs of extant taxa (Arnegard et al. 2010;Carlson et al. 2011), or by describing the net change per unit time along ancestor-descendent lineages (Polly 2008; see also Gingerich 2009). However, these approaches do not provide a single estimate of the evolutionary rate for an entire clade of species in a phylogenetic context, as is found in the univariate case ( 2 ). Alternatively, one may characterize the tempo of evolution for multivariate data using the evolutionary rate matrix (R), which is the algebraic extension of 2 (Revell and Harmon 2008). However, while this approach extends the concept of phylogenetic evolutionary rates to a multivariate context, several previously unappreciated challenges exist when implementing rate matrix methods with highdimensional data. For instance, when the number of trait dimensions equals or exceeds the number of taxa in the phylogeny (p ≥ N), likelihood tests based on R cannot be accomplished, because the matrix computations used to obtain the likelihood are singular (specifically, (R⊗C) −1 is singular). Thus, while evolutionary biologists increasingly characterize organismal phenotypes using highly multivariate data, likelihood methods cannot be used to evaluate rates of evolution for these data sets, as the likelihood of both the null model and its alternative cannot be determined. Second, as I demonstrate below, Type I error rates of likelihood tests based on R increase as the number of trait dimensions (p) increases, and become unacceptably large when only a few trait dimensions are considered simultaneously. This severely limits the utility of rate matrix methods when evaluating evolutionary hypotheses with multivariate data, as it is not known whether a significant finding represents actual differences in evolutionary rate matrices or a Type I error. Clearly, for high-dimensional phenotypic data like shape, an alternative procedure is required.
In this article, I propose a method for estimating the evolutionary rate of change of high-dimensional phenotypic data in a phylogenetic context. The approach is developed for phenotypic data derived under evolutionary models commonly used in phylogenetic comparative biology (described below), and is found from the equivalency between statistical methods based on covariance matrices and those based on distance matrices. Consistent with this equivalency, the distance-based approach provides estimates of evolutionary rates that are numerically identical to those obtained using standard covariance-based implementations when implemented on univariate data. Further, when used on multivariate data, the approach is not constrained by trait dimensionality as is the case with rate matrix methods. I provide a hypothesis-testing procedure for comparing multivariate evolutionary rates on a phylogeny, and through simulation show that this procedure has acceptable Type I error and appropriate statistical power for detecting differences in rates among groups, under both isotropic and non-isotropic conditions. Finally, a biological example is presented which demonstrates the utility of the approach. Computer code written in R for implementing the procedures is also provided.

MODELING PHENOTYPIC EVOLUTION ALONG A PHYLOGENY
Statistically evaluating phenotypic trends on a phylogeny first requires that one specify the underlying model that describes the evolutionary process that generated the data. For continuous data, the most commonly utilized model is Brownian motion (BM; Edwards and Cavalli-Sforza 1964;Felsenstein 1973;Felsenstein 1981;Felsenstein 1985;Felsenstein 2004). Here, phenotypic changes are assumed to be independent from time step to time step, with a mean displacement of zero and a variance ( 2 ) proportional to time (Felsenstein 1973;Felsenstein 2004). Thus, for a set of taxa related by a phylogeny, the expected net evolutionary change in a phenotypic trait from the ancestral state is zero, and the expected phenotypic variance among taxa increases linearly with time ( 2 t). Because 2 describes the rate at which phenotypic variation accumulates per unit time (Felsenstein 1985), it is commonly referred to as the evolutionary rate parameter (Martins 1994;Garland and Ives 2000;O'Meara 2012). When multivariate data are considered, each trait dimension has an expected net change of zero, but 2 may differ for each trait dimension (Felsenstein 1988). Additionally, changes at each time step may be correlated across trait dimensions. Thus, for multivariate data, the BM process is described by a covariance matrix (R), whose diagonal elements represent the evolutionary rate for each trait dimension ( 2 ), and whose off-diagonal elements express the covariation in changes among dimensions (Felsenstein 1988). Because R is the multivariate generalization of 2 , it is commonly referred to as the evolutionary rate matrix (Revell and Harmon 2008;Revell and Collar 2009 Given an evolutionary model, it is possible to estimate the net rate of phenotypic evolution over time along a phylogeny for a set of related taxa. For BM models, several implementations have been proposed. First, one can utilize phylogenetic generalized least squares to estimate 2 (PGLS: Grafen 1989;Martins and Hansen 1997). Here, the least-squares estimate is equivalent to the maximum-likelihood estimate (Felsenstein 1973;Garland and Ives 2000;O'Meara et al. 2006), and is found as where Y is an N ×1 vector of phenotypic values for the N species, E(Y) is an N ×1 vector containing the phylogenetic mean at the root of the phylogenŷ a = (1 t C −1 1) t (1 t C −1 Y), and C −1 is the inverse of the N ×N phylogenetic covariance matrix. Time is included in this calculation via the matrix C, whose diagonal elements contain the phylogenetic distance from the root of the tree to the tips, and the off-diagonals contain the phylogenetic distances from the root of the tree to the most recent common ancestor for each pair of species (Martins and Hansen 1997;Garland and Ives 2000;Rohlf 2001). Thus, one can view 2 as a phylogenetically standardized variance. Note that the unbiased estimate of 2 may be obtained by dividing by (N −1) rather than N (Garland and Ives 2000;O'Meara et al. 2006). For a set of traits treated simultaneously, equation (1) results in an evolutionary rate matrix (R), as described above. A second approach for estimating evolutionary rates is based on independent contrasts (Felsenstein 1985, Garland 1992Garland and Ives 2000;Ackerly 2009). Here, 2 is obtained from the cross product of the vector containing the phylogenetically independent contrasts among taxa on a phylogeny, divided by N (Blomberg et al. 2003). This approach yields numerically identical estimates to those found using equation (1) above (Garland and Ives 2000).
Evolutionary rates can also be obtained via phylogenetic transformation (Blomberg et al. 2003). Here, the phenotypic data are first transformed by the phylogeny as where Y and E(Y) are N ×1 vectors as described above, and D is obtained from an eigen-decomposition of phylogenetic covariance matrix C (DCD = I : Garland and Ives 2000; Blomberg et al. 2003). Equation (2) rotates the phenotypic data to the principal eigenvectors of the phylogenetic covariance matrix, and describes them by a set of scores projected on those axes. From the transformed data, 2 is then estimated as Notice that equations (2) and (3) differ slightly from what was originally presented in Blomberg et al. (2003).
The reason is that the phylogenetic mean (â) must be used in the transformation step of this procedure (equation (2)), not in the variance-estimation step as was originally presented (Ives T., personal communication). Importantly, for the same phenotypic data, phylogeny, and evolutionary model (BM), identical estimates of 2 are obtained using PGLS, phylogenetically independent contrasts, and phylogenetic transformation (see, e.g., Garland and Ives 2000).

AN EVOLUTIONARY RATE FOR HIGH-DIMENSIONAL TRAITS
Evolutionary rates quantify the accumulation of variance in a trait over time while accounting for the phylogenetic relationships among species. As such, the variance-based equations described above represent logical procedures for estimating parameters like 2 . Statistically, methods that summarize data matrices by their variances and covariances are referred to as R-mode techniques (sensu Legendre and Legendre 1998). For any centered data matrix, Y, the covariance matrix may be found from the inner product of (Y t Y) divided by N-1 (Krzanowksi 1993;Rencher 1995). However, data matrices can also be summarized by the matrix of pairwise distances among objects (Qmode techniques), which may be derived from the outer product of Y: (YY t : Krzanowksi 1993, Rencher 1995. Importantly, for Euclidean data there exists a relationship between statistical approaches based on covariance matrices and those based on distance matrices, such that empirical results from Q-mode and R-mode techniques may be numerically identical. One demonstration of this property is found with multivariate ordination techniques. Here it can be shown that a principal coordinates ordination based on the centered distance matrix between objects is identical to that obtained from a principal components analysis based on the covariance matrix of the same Euclidean data (Gower 1966; see also Krzanowksi 1993;Legendre and Legendre 1998). In addition, it has also been shown that sums of squares obtained for factors in linear models are identical when obtained using standard covariance-based methods (ANOVA and regression), or when derived from distance-based procedures (for a mathematical derivation of this property and examples, see Anderson 2001;McArdle and Anderson 2001). Thus, when implemented properly, statistical summaries based on covariances and those based on distances can yield identical results for Euclidean data sets.
In light of these observations, I propose a Q-mode approach to estimate multivariate evolutionary rates ( 2 mult ) based on distances rather than covariances. The method is suitable for high-dimensional multivariate data sets whose evolution is best considered as occurring in a multivariate trait space (Fig. 1). The method assumes that variation among species accumulates over time following a BM model of trait evolution as described above (Felsenstein 1973 (Butler et al. 2000;Blomberg et al. 2003).
To obtain this estimate, the phenotypic data are first transformed by the phylogeny (sensu Garland and Ives 2000;Blomberg et al. 2003): In this case, Y is an N ×p matrix of phenotypic trait values (N species by p dimensions), and E(Y) is an N ×p matrix of multivariate phylogenetic means. Next, the Euclidean distance between each phylogenetically transformed species (the rows of U Y ) and the origin of the data space is calculated, and these are concatenated into an N ×1 vector containing the Euclidean distances from each species mean to the origin (PD U,0 ). Finally, an estimate of the evolutionary rate is obtained from the sum of squared distances between the phylogenetically transformed data and the origin, PD U,0 . This is found by re-expressing equation (3) as This value is then divided by the number of trait dimensions (p) so that 2 mult is expressed as the average rate across trait dimensions.
For univariate data, it can be shown that estimates of evolutionary rates obtained using this Q-mode procedure are numerically identical to those obtained from the variance-based methods typically used. A demonstration of this property is found in Appendix 1. Thus, for univariate data, the distance-covariance equivalency has been preserved, since estimates of 2 and 2 mult are identical. When used on multivariate data, the Q-mode formulation yields a single evolutionary rate that describes the net overall tempo of phenotypic evolution per unit time in the multi-dimensional trait space. This is conceptually similar to methods that use Mahalanobis distances (Polly 2008;Arnegard et al. 2010;Carlson et al. 2011), but provides a single estimate of the evolutionary rate of change for the entire clade of species along their phylogeny, as in 2 .
Further, the Q-mode approach alleviates the analytical challenges found when methods based on the evolutionary rate matrix are implemented on high-dimensional data. Specifically, when the number of trait dimensions equals or exceeds the number of taxa in the phylogeny (p ≥ N), likelihood calculations based on R will be singular. Indeed, such situations are likely to be common when examining phenotypic traits like shape, as the number of dimensions required to adequately represent the phenotypic trait can be quite large, and can easily exceed the number of species in a lineage (e.g., McPeek et al. 2008;Klingenberg and Gidaszewski 2010). In contrast, 2 mult and hypothesis tests based on it (below) can always be performed because the number of taxa, and thus the number of distances to be calculated, remains the same regardless of trait dimensionality. Therefore, the Q-mode approach proposed here provides a means of testing evolutionary hypotheses of rate shifts for high-dimensional phenotypic data that are not possible to evaluate if methods based on the evolutionary rate matrix are utilized.

HYPOTHESIS-TESTING PROCEDURES
Having developed a means of estimating multivariate evolutionary rates, it is now of interest to determine whether evolutionary rates differ for two or more groups of taxa. Addressing this question requires two things; a test statistic for comparing evolutionary rates and a procedure for statistically evaluating it. In phylogenetic comparative biology, a number of statistical procedures have been implemented for comparing evolutionary rates and rate matrices. Frequently, likelihood methods are used, where the fit of the data to the phylogeny under a model with a single evolutionary rate is compared with the fit of a model containing multiple evolutionary rates (e.g., O'Meara et al. 2006;Thomas et al. 2006;Revell 2008). Bayesian MCMC methods (Eastman et al. 2011;Revell 2012), and modified t-tests of phylogenetically independent contrasts (Garland 1992;McPeek 1995) have also been proposed. More broadly, evolutionary hypotheses can be evaluated using a null distribution of possible outcomes generated by phylogenetic permutation (Blomberg et al. 2003;Klingenberg and Gidaszewski 2010) or phylogenetic simulation (Garland et al. 1993;Boettiger et al. 2012).
Here I use phylogenetic simulation to test hypotheses that compare multivariate evolutionary rates among groups.
To determine if two groups of species on a phylogeny display distinct evolutionary rates, the following procedure is used. First, the phenotypic data for all species are transformed by the phylogeny using equation (4), and the Euclidean distances of each species to the origin are obtained (PD U,0 ). Next, multivariate evolutionary rates ( 2 mult ) are estimated for each group of taxa with equation (5), using the subset of PD U,0 corresponding to each group divided by the number of taxa in that group. The rates are then rank-ordered, and the ratio between them is obtained ( 2 mult.A / 2 mult.B ) such that: 2 mult.A > 2 mult.B . This ratio represents the relative difference in evolutionary rates between groups and serves as the test statistic for hypothesis testing. The observed ratio of evolutionary rates is then statistically compared with a distribution of possible ratios obtained under the null hypothesis of no rate difference between groups (i.e., under the hypothesis that a single evolutionary rate exists for the entire phylogeny). To obtain this distribution, the multivariate evolutionary rate for the entire phylogeny ( 2 mult.Tot ) is found for the observed data. Next, phenotypic data sets are simulated along 170 SYSTEMATIC BIOLOGY VOL. 63 the phylogeny, where 2 mult.Tot is used as the input rate for stochastic simulations. For each simulated data set, 2 mult is estimated for both groups, and the ratio between them is then obtained. The proportion of simulated ratios greater than the observed is the probability that the two evolutionary rates differ statistically from one another (for implementation details, see Appendix 2 and the Supplemental Material available on Dryad (http://dx.doi.org/10.5061/dryad.41hc4) upon publication). If more than two groups of taxa are to be compared, the method can be generalized in the following manner. First, 2 mult is obtained for all groups of taxa, and the N(N −1)/2 ratios of evolutionary rates are obtained for each pair of groups. Next, the average ratio is obtained and treated as the observed test statistic for the data set. Phylogenetic simulations proceed exactly as described above to generate a distribution of possible average evolutionary rate ratios. The proportion of simulated ratios greater than the observed is the probability that the evolutionary rate of at least one group differs from the rest. Thus, for three or more groups, a statistically significant result is interpreted in a manner equivalent to that found in ANOVA. Computer code in R for implementing these procedures is found in Appendix 2.

STATISTICAL PERFORMANCE
To evaluate the statistical performance of the hypothesis-testing procedure proposed above, I executed a series of computer simulations. The multivariate data used in these simulations were generated using two different patterns of error covariance: isotropic error and non-isotropic error. Initial simulations were conducted on a balanced phylogeny containing 32 species, divided into two monophyletic groups (subclades) of 16 species each. For each simulation, the number of trait dimensions was first selected (P = 2, 3, 4, 5, 7, 10). Next, input error covariance matrices of dimension p×p were constructed for each group, which were used to generate the phenotypic data. For simulations assuming isotropic error, the evolutionary rate for the first group was set to 2 1 = 1.0 for each trait dimension. For the second, group, the evolutionary rate for all trait dimensions was set to a fixed proportional difference relative to that in the first group. These varied across simulations such that the initial rate difference was known between groups ( 2 2 = 1.0, 1.5, 2.0, 3.0, 4.0). For simulations assuming non-isotropic error, the evolutionary rate for the first group was drawn from a normal distribution ( = 1.0; = 0.1) for each trait dimension. The evolutionary rates for all traits for the second group were then obtained by multiplying these values by a constant (k = 1.0, 1.5, 2.0, 3.0, 4.0) to obtain a known initial rate difference between groups.
From these initial covariance matrices, 1000 phenotypic data sets were obtained for each monophyletic group by evolving multi-dimensional traits along the phylogeny following a BM model of evolution. For each data set, evolutionary rates ( 2 mult ) were then estimated for each phylogenetic group, and compared with one another using the hypothesistesting procedure described above. The proportion of significant results (out of 1000) was then treated as an estimate of the Type I error (when 2 2 = 2 1 ) or statistical power (when 2 2 > 2 1 ) of the test. Finally, for simulations evaluating Type I error ( 2 2 = 2 1 ), rate comparisons were also performed using likelihood methods based on the evolutionary rate matrix (R) to evaluate the effect of trait dimensionality on the performance of this test.
Simulations were also performed across a wider set of conditions to evaluate the robustness of the method proposed here. These simulations evaluated the effect of the number of taxa in the phylogeny (N = 16, 32, 64, 128), as well as the effect of randomly generated phylogenies on statistical performance (N = 16, 32, 64, 128). In addition, simulations were also performed on random phylogenies containing three groups of taxa, to evaluate the performance of the method for comparing more than two groups. Additional implementation details and results from all simulations are found in the Supplementary Material, available at http://dx.doi.org/10.5061/dryad.41hc4.
Results: For all simulations, hypothesis tests based on 2 mult displayed appropriate Type I error rates near the nominal value of = 0.05 when initial error covariance matrices for the two groups were identical ( 2 2 = 2 1 ). This pattern remained consistent across the range of trait dimensionality examined in this study (Fig. 2). In contrast, tests based on the evolutionary rate matrix (R) had unacceptably high Type I error rates under these FIGURE 2.
Simulation results evaluating the Type I error of hypothesis-testing procedures that compare evolutionary rates on a phylogeny. Data were simulated under isotropic conditions on a balanced phylogeny containing 32 species. Evolutionary rates between two groups within the phylogeny were compared. Type I error rates for methods based on the evolutionary rate matrix R and 2 mult are shown relative to increasing levels of trait dimensionality. Simulation results evaluating the statistical power of hypothesis-testing procedures that compare evolutionary rates on a phylogeny. Data were simulated on a balanced phylogeny containing 32 taxa, and evolutionary rates between two groups within the phylogeny were compared. a) Statistical power curves for methods based on 2 mult for data generated under isotropic conditions. Relative rate differences in initial rate matrices between groups are shown on the abscissa, and power is on the ordinate. Curves for increasing levels of trait dimensionality are shown. b) Statistical power curves for data generated under non-isotropic conditions. All other variables as described in panel a).
conditions, and Type I error of these tests increased as a function of trait dimensionality (Fig. 2). With the latter approach, Type I error rates were over 8% for tests on bivariate traits, nearly 14% when p = 4, and exceeded 50% when p = 10. Prior individual-based simulations evaluating this approach yielded Type I error rates significantly higher than = 0.05 (Revell and Harmon 2008), though those simulations did not investigate data sets as highly dimensional (p = 10) as the ones examined here. It should also be recalled that when the number of trait dimensions equals or exceeds the number of taxa in the phylogeny (p ≥ N), likelihood tests based on R cannot be executed because the matrix computations for estimating the likelihood are singular. Together, these findings reveal that methods based on the evolutionary rate matrix may not be the optimal approach to utilize for detecting shifts in evolutionary rates in highly dimensional multivariate data sets.
When 2 2 > 2 1 , the statistical power of tests based on 2 mult increased rapidly as the relative difference in evolutionary rates between groups increased. This was the case when data were simulated under both isotropic (Fig. 3a) and non-isotropic conditions (Fig. 3b). Statistical power also increased with increasing trait dimensionality, attaining very high power (>0.8) when p = 5, even for relatively small relative rate differences ( 2 2 / 2 1 = 2.0). Similar results were obtained on phylogenies with different numbers of taxa, on random phylogenies, and for the comparison of evolutionary rates for three groups of species (Supplementary Material, available at http://dx.doi.org/10.5061/dryad.41hc4.). Power curves were not generated for likelihood tests based on R because power estimates are not interpretable for methods with high Type I error rates. Overall these simulations reveal that hypothesis tests based on 2 mult provide a powerful means of detecting differences in evolutionary rates for high-dimensional phenotypic traits under both isotropic and non-isotropic conditions, and that 2 mult represents an important analytic advance over existing approaches when used on multivariate data sets.

A BIOLOGICAL EXAMPLE
To illustrate the method described above, I provide a biological example comparing evolutionary rates for a high-dimensional trait (head shape) in a lineage of Plethodon salamanders. Salamanders of the genus Plethodon are long-lived terrestrial species found in North American forests (Highton 1995). Considerable ecological and behavioral work has demonstrated that interspecific competition plays an important role in dictating species ranges and community composition at local scales (e.g., Jaeger 1970;Jaeger 1971;Hairston 1980;Anthony et al. 1997), and affects community structure at regional scales ). Further, competitive interactions often result in evolutionary shifts in morphology within and between species, particularly in head shape (e.g., Adams et al. 2007;Arif et al. 2007;Adams 2010;Deitloff et al. 2013). Head shape also displays a strong genetic component (Adams 2011), enabling selection to generate heritable microevolutionary changes. Within Plethodon, the P. cinereus subclade comprises several small-bodied species, distributed throughout eastern North America (Fig. 4a). Several of these species are threatened or endangered, and have restricted distributions that limit them to high-elevation habitats (e.g., P. hubrichti, P. nettingi, P. shenandoah, P. sherando : Highton 1999;Highton 2004). Work in other taxa has shown that species with restricted distributions can display elevated 172 SYSTEMATIC BIOLOGY VOL. 63 FIGURE 4. a) Geographic distributions of species in the P. cinereus subclade (data from Global Amphibian Assessment). b) Positions of 11 anatomical landmarks used to quantify head shape in Plethodon salamanders (image from Adams et al. 2007). c) Fossil-calibrated molecular phylogeny displaying the estimated phylogenetic relationships among species in the P. cinereus subclade (from Wiens et al. 2006). d) Plot of phylomorphospace viewed as the first two principal component axes of tangent space. rates of morphological evolution (e.g., Millien 2006;Harmon et al. 2008). Thus it is of interest to determine whether these species display elevated rates of morphological evolution as compared with other species in the lineage.
To test this prediction I quantified head shape from 478 adult specimens, representing all species in the P. cinereus subclade (data from previously published studies). Head shape was quantified using geometric morphometric methods (Bookstein 1991;Adams 2013). First, 11 two-dimensional landmark coordinates were digitized from images of the left-lateral side of each head (Fig. 4b), and variation in the position of the jaw relative to the skull was standardized mathematically by rotating the jaw so that the articulation angle was invariant among specimens. A Generalized Procrustes analysis was then performed to align the specimens to a common coordinate system (Rohlf and Slice 1990), and to obtain a set of Procrustes tangent coordinates.
These coordinates describe the shape of each specimen, which is represented as a point in a high-dimensional trait space. For each species, the mean head shape was obtained, and the rate of head shape evolution was then evaluated using a time-calibrated molecular phylogeny for Plethodon, which contained 9 of the 10 members of this subclade (Wiens et al. 2006: Fig. 4c). Using this phylogeny, I obtained two evolutionary rates: one for the geographically restricted species (not including P. sherando which was not represented in the phylogeny), and a second for the remaining species in the clade. The ratio of the two evolutionary rates was then obtained as describe above, and evaluated using phylogenetic simulation with 9999 iterations. In addition, because several alternative phylogenetic hypotheses exist for the P. cinereus subclade (Highton 1999;Sites et al. 2004), evolutionary rates for each group of taxa were estimated and evaluated using these alternative phylogenetic hypotheses for comparison. Finally, phenotypic evolution was visualized in phylomorphospace (Sidlauskas 2008), where the extant taxa and the phylogeny were projected into the morphological trait space, and visualized along the first two axes of this space using principal components analysis. All analyses were performed in R 3.0.1 (R Core Development Team 2013) using routines in the library geomorph (Adams and Otárola-Castillo 2012;Adams and Otárola-Castillo 2013), and new routines written by the author (Appendix 2).
Results: Phylogenetic simulation revealed that the net evolutionary rates over time for the two groups differed significantly ( 2 mult.A / 2 mult.B = 1.8381; p sim = 0.0052). Specifically, the rate of head shape evolution was nearly two times faster in the geographically restricted species than in the remaining species in the clade ( 2 restricted = 7.268×10 −5 ; 2 rest = 3.954×10 −5 ). Concordant patterns were obtained using alternative phylogenetic hypotheses for the group (results not shown). When viewed in phylomorphospace, patterns of head shape appeared to evolve in a diversifying fashion, where the morphology of extant taxa emanated from a central morphological point as exemplified by the hypothesized ancestors (Fig. 4d). Further, there was little evidence of repeated or convergent evolution in head shape, as phylogenetic branches did not cross one another, and more distantly related taxa were not found in close proximity when viewed in the multivariate trait space. Biologically, these observations suggest that the geographically restricted species have elevated rates of morphological evolution as compared with their more wide-ranging congeners. This finding is consistent with the hypothesis that morphological evolution may be more accelerated in populations and species with reduced ranges (see, e.g., Millien 2006;Harmon et al. 2008). Finally, it should be noted that for this example the number of species (9) was less than the number of trait dimensions (18). Thus, this interesting biological observation could only be obtained by using the Q-mode rate method ( 2 mult ) described here.

DISCUSSION
A common feature in macroevolution is that rates of phenotypic evolution differ among clades and among traits. In this article, I described a Q-mode approach for estimating phylogenetic evolutionary rates for multivariate data ( 2 mult ) using distances, and proposed hypothesis tests for comparing two or more multivariate evolutionary rates on a phylogeny. The approach complements existing phylogenetic comparative methods by allowing high-dimensional data sets like shape to be evaluated in a phylogenetic framework. Using simulations, I demonstrated that the method has appropriate Type I error rates and high statistical power for detecting known differences in 2 mult among groups, for data simulated under both isotropic and non-isotropic conditions. Further, I showed that likelihood methods based on the evolutionary rate matrix display unacceptably high Type I error rates for these same data conditions, and in some circumstances cannot be analytically completed (when p ≥ N). These results have several important implications for the study of evolutionary rates and how they should be evaluated. First, the evolutionary rate parameter developed here greatly expands the scope of phylogenetic comparative biology by allowing a much wider class of phenotypic attributes to be investigated. Because current methods were largely developed in a univariate context, nearly all studies of phylogenetic evolutionary rates for clades have been restricted to the analysis of single univariate traits (e.g., Collar et al. 2009;Thomas et al. 2009;Mahler et al. 2010;Price et al. 2010;Valenzuela and Adams 2011), or in a few cases, sets of traits (Revell and Collar 2009;Kozak and Wiens 2010;Rabosky and Adams 2012; see also Martin and Wainwright 2011). However, many evolutionarily important phenotypic attributes are also highly dimensional; yet evolutionary rates for these traits cannot be reliably compared with standard approaches, because hypothesis tests based on the evolutionary rate matrix have unacceptably high Type I error rates or cannot be algebraically computed. In contrast, hypothesis tests based on 2 mult are not adversely affected by trait dimensionality, and have appropriate power to detect differences in evolutionary rates for such traits when such differences are present. Therefore, using 2 mult allows evolutionary rates for traits such as shape, coloration, thermal performance curves, ontogenetic growth trajectories, and other highly multivariate data sets to be examined and compared in a manner analogous to univariate data sets.
Second, deriving 2 mult from a distance-based perspective illustrates that the distance-covariance statistical equivalency also holds for analytical methods used for evolutionary questions. Typically in the biological sciences, this mathematical property is leveraged for analytical tools that assess patterns in contemporary data sets, and for methods that do not account for shared evolutionary history (for example, see Legendre and Legendre 1998;McArdle and Anderson 2001). By applying this equivalency in an evolutionary context, the approach developed here provides a framework for merging macroevolutionary theory with the analysis of high-dimensional data sets like shape. As such, it is possible that other evolutionary methods developed for univariate data may also be generalized for the analysis of high-dimensional phenotypic data.
Finally, the approach developed here could be additionally extended (in theory) to provide a framework for comparing evolutionary rates among traits that differ in their dimensionality. Currently, likelihood methods may be used to compare evolutionary rates between two or more traits on the same phylogeny (Adams 2013), but this approach is only implemented for univariate data. Unfortunately, many comparisons evolutionary biologists wish to make are 174 SYSTEMATIC BIOLOGY VOL. 63 between traits that could be quantified using a different number of dimensions (e.g., does body size evolve at a faster rate than body shape?). The methods developed here could be extended to address such questions, because the ratio between evolutionary rates for size and shape could be quantified, and then statistically evaluated via simulation. Thus, using methods based on 2 mult may provide a way forward in testing this, and other long-standing evolutionary hypotheses, so that we may better decipher how changes in the pace of evolution translate into large-scale patterns of phenotypic diversification within and among lineages. Under this model of evolution, the expected value at the root of the phylogeny (or phylogenetic mean) is estimated as a = (1 t C −1 1) −1 (1 t C −1 Y) where 1 is a vector of ones, and C and Y are as defined above (Rohlf 2001;Revell and Harmon 2008). In this case, a = 3.684211. The maximum-likelihood estimate of the evolutionary rate parameter under a BM model of evolution (O'Meara et al. 2006) may then be found using the standard variance-based equation: For this example, the estimate of 2 = 0.877193. To obtain the distance-based estimate ( 2 mult ), Y is first transformed by the phylogeny as where E(Y) is the phylogenetic mean and D is found from an eigen-decomposition of phylogenetic covariance matrix (Garland and Ives 2000;Blomberg et al. 2003). The distances between U Y and the origin are then obtained, which for this example are: Finally, 2 mult is obtained as 2 mult = PD t U,0 PD U,0 N .
For this example, the distance-based evolutionary rate is estimated as 2 mult = 0.877193, which identical to the univariate value obtained above.

APPENDIX 2
Computer Code for R The function below estimates the phylogenetic evolutionary rate parameter ( 2 mult ) for multidimensional traits for two or more groups as specified on a phylogeny. For two groups a test-ratio is obtained, while for more than two groups the average test-ratio is utilized as the test statistic. This value is then evaluated statistically using phylogenetic simulations based on a single evolutionary rate for all taxa in the phylogeny.