Multivariate comparison of variance in R

Comparisons of pattern and magnitude of phenotypic variation are central to many studies in evolution and ecology, but a meaningful comparison of multivariate variance patterns can be challenging. Here, we review an effective exploratory strategy, relative principal component analysis (relative PCA), for the comparison of variance–covariance matrices based on their relative eigenvalues and eigenvectors. Relative PCA allows for the identification of multivariate traits (linear combinations of variables) with maximal or minimal variance ratios between two groups. It can be used to explore the generation and canalization of phenotypic variance throughout ontogeny and phylogeny. Relative PCA also gives rise to a natural metric for the ordination of a sample of covariance matrices. We present a novel biometric justification of these approaches and discuss numerical difficulties as well as strategies for statistical inference, along with a new implementation of these methods in R. In an application of relative PCA to geometric morphometric data on cichlid body shape, we found that the phenotypic variance–covariance structure differs between males and females as well as between allopatric and sympatric populations of Tropheus moorii. Divergent selection in these populations mainly affected shape features related to swimming ability, whereas tropic morphology appears to be under stabilizing selection. In biology and biomedicine, individual variation is a key signal. Standard ordination methods as well as scalar summary statistics of multivariate variation often pool different biological factors with heterogeneous variational dynamics and thus conceal differences in variance–covariance pattern among groups. Relative PCA complements existing tests for the proportionality of covariance matrices and is an effective exploratory method to identify multivariate traits that differ in variance between populations, age groups, or treatment groups. By comparing within‐ and between‐population covariance matrices, relative PCA can reveal traits under divergent or stabilizing selection.

A biologically meaningful analysis of multivariate variance patterns is much more challenging than the analysis of averages.
Standard multivariate tools, such as principal component analysis, do not necessarily permit biological insights into the phenomena producing or canalizing phenotypic variation (see below and Figure 1). Scalar summary statistics of the magnitude of multivariate variation within a group (or of the difference in multivariate variation between groups) lump different factors with different biological dynamics; they rarely are biologically interpretable (Bookstein & Mitteroecker, 2014). For example, two groups, A, B, may have the same amount of total variance (sum of the variances for all variables), but the variance of the first variable may be higher in group A than in group B, and vice versa for the second variable. The finding of equal total variance would thus be misleading and conceal the different variational properties of the variables.
Here, we review an effective exploratory strategy, relative principal component analysis (relPCA), for the comparison of variance-covariance patterns based on their relative eigenvalues and eigenvectors (Bookstein & Mitteroecker, 2014;Flury, 1985). We present a new implementation of these methods in R, with an application to geometric morphometric data on cichlid body shape.

| MATERIAL S AND ME THODS
For a single continuous quantity x, the variation -or spread of the distribution -is commonly quantified by the sample variance the average squared deviation from the sample mean. Of course, two distributions can also differ in other properties than the mean and the variance: a distribution can be asymmetric (skewed), multimodal, lepto-or platykurtic (fat or thin tails), can have outliers, etc.
Meaningful comparison of variation across two or more samples in terms of the variance (Equation 1) requires that the distributions differ only (or at least primarily) in mean and variance, while other properties of the distributions are the same.
A comparison of variance across groups also requires that the variance can -in principle -change independently of the mean (or that the groups do not differ in mean values). In biology, however, the standard deviation (square root of the variance) of many size variables scales with the mean of that variable (large structures vary more than small ones), so that the variance needs to be interpreted relative to the mean. Such variables are thus often log-transformed because the resulting variance is unaffected by linear scaling. Other researchers prefer to divide the variable through its mean, which leads to the coefficient of variation (standard deviation expressed as a multiple of the mean) or its square (Houle, Pélabon, Wagner & Hansen, 2011). For small variations relative to the mean value, the mean-standardized variable is approximately proportional to the natural logarithm of this variable and lead to the same variance estimate. For larger ranges of variation, however, nonlinear transformations, such as the logarithm, severely distort the data space and affect most statistics (Huttegger & Mitteroecker, 2011). Ultimately, it depends on the biological and biometric context whether an association of mean and variance is considered an empirical result or an artefact that requires standardization of the variable.
Furthermore, not every measured variable has a biologically meaningful variance, i.e. it may not estimate an interpretable population property. For instance, consider a study of swimming speed in differently sized fish. If the fish were sampled randomly from a population, the sample variance in the length of the fish is an estimate of the population variance in length, which may have a biological interpretation. If the study design was based on an equal number of fish in, say, five different predefined size classes, the variance in fish length cannot be interpreted; there is no corresponding population property. Of course, the variance can still be computed, and a regression slope of swimming speed on fish length (which equals their covariance divided through the variance of fish length) may be interpretable, but the variance of fish length per se is not. Even if individuals are sampled randomly, variances may not necessarily be directly interpretable. For instance, if we sample fish of different ages, the variance of body length at a given age (a conditional variance) may be interpretable, but the variance of the full sample is not (because age may not have an interpretable variance). The same applies to samples of fossil specimens of different geological age.
Clearly, for a variance to be interpretable, the variable must also have at least an interval scale, enabling the meaningful computation of differences (Houle et al., 2011). On the other hand, a ratio scale, comprising a meaningful origin (zero value) of the variabe, is not necessary for computing the variance. But in some situations the (1) average deviation of individuals from the origin can be more insightful than the deviation from their mean. For instance, leg length has a ratio scale (no leg can be shorter than zero), but this origin is far outside of the actual sample distribution. Variation in leg length thus is measured by the variation around the sample mean (Equation 1).
But the difference in length between left and right leg, as a measure of asymmetry, has a meaningful origin, namely perfect symmetry, which is close to or within the sample distribution. Individual variation of leg asymmetry is measured best just by the mean sum of squares (i.e. the variance around the origin); the sample variance of leg asymmetry, often referred to as fluctuating asymmetry, neglects the average magnitude of asymmetry in the sample (directional asymmetry  (Bookstein, 2014;Houle et al., 2011;Huttegger & Mitteroecker, 2011). Even if all variables share the same unit, they (2) log Var(x B ) .
F I G U R E 1 Consider the two variables body height and body mass, measured for the three groups A, B, and C. Panel (a) shows a scatter plot of the log-transformed variances for both variables, measured in cm and kg, whereas panel (b) shows the log variances for the variables measured in m and g, instead. The change in unit (a linear transformation of the variables) has changed all variances, hence the different scales in the two diagrams, but the distances between the three points in the diagram have not changed. In other words, as differences between log-transformed variances are scale-invariant, the Euclidean distance in the space of log variances (Equation 4) is a scale-invariant metric for multiple variances. In (c), the variances are plotted on the original scales (without log transformation). A change of unit would alter this plot drastically, including the absolute and relative distances between the three points. For example, here the groups B and C appear to be more similar in variance one relative to the other than both are relative to group A. However, this statement depends on the -arbitraryunit of measurement. Expressing body mass in g instead of kg would lead to a different conclusion may not be considered equally 'important' in a specific biological context, especially if different signals in the data (e.g. different functional anatomical units) are covered by different numbers of variables. For a meaningful interpretation of a multivariate summary statistic, one may thus wish to weight these variables differently.
However, variance ratios are unaffected by any such linear weighting, differences of variances imply an (perhaps unwarranted) equal weighting of all variables (Huttegger & Mitteroecker, 2011).
The multivariate biometrics of variance should thus be based on variance ratios rather than differences between variances. However, many geometric concepts in multivariate statistics employ the notion of a 'distance' that gives rise to a metric data space or parameter space (Mitteroecker & Huttegger, 2009;Stadler, Stadler, Wagner & Fontana, 2001). The most common metric employed is the Euclidean distance, the square root of the summed squared differences for all measured variables. As the log variance ratio equals the difference of the logged variances, we may thus compare multiple groups on the basis of their log variances. In a space of log variances, the Euclidean distances between groups are scale-invariant, giving rise to a natural multivariate metric for variances: the square root of the summed squared differences between log variances, or equivalently, of the summed squared log variance ratios: where x i,A is the ith variable for group A, and p is the number of variables ( Figure 1).
The p variances are not sufficient to describe the multivariate variance-covariance pattern as they do not capture all the pairwise covariances between the variables. Only if the variables were all mutually uncorrelated, p variances would suffice for this purpose.
Similarly, the p variance ratios are not sufficient to describe all the differences between two multivariate variance patterns: only if the variables were mutually uncorrelated within both groups (i.e. all variables are uncorrelated among the specimens of group A and also among those of group B, but not necessarily in the pooled sample), p variance ratios would suffice to describe all the differences in the multivariate variance patterns.
It is common in multivariate statistics to transform the measured variables into a set of one or more linear combinations (weighted sums) of these variables. Geometrically, this equals a projection of the data points onto one or more new coordinate axes in the data space.
It is a fundamental theorem that every set of variables can be transformed linearly into another set of variables that are all uncorrelated and for which the coordinate axes remain orthogonal: these variables are called principal components (PCs; e.g. Jolliffe, 2002;Mardia, Kent & Bibby, 1979;Mitteroecker & Bookstein, 2011). On the basis of these principal components, the variance pattern can be described sufficiently by the p variances for these components; all covariances between them are zero. The orthogonality of the PC axes implies that, geometrically, principal component analysis is a mere rotation of the data space and thus leaves the variance pattern unchanged.
Although less familiar than principal component analysis, it is also possible to transform a set of variables into a set of linear combinations that are uncorrelated separately within two groups, the so-called relative principal components (relative PCs; Bookstein & Mitteroecker, 2014;Flury, 1985). Hence, on the basis of these relative PCs, the difference in multivariate variance-covariance pattern between two samples can be described sufficiently by the corresponding p variance ratios. In contrast to ordinary PCA, however, relative principal component analysis is a non-orthogonal transformation that alters each variance but leaves the variance ratios unchanged (as they are invariant to linear transformation; see Equation 2).
Let X A , X B be the mean-centred data matrices of groups A and B with dimensions n A × p and n B × p, respectively. The p×p covariance ma- thus are symmetric (i.e. they remain the same if exchanging group A with group B, cf. Equation 3) and give rise to a fully multivariate metric for variance-covariance matrices: where i are the eigenvalues of either S −1 This metric is equivalent to the Fisher information metric for multivariate normal distributions with common mean vectors. It is also the Riemannian metric on the space of square symmetric positive definite matrices, a manifold of dimension p(p + 1)/2 with the form of a convex cone in the vector space of symmetric matrices (Dryden, Koloydenko & Zhou, 2009;Förstner & Moonen, 1999;Lenglet, Rousson, Deriche & Faugeras, 2006;Mitteroecker & Bookstein, 2009).
Relative PCA and the metric in Equation 5 apply only to two groups with distinct variance-covariance matrices. For a set of m covariance matrices, the m(m − 1)/2 pairwise Riemannian distances (Equation 5) can be used to quantify the heterogeneity of all the variance-covariance patterns. Because these distances span a highdimensional non-Euclidean space, they cannot be directly plotted, but they can be represented by a principal coordinate analysis (also referred to as metric multidimensional scaling), which approximates (in a least squares sense) the Riemannian distances by the Euclidean distances in a low-dimensional diagram (Mitteroecker & Bookstein, 2009).
Relative eigenvalues can also be used to express deviations in the magnitude of multivariate variation (rather than in the full variancecovariance pattern) between two groups. The product of relative eigenvalues equals the ratio of the generalized variances (determinant of the covariance matrix) of the two groups (Mitteroecker & Bookstein, 2009): As relative eigenvalues are invariant to linear scaling of the variables, also ratios of generalized variances are unaffected by the scale of variables (Huttegger & Mitteroecker, 2011). Plots of log generalized variances are thus useful to compare the magnitude of multivariate variation across groups.

| Numerical aspects
The computation of relative eigenvalues requires the inversion of a covariance matrix, which imposes a certain computational burden: only matrices of full-rank can be inverted. For a covariance matrix of full rank, all principal components have non-zero variance, which requires an excess of cases over variables as well as linear independence of the variables. Linear dependences typically arise in the course of standardization of variables. Procrustes shape coordinates, for instance, are standardized for position, scale, orientation. Due to the resulting loss in degrees of freedom, the last four PCs (for 2D landmarks) or seven PCs (for 3D landmarks) have zero variance, and the corresponding covariance matrix cannot be inverted (Dryden & Mardia, 1998;Mitteroecker & Gunz, 2009).
Even for full-rank data, the relative eigenvalues are affected by the number of variables (p) relative to the number of cases (n). The larger the ratio p/n, the larger are the first relative eigenvalues and the more unstable are the corresponding eigenvectors (Bookstein, 2017;Mitteroecker & Bookstein, 2011).

| Euclidean approximation
Another way to circumvent numerical challenges is to replace the  Figure 1).

| Statistical inference
Like ordinary PCA, relative PCA is primarily an exploratory technique that is applied when no specific hypotheses could have been formulated a priori. As in related multivariate settings, it is useful to start with a test of the proportionality of covariance matrices: This test is equivalent to a test of the equality of all the relative eigenvalues of S A and S B : For two multivariate normal distributions, the corresponding maximum likelihood (ML) test is based on the quantity The parameters a and g are the arithmetic and geometric means of the relative eigenvalues, and n is the sample size per group. In the case of unequal sample sizes, n can be estimated by the harmonic mean of the two sample sizes. In the limit of n ≫ p, the log likelihood ratio (Equation 10) is distributed approximately as a 2 with (p − 1) (p + 2)/2 degrees of freedom, where p is the number of eigenvalues being compared (Anderson, 1963;Bookstein & Mitteroecker, 2014).
The arithmetic mean of the relative eigenvalues also serves as the maximum likelihood estimate of the scaling factor (the parameter k in Equation 8) for two proportional covariance matrices (Mardia et al., 1979).
After rejecting the hypothesis of proportionality, relative PCA can be used for identifying the linear combinations that drive this deviation from proportionality. In some situations, it can be interesting to test if the first (or last) two relative eigenvalues differ and thus support a separate interpretation of the corresponding relative eigenvectors. In this case, the test statistic in Equation 10 reduces to (8) H 0 :S B = kS A .
F I G U R E 2 (a) Consider two groups, A and B, and two variables, x 1 , x 2 , with covariance matrices visualized by the solid and dashed equal frequency ellipses, respectively. The dashed lines represent the first principal component (PC) for each group, i.e. the major axes of the two ellipses. The solid line is the first PC of the pooled data (assuming equal sample size), hence ignoring the group structure in the data. (b) The axes are rotated to the PCs of the pooled data. The resulting PC scores are uncorrelated for the pooled data but not within each group (the two ellipses are not oriented along their major axes). (c) When rotated to the PCs of group A, the scores are uncorrelated in A but not in group B (and vice versa if rotated to the PCs of B). The two variance patterns thus cannot be compared sufficiently by the variance ratios for any of these PCs because they would ignore the covariances between the scores. (d) The lines here represent the relative PCs (the directions with maximal and minimal variance in group A relative to B). In contrast to ordinary PCs, the relative PCs are not orthogonal. (e) Projecting the data on these relative PC axes (corresponding to a non-orthogonal rotation) yields scores that are uncorrelated within both groups: the two ellipses are aligned along their major and minor axes. On this basis, the two variance ratios sufficiently describe deviations in variancecovariance structure between the groups

| PACK AG E DE SCRIP TI ON AND I N S TA LL ATI O N
The vcvComp package is written in the R scientific computing language (R Development Core Team, 2018). It provides functions for the comparison of variance-covariance patterns based on relative eigenanalysis ( and T. polli, collected from eight locations of Lake Tanganyika (Herler, Kerschbaumer, Mitteroecker, Postl & Sturmbauer, 2010;Kerschbaumer, Mitteroecker & Sturmbauer, 2014). The data are also available via the Dryad repository (Kerschbaumer, Mitteroecker & Sturmbauer, 2013). we investigated if and how they differ in phenotypic variance-covariance structure. We also explored differences in variance pattern between female and male specimens, because these populations show significant sexual dimorphism in mean head shape (cichlids are maternal mouthbrooders; Herler et al., 2010;Kerschbaumer et al., 2014). Finally, we searched for signs of stabilizing and divergent selection among the six Tropheus populations by contrasting withinand between-group covariance matrices.
First, we loaded the vcvComp package and the data.
Five specimens are outliers for landmark 2 and were excluded from the sample. After selecting the subsample, we created a new variable combining population and sex. scaling.BW Computes the ML scaling factor between two covariance matrices (k in Equation 8) The landmark coordinates were extracted to create a matrix.
Then, we performed a generalized Procrustes superimposition (Rohlf & Slice, 1990) of the landmark coordinates using the function gpagen of the geomorph package.
We reduced the Procrustes shape coordinates to the first five principal components to avoid collinearities and to guarantee a sufficient excess of cases over variables in further analyses.

| Population comparison
Because the samples were not balanced regarding sex, we computed the pooled population covariance matrices as unweighted averages of the corresponding male and female covariance matrices.
To explore the heterogeneity of variance-covariance structure in body shape across populations, we performed an ordination analysis of the six pooled within-sex covariance matrices. F I G U R E 4 Eigenanalysis of the population IKA1 relative to IKS5 (relative PCA). (a) Relative eigenvalues (on a log scale), which are equal to the ratios of variance for the corresponding relative PCs. (b) Visualization as thin-plate spline (TPS) deformation grids (Bookstein, 1991) of the shape pattern corresponding to the first relative PC, which has the maximal excess of variance in IKA1 relative to IKS5

(a) (b)
The first three principal coordinates together accounted for 88% of total variance (Figure 3a); this also equals the fraction of summed squared Riemannian distances explained by the summed squared Euclidean distances within the first three principal coordinates. The populations living in sympatry (IKS3, IKS4, IKS5) were separated from the allopatric populations (IKA1, IKA2, IKA3) along the third principal coordinate (Figure 3b).
To investigate the actual differences in variance-covariance pattern between sympatric and allopatric populations, we compared the populations IKA1 and IKS5, but other pairs of sympatric and allopatric populations led to very similar results. The ML test (Equation   10) indicated that the covariance matrices of IKA1 and IKS5 deviate from proportionality at p = .02.
The generalized variance of IKA1 was only 18% less than that of IKS5, but the relative PCA showed that the various shape features deviate strongly in their variational properties across populations ( Figure 4a). The first relative PC was roughly twice as variables in IKA1 than in IKS5 (first relative eigenvalue was 2.3), whereas the variance of the last relative PC in IKA1 was only half of that in IKS5 (last relative eigenvalue was 0.49).
The shape patterns depicted by each relative eigenvector can be visualized by deformations of the average shape along the positive and the negative directions of the corresponding vector ( Figure 4b). Note that when the initial variables were reduced to the first principal components, as was the case here, the loadings of the eigenvectors must be multiplied by the loadings of the principal components to get shape patterns.
The shape features captured by relative PC 1 were head shape, relative eye size, and body depth (maximum distance between dorsal and ventral parts); these were the features with maximal excess of variance in allopatric populations relative to sympatric populations ( Figure 4b). Or in other words, these were the shape features maximally canalized in the populations living in sympatry.
In Lake Tanganyika, allopatric populations of Tropheus moorii live in the whole water column, whereas populations in sympatry with T.
polli usually are forced to live at greater water depth (Kerschbaumer et al., 2014). The broader trophic niche and the larger environmental heterogeneity in allopatric populations may account for the larger variance in body depth and head shape, but the higher competition and harder living conditions in sympatric populations may also impose a stronger stabilizing selection regime than in allopatric populations.

| Comparison between sexes
We separated males and females and performed a principal coordinate analysis of the 12 sex-specific covariance matrices in order to investigate deviations in variance-covariance structure between the sexes.
The first two components together accounted for 62% of total variance ( Figure 5) and showed that for all populations, except IKA3, males had higher values than females for the first principal coordinate ( Figure 6).
To explore this shared sex difference in variance-covariance pattern, we performed a relative PCA of the females relative to the males of IKA1.
Overall, females were approximately half as variable as males (ratio of generalized eigenvalues was 0.57), but pooling over all dimensions was again misleading here. In fact, the first relative PC was 4.6 times more variable in females than in males, whereas the other dimensions were all more variable in males ( Figure 8).
The first relative PC mainly corresponded to the relative size of the head (Figure 8a), whereas the last three relative PCs were all related to the shape and orientation of the head and mouth ( Figure 8b).
F I G U R E 5 Fraction of variance explained by each principal coordinate of the 12 sex-specific covariance matrices Cichlids are mouth brooders and females typically have a larger head and mouth than males. This pattern of sexual dimorphism in body shape was also found for the present Tropheus moorii sample (Herler et al., 2010;Kerschbaumer et al., 2014). The increased variance in relative head size likely is a direct consequence of the enlarged head in females, whereas other aspects of head morphology, such as the relative position and orientation of the mouth, seems to be more canalized in females than in males.

| Stabilizing versus divergent selection of cichlid body shape
Under idealized assumptions, the expected amount of phenotypic change due to genetic drift is proportional to the amount of additive genetic variance within the population. Extending this model of neutral evolution to multiple traits leads to the expectation that the between-group covariance matrix for a set of related populations is proportional to the additive genetic covariance matrix of their common ancestral population (Lande, 1979). This rational has inspired statistical tests for natural selection by contrasting the covariance matrix of population means with the pooled phenotypic withinpopulation covariance matrix (as a surrogate of the ancestral genetic covariance matrix; e.g. Ackermann & Cheverud, 2004;Cheverud, 1988;Marroig & Cheverud, 2004;Martin, Chapuis & Goudet, 2008): deviations from proportionality are signs of stabilizing or divergent selection. Most of these approaches, however, only rely on statistical significance tests of proportionality of the between-and withinpopulation covariance matrices (B and W, respectively). Relative PCA ideally complements these approaches as an exploratory tool to identify the specific trait combinations that deviate from the null model of neutral evolution (Bookstein & Mitteroecker, 2014). If both divergent and stabilizing selection acted in a set of populations, the first relative PCs of B with respect to W will reveal the trait combinations that were affected by divergent selection (the features with maximal between-population variance relative to within-population variance), and the last relative PCs will show the trait combinations under stabilizing selection (least between-population variance relative to within-population variance).
We computed B and W (pooled by sex) for the Tropheus populations based on the first five PCs of the Procrustes coordinates. As we have only six populations in this example, B is estimated with great uncertainty and also the chi-square approximation of equations 10 and 11 is critical; results have to be interpreted with care.
The ML test suggested a deviation from proportionality between B and W (p=.034) and thus the action of selective forces.
F I G U R E 6 Principal coordinates ordination of the 12 sex-specific covariance matrices. Males in blue, females in red. Populations living in sympatry with Tropheus polli in dark colors F I G U R E 7 Relative eigenvalues (maximal ratios of variance) of females relative to males for the population IKA1 F I G U R E 8 Visualization as thin-plate spline (TPS) deformation grids (Bookstein, 1991) of the shape patterns corresponding to (a) the first relative PC, which has the maximal excess of variance in females relative to males for the population IKA1, and (b) the last relative PC, which has the maximal excess of variance in males relative to females However, this test does not specify the magnitude of deviation from proportionality, and especially with the small number of populations in this example, the interpretation of the p-value alone is not sufficient. We thus performed an ordination of the six population covariance matrices (pooled by sex), together with W and B (scaled to fit W using the mean of their relative eigenvalues). Figure 9 shows that relative to the heterogeneity of population covariance matrices, B clearly deviates from W along the first principal coordinate. The interpretation of the relative PCs of B and W thus seems warranted. Figure 10 shows the relative eigenvalues of B with respect to W. The first relative eigenvalue is more than 5 times larger than the second one, which is significant at p < .05, and similarly for the last relative eigenvalue (Equation 11). This supports an interpretation of the first and last relative PC, even though the absolute relative eigenvalues cannot be evolutionarily interpreted without knowing the variance ratio expected under neutral evolution, which depends on the number of generations since divergence and effective population size (Lande, 1979). One way to estimate this expected variance ratio is based on genetic data using the F ST statistic (Holsinger & Weir, 2009). Under pure genetic drift, Lynch & Walsh, 1998;Martin et al., 2008). Kerschbaumer et al. (2014) reported F ST values for these F I G U R E 9 Principal coordinates ordination of the six populations (males and females pooled), along with their between-group (B) and their within-group (W) covariance matrices F I G U R E 1 0 Relative eigenvalues of the between-group covariance matrix versus the within-group covariance matrix for the six Tropheus populations F I G U R E 11 Visualization of the shape patterns corresponding to (a) the first relative PC, which has the maximal excess of variance between the six Tropheus populations relative to the variance within populations, and (b) the last relative PC, which has the maximal excess of variance within populations relative to that between populations populations ranging from 0.033 to 0.085, which translates into ratios of between-to within-population variance of 0.034-0.093.
The first relative eigenvalue (2.06) clearly exceeded this threshold and suggests strong divergent selection. The last relative eigenvalue (0.025) may indicate weak stabilizing selection.
The first relative PC corresponded to overall body depth and the positions of caudal, dorsal and ventral fins (Figure 11a). These features, which determine the hydrodynamics and swimming ability of the fish, were under strong divergent selection in the studied Tropheus populations: their inter-population variance strongly exceeded the variation expected for neutral evolution.
The shape pattern corresponding to the last relative eigenvector mainly involved the shape of the head, especially the position and orientation of the mouth (Figure 11b). These features, which crucially affect feeding performance, likely were under stabilizing selection. Canalization of trophic head morphology is also supported by our finding that females show little variation in mouth position and orientation despite their large and variable head size (Figure 8).

| CITATI ON OF VC VCO MP
Scientists using vcvComp in a published paper should cite this article or the vcvComp package directly. Citation information can be obtained by typing:

DATA AVA I L A B I L I T Y S TAT E M E N T
The example dataset used in this manuscript is packaged with vcv-Comp and is publicly available via CRAN (https ://CRAN.R-proje ct.org/packa ge=vcvComp). The data are also available via the Dryad repository (Kerschbaumer et al., 2013) https ://doi.org/10.5061/ dryad.fc02f.