Evaluating the Predictive Power of Ordination Methods in Ecological Context

: When striving for the ordination methods best predicting independently measured site factors, the following questions arise: does the optimal choice depend on the kind of biological community analysed? Are there different ordination methods needed to address different site factors? Simultaneously, I explore alternative similarity approaches of entire ordinations, as well as the role of the transformations applied to the scale used in measuring species performance. The combination of methods and data transformations results in 96 alternative solutions for any one data set. These are compared by a graphical display, that is, an ordination of ordinations. The goodness-of-ﬁt of independently measured site factors is assessed by two alternative methods. The resulting 96 performance values serve as independent variables in trend surfaces overlaid to the ordination of ordinations. The results from two real-world data sets indicate that some ordination methods greatly vary with data transformation. Scores close to a binary scale perform best in almost all ordination methods. Methods that intrinsically constrain the range of species scores, such as principal components analysis based on correlation, correspondence analysis (including its detrended version), nonmetric multidimensional scaling, as well as principal coordinates analysis based on the Bray-Curtis distance, always ﬁgure among the most successful methods, irrespective of data used.


Introduction
Plant species and species groups tend to re-occur under similar natural conditions. Therefore, knowing the ecological factors modulating vegetation allows for the prediction of the composition of plant cover in probabilistic terms, and vice versa. Data analysis in vegetation ecology is, by and large, recognition of vegetation patterns and relating these to the same found in ecological space. The pattern is typically assessed in two steps: the first being a pairwise comparison of all sampling units to locate these in a similarity space; the second, searching for regularities therein [1].
A graph depicting multivariate observations as points in two-or three-dimensional dissimilar (distance) spaces is an ordination, which is the target of this paper. According to the similarity theory [2], interactions between living organisms and environmental conditions cause similar point patterns in their respective data spaces, a hallmark of past, present, or future ecological processes. For example, an environmental gradient from cold to warm conditions likely causes a parallel trend in species occurrence, from cold-tolerant to cold-intolerant ones. An ordination is considered useful if it allows the analyst to interpret emerging vegetation trends in environmental terms.
Finding good ordination methods was more elusive than expected [3]. None of the existing methods really fulfil all the theoretical requirements, as pointed out by Austin [4]. One would be to appropriately handle nonlinearity in species' response to the environment [5]. Given continuously varying environmental factors, species will perform best within an optimal range, decreasing above,

Ordination Methods Compared
The starting point of the analysis is a panel of ordinations, all derived from the same data set [9]. The evaluation focuses on some of the most popular ordination methods used in vegetation ecology, that is, PCA, PCOA, CA, DCA, and nonmetric multidimensional scaling (NMDS) (Figure 1, [3]). Two versions of PCA are involved in the evaluation. The first uses covariance to compare the species vectors without further transformation (1). The second is based on the correlation-based algorithm with intrinsic standardization ((2), see Supplementary Information S1).
Gower [15] devised a method to apply the principle of PCA to almost any measure of dissimilarity, whether metric or non-metric, named principal coordinates analysis (PCOA). I omit the version using Euclidean distance, as the outcome is the same as in PCA based on a variance-covariance matrix [7], but includes PCOA versions with Manhattan distance (3), Bray-Curtis distance (4), Canberra distance (5), correlation transformed to distance (6), and chord distance ( (7), see Supplementary Information S1).
Biologists also learned that they might interpret their data matrices as contingency tables for which Fisher's contingency table analysis [16] offered another alternative for ordination, known today as correspondence analysis CA (8), or its modified variant, detrended correspondence analysis (DCA (9), [17]). In CA, the species are considered as frequency counts even though they are usually of different origin (biomass records or cover percentage, for instance). DCA post-processes a CA ordination by altering the location of data points, that is, stretching horseshoe-shaped point clouds (p. 482, [8]). The implementation I use here is the same described by Oksanen and Minchin [18]. The result of DCA will depend on some of the input parameters chosen, or namely, the number of segments (I chose the default, which is 26 in the R implementation).
In NMDS [18], a given ordination is iteratively altered to minimize the rank order of distances between sampling units in the original (full) data space, and in the reduced one-, two-, or three-dimensional ordination. The algorithm presently used in R is described in more detail by Venables and Ripley [19] (p. 305). The versions included here are NMDS based on Euclidean distance (10), the same using Bray-Curtis distance (11), and NMDS based on correlation transformed to distance (12).

Scaling Species Performance
The proportion of the surface a species covers within a plot in the field is unstable and seasonally changing; therefore, plant ecologists often omit precise measurements, using course visual estimations instead. Quantities are frequently even disregarded, and only the presence or lack of the species is considered, as reflected when using a binary variable with zeros and ones. Van der Maarel [20] suggested the use of a rank scale for this estimation, with ranks from zero (species absent) to six (for species covering close to 100% of the plot). Using simple power transformation, x = x 2.5 this can be adjusted to mimic covering percentages on the one hand (0 ≤ x ≤ 100), and presence-absence on the other (0 ≤ x ≤ 1). As shown in Table 1, there is a convention to use symbols ("−", "+") for absence and occurrence in small quantities, respectively, followed by ranks for covering classes of uneven widths. Van der Maarel [20] suggested replacing these by a rank scale to simplify transformation and subsequent analysis. Even though Podani [21] argues that, from a formal point of view, species scores should be recognized as nominal-type data, the strategy of Maarel [20] is adopted here, as it has been used almost exclusively to date. Table 1. Transforming cover-abundance values in vegetation ecology. x 1.0 is the vector of rank scores replacing the conventional code in column one. x = x 2.5 is approximating the cover percentage [9,20].

Comparing Ordinations
A first possible method is the Mantel correlation [7,22]. In this Euclidean distance, matrices are computed from any two ordinations, and their corresponding elements are compared by correlation. Randomization of elements is used to assess error probability. This method of relating multiple sets of multivariate data was used by Ardisson et al. [23] to reveal the temporal pattern of epibenthic communities. As an alternative to this, Peres-Neto and Jackson [24] suggest the use of Procrustes superimposition introduced by Gower [25]. This aligns any two ordinations as far as is feasible by standardization of ordination scores, mirror reflection if required, and rotation, to minimize the squared differences between homologous points. In the present context, the sum of squared differences serves as a measure of ordination dissimilarity.
The ordinations of ordinations are obtained through PCOA by use of distance matrices comparing all individual ordinations. When using Mantel correlation, the elements have to be transformed according to: where r M is the Mantel correlation, and d M is the Mantel distance (Legendre and Legendre [8], p. 510). When using Procrustean analysis, no further transformation is required as the sum of squared differences are dissimilarities, and therefore this should be expanded upon PCOA.

Measuring Fit of Site Factors with Ordinations
Given a distance matrix derived from any one ordination (three-dimensional in the present case), a linear approach to measure fit with one site factor offers non-parametric multivariate analysis of variance (NP-MANOVA, [26]). This yields the proportion of the latter in explaining variance within the vegetation distance matrix, or as a percentage when multiplied by 100.
A more flexible method, although constrained to the first two ordination axes, is fitting a trend surface to the ordination. Various methods can do this, such as those compared by Venables and Ripley [19], for instance. One computer program suitable for the purpose is function ordisurf, found in R package vegan (version 3.4.3, R core development team 2017). This fits a generalized additive model (GAM) to obtain a smooth, non-linear surface. To measure the goodness of fit, I used Akaike's criterion (AIC) to express the residual deviance (Crawley [27], p. 208).
Irrespective of the method used, the fit of ordinations with a site factor is either given in tabular form or visualized as a trend surface across the ordination of ordinations.

Software Implementation
All computations shown below are implemented in a single R function, ordpanel(), included in version 2.1, and later of R package dave [9] (version 3.4.3, R core development team 2017) . The use of this function is shown in the Supplementary Information, S2. For ease of use, the function ordpanel offers four models, according to the methods used for the ordination of ordinations and the way the fit of site factors is measured: Model 1: Ordination by Mantel distance; fit measured by AIC criterion of the trend surface. Model 2: Ordination by Mantel distance; fit is variance explained in NP-MANOVA. Model 3: Ordination by sum of squares distances from Procrustes analysis; fit measured by AIC criterion of the trend surface. Model 4: Ordination by sum of squares distances from Procrustes analysis; fit is variance explained in NP-MANOVA.

Data
I present the results of two different case studies, explained in Wildi [9] and published in the related R package, dave 2.0 (version 3.4.3, R core development team 2017). Both rely on systematic plot sampling represented in a pair of multivariate data matrices, one assessing vegetation and the second a set of environmental factors. The first, named "sveg", describes a small wetland in the Pre-Alps, where the sample size is n = 63. The number of variables (species) in the vegetation matrix is p = 119, and in the environmental matrix, p = 20 (site factors). The main trend inherent in the data is caused by a difference in nutrients and water supply. The second, called "ws200", originates from a grid of Swiss forest samples with grid width of 4000 m. Sample size is n = 726, the number of species p = 1160, and the same of site factors p = 24. This data set is characterized by a strong gradient in elevation and factors related to this, like temperature and precipitation. Note that in both cases, the number of variables typically exceeds the same of sampling units, which, therefore, are over-determined.

Similarity and Performance of Ordinations
A graphical panel of ordinations which result when using data set "sveg" is shown in Figure 1, which is a slightly modified version of Figure 6.15 in Wildi [9]. The third dimension is visible too, as the circle diameter is taken proportional to the scores on axis number 3. The five groups of points highlighted by different colours result from cluster analysis, and serve as an ease of visual comparison of ordinations only. Figure 1 confirms some well-known rules mentioned above. First, there are ordinations strongly responding to the transformation of species scores, such as PCA using covariance in row (2), or NMDS using Euclidean distance in (10). In contrast to this, in row (5), the Canberra distance almost completely suppresses any rank difference, such that all ordinations in this row are quite similar.
On the one hand, ordinations tend to differ considerably when using species scores reflecting cover percentage (first column in Figure 1). On the other hand, when scores approach presence-absence, the ordination patterns tend to converge (last column in Figure 1).
Measuring pairwise dissimilarity of all ordinations in Figure 1 and assembling these to an ordination of ordinations results in Figure 2. Panel (a) is a full display of the Mantel distance-based ordination, with isoclines showing the fit of site factor "pH.peat" (a variable in data frame "ssit") to the individual ordinations of vegetation. The explicit values used for fitting the trend surface are given in Table 2. Symbols and colours distinguish ordination methods, and the transformation of data used within each method is expressed by symbol size (the largest corresponding to x = x 2.5 , and the smallest to x = x 0.01 ). Light lines connect identical methods. Panel (b) is an enlargement of the central area of panel (a). Panels (c) and (d) are ordinations in which ordination dissimilarity (distance) is measured as the sum of squared differences of analogue point locations obtained by procrustean analysis. Although they differ from panels (a) and (b), it is easy to see that the overall configurations are rather similar.
The content of Figure 2 can be summarized as follows: 1.
All methods relying on Euclidean or Manhattan distance, including PCA based on covariance (see Supplementary Information for details of these), strongly respond to transformation of species scores, as shown in Table 1. They extend over a wide range within all graphs.

2.
All methods using distance measures with intrinsic transformations (e.g., by using correlation which relies on standardisation of vectors), including the Bray-Curtis distance, Canberra distance, chord distance, and correlation used as a distance, but also PCA based on correlation, exhibit restricted response to species transformations. Their range in the ordination of ordinations is small. 3.
All ordinations tend to converge in terms of similarity patterns towards solutions based on presence-absence types of species scores (in the present case, x = x 0.01 ). These concentrate on the left-hand side of the graphs. 4.
The best fit to ordinations, irrespective of method, is achieved by transformations close to about x = x 0.5 , that is, approximately the square-root transformation of the rank scale, as proposed by van der Maarel [20], Table 1.

5.
Using the sum of squared distances of analogue point locations from procrustean analysis leads to ordinations of ordinations, and hence similar conclusions. 6.
Correspondence analysis (CA) is different from any other method in many ways, whereas its detrended version, DCA, is not.
Mathematics 2018, 6, x FOR PEER REVIEW 6 of 15 3. All ordinations tend to converge in terms of similarity patterns towards solutions based on presence-absence types of species scores (in the present case, x' = x 0.01 ). These concentrate on the left-hand side of the graphs. 4. The best fit to ordinations, irrespective of method, is achieved by transformations close to about x' = x 0.5 , that is, approximately the square-root transformation of the rank scale, as proposed by van der Maarel [20], Table 1. 5. Using the sum of squared distances of analogue point locations from procrustean analysis leads to ordinations of ordinations, and hence similar conclusions. 6. Correspondence analysis (CA) is different from any other method in many ways, whereas its detrended version, DCA, is not. Figure 1. The panel of ordinations obtained with the data frame "sveg", n = 63. Variation of ordination in vertical direction is caused by changing methods and in the horizontal direction by alternative data transformations. The classification of data points (colours) serves as an ease of visual comparison only, and has no effect on the results. Table 2. The efficiency of ordinations of "sveg" data ( Figure 1) measured by fitting trend surfaces of site factor "pH.peat" to the first two ordination axes. The elements of the table are Akaike's criterion (AIC) values expressing residual deviations with a correction for degrees of freedom [27]. Extensions in the method names signify: "cor" for correlation-based, "cov" for covariance-based, "m" for using and has no effect on the results. Table 2. The efficiency of ordinations of "sveg" data ( Figure 1) measured by fitting trend surfaces of site factor "pH.peat" to the first two ordination axes. The elements of the table are Akaike's criterion (AIC) values expressing residual deviations with a correction for degrees of freedom [27]. Extensions in the method names signify: "cor" for correlation-based, "cov" for covariance-based, "m" for using Manhattan distance, "e" for Euclidean distance, "bc" for Bray-Curtis distance, "chord" for chord distance, "can" for Canberra distance, and "cd" for correlation used as distance. Manhattan distance, "e" for Euclidean distance, "bc" for Bray-Curtis distance, "chord" for chord distance, "can" for Canberra distance, and "cd" for correlation used as distance.

Response to Different Site Factors
A site factor almost orthogonal to "pH.peat" in the "ssit" data is the average water level, "Waterlev.av", a proxy for water supply [9]. In all four ordinations of ordinations of Figure 3, trend surfaces are fitted to the efficiency scores of the individual ordinations in accounting for this site factor.
A comparison of Figure 2a, fitting pH of peat and Figure 3a, and fitting the average water level to the same ordination, suggests that different site factors in this specific case show the same overall pattern of agreement; that is, the fit of ordinations to site factors is best on the left-hand side where species performance scores approach presence-absence. Using NP-MANOVA instead to measure fit reveals a rather similar pattern (Figure 3b). In agreement with results from a variety of other methods [9], the water level response is generally weaker than the same of peat (Figure 2). This is a property inherent in the natural system analysed. As shown in Figure 3c,d the patterns differ when using sum of squared differences from Procrustes analysis for the computation of ordination of ordinations, but the overall trends in performance remain unchanged.   Figure 1. Light lines connect identical methods. Isoclines show the fit of site factor "Waterlev.av" (a variable in data-frame "ssit" measuring water supply) to the individual ordinations derived from vegetation data-frame "sveg". In panels (a,c) the fit of trend surfaces adjusted to the individual ordinations is used; in panels (b,d) it is variance explained by NP-MANOVA.
Surface fitting offers a quick interpretation of the ordination of ordinations, as shown in Figures  2 and 3. For a more detailed view, the explicit data for "pH.peat" is displayed in Table 2 (that is, the performance of ordinations as a function of method chosen and data transformation). These are displayed graphically in Figure 4a-d. Considering site factor "pH peat" and measuring the fit by trend surface analysis (panel (a)) shows CA being the best-performing ordination method, but when using NP-MANOVA instead, it is the worst (panel (b)), whereas DCA now performs best. Proceeding to site factor "average water level", panels (c) and (d) reveal that the same finding is valid here too, except that CA now always performs best. From the scaling of the y-axis, it can be seen that the fit of water level is poorer than the same of pH peat: the AIC values in panel (c) are  (a,b), and sum of squared differences obtained in procrustean analysis (c,d). Symbol size corresponds with the transformations shown in Figure 1. Light lines connect identical methods. Isoclines show the fit of site factor "Waterlev.av" (a variable in data-frame "ssit" measuring water supply) to the individual ordinations derived from vegetation data-frame "sveg". In panels (a,c) the fit of trend surfaces adjusted to the individual ordinations is used; in panels (b,d) it is variance explained by NP-MANOVA. Figures 2  and 3. For a more detailed view, the explicit data for "pH.peat" is displayed in Table 2 (that is, the performance of ordinations as a function of method chosen and data transformation). These are displayed graphically in Figure 4a-d. Considering site factor "pH peat" and measuring the fit by trend surface analysis (panel (a)) shows CA being the best-performing ordination method, but when using NP-MANOVA instead, it is the worst (panel (b)), whereas DCA now performs best. Proceeding to site factor "average water level", panels (c) and (d) reveal that the same finding is valid here too, except that CA now always performs best. From the scaling of the y-axis, it can be seen that the fit of water level is poorer than the same of pH peat: the AIC values in panel (c) are much higher than the same in panel (a), and the portion of explained variance in panel (d) is lower than the same in panel (b).  In the data sets "sveg" and "ssit", one gradient dominates the similarity pattern, resulting in a typical horseshoe-shaped point cloud [5]. In Figure 5, the same procedure is applied to the much larger data sets "ws200" and "wssit" to probe the findings in a case where more complex patterns can be expected. Unlike the panel in Figure 1, the point pattern is more difficult to interpret visually. Most striking is the tendency in PCA based on correlation (row (1)) and CA (row (8)), to mimic outlier observations-that is, data points located outside the main sample. This results from transforming vectors of rare species in which the majority of elements are zero scores, which is a known issue explained, for example, in Legendre and Legendre [8].  peat (a,b) and water level (c,d) in data sets "sveg" and "ssit". This performance is measured by an AIC value in the course of trend surface fitting (a,c), and alternatively, by explained variance in NP-MANOVA (b,d).

Surface fitting offers a quick interpretation of the ordination of ordinations, as shown in
In the data sets "sveg" and "ssit", one gradient dominates the similarity pattern, resulting in a typical horseshoe-shaped point cloud [5]. In Figure 5, the same procedure is applied to the much larger data sets "ws200" and "wssit" to probe the findings in a case where more complex patterns can be expected. Unlike the panel in Figure 1, the point pattern is more difficult to interpret visually.
Most striking is the tendency in PCA based on correlation (row (1)) and CA (row (8)), to mimic outlier observations-that is, data points located outside the main sample. This results from transforming vectors of rare species in which the majority of elements are zero scores, which is a known issue explained, for example, in Legendre and Legendre [8].
Mathematics 2018, 6, x FOR PEER REVIEW 11 of 15 Figure 5. The panel of ordinations, obtained with the data frame "ws200", n = 726. Variation of ordination in the vertical direction is caused by changing methods, and in the horizontal direction by alternative data transformations. The classification of data points (colours) serves as an ease of visual comparison only, and has no effect on the results. Figure 6 shows the analysis of the graphical display. Panel (a) is the fit of site factor elevation with all 96 ordinations, and panel (c) the same of site factor precipitation. The units of the y-axes are AIC values. These are smaller in panel (a) compared to (b), indicating that the fit of elevation is superior. CA and DCA outperform other methods even more clearly than in the smaller data set, "sveg". Panels (b) and (d) are the ordinations of ordinations based on Mantel distance. As in the previously analysed data (Figures 2 and 3), ordination of presence-absence-type scores are mainly concentrated on the left-hand side, but the overall similarity patterns differ widely from the wetland sample (Figures 1-4). CA and PCA based on correlation show unique patterns, whereas all types of NMDS are clustered, indicating that they are rather similar. Figure 5. The panel of ordinations, obtained with the data frame "ws200", n = 726. Variation of ordination in the vertical direction is caused by changing methods, and in the horizontal direction by alternative data transformations. The classification of data points (colours) serves as an ease of visual comparison only, and has no effect on the results. Figure 6 shows the analysis of the graphical display. Panel (a) is the fit of site factor elevation with all 96 ordinations, and panel (c) the same of site factor precipitation. The units of the y-axes are AIC values. These are smaller in panel (a) compared to (b), indicating that the fit of elevation is superior. CA and DCA outperform other methods even more clearly than in the smaller data set, "sveg". Panels (b) and (d) are the ordinations of ordinations based on Mantel distance. As in the previously analysed data (Figures 2 and 3), ordination of presence-absence-type scores are mainly concentrated on the left-hand side, but the overall similarity patterns differ widely from the wetland sample (Figures 1-4). CA and PCA based on correlation show unique patterns, whereas all types of NMDS are clustered, indicating that they are rather similar. Figure 6. Measuring the efficiency of ordinations in explaining elevation (a) and precipitation (c) in the data sets "ws200" and "wssit", respectively. The ordination of ordinations include trend surfaces from measuring the explanation power of elevation (b) and precipitation (d). Both are based on Mantel distance.

Discussion
The procedure devised in this paper to evaluate the performance of ordinations of vegetation data is based on independent site factors measured within the same ecosystem. Its application is limited to cases where are exist linear or non-linear relationships between the living part (plants, in this case) and the environment (site factors). This is usually the case in ecosystems dominated by one or several site factors-pH in the first example presented, and elevation in the second example. In contrast, this cannot be expected in highly disturbed systems carrying ruderal vegetation. The reliability of the results depends on the strength of this relationship. If it is weak, then random variation tends to blur interactions, and any ordination method will fail to predict environmental conditions. One of the consequences of this is that any evaluation of this kind will apply to the example under investigation only. The two case studies presented here confirm that the similarity pattern generated by alternative ordination methods (the ordinations of ordinations) is Figure 6. Measuring the efficiency of ordinations in explaining elevation (a) and precipitation (c) in the data sets "ws200" and "wssit", respectively. The ordination of ordinations include trend surfaces from measuring the explanation power of elevation (b) and precipitation (d). Both are based on Mantel distance.

Discussion
The procedure devised in this paper to evaluate the performance of ordinations of vegetation data is based on independent site factors measured within the same ecosystem. Its application is limited to cases where are exist linear or non-linear relationships between the living part (plants, in this case) and the environment (site factors). This is usually the case in ecosystems dominated by one or several site factors-pH in the first example presented, and elevation in the second example. In contrast, this cannot be expected in highly disturbed systems carrying ruderal vegetation. The reliability of the results depends on the strength of this relationship. If it is weak, then random variation tends to blur interactions, and any ordination method will fail to predict environmental conditions. One of the consequences of this is that any evaluation of this kind will apply to the example under investigation only. The two case studies presented here confirm that the similarity pattern generated by alternative ordination methods (the ordinations of ordinations) is case-dependent. Considering the performance pattern, however, there appears to exist some kind of tendency. In both data sets, correspondence analysis, as well as detrended correspondence analysis, frequently allows to explain a larger proportion of variance than is the case in other methods. This agrees with the findings of Ruokolainen and Salo [14], that: "In general, metric scaling methods, particularly CA and DCA, were far better in reflecting the main gradient in numerical terms, as compared with NMDS". It remains unclear whether this would be changed in cases where there is no dominant gradient, such as pH in the first example and elevation in the second one. Ruokolainen and Salo [14] also found NMDS to be pleasing "in graphical terms", and in the present evaluations, this group of methods always ranks among the best. In summary, they conclude that "none of the compared methods is perfect in reflecting a complex vegetation gradient".
The fact that the performance of ordination methods is case-dependent would explain why there is currently no final ranking described in the literature. Lepš and Šmilauer [28] concentrate on the use of program package CANOCO, and within this PCA and CA, both also providing constrained ordinations. In their introduction, they declare to abstain from the controversy between proponents of various approaches to analyse multivariate data. They simply argue that "the solutions presented are not the only ones, but they work for us, as well as many others". In light of the present evaluation, their choice of methods, PCA (based on correlation), CA, and DCA, is among the most favourable. Surprisingly, PCA, despite being a linear method, in Figure 2 (the "sveg" example) and Figure 6 ("ws200" example) is also among the better-performing methods.
For many vegetation ecologists, the overall performance of DCA will certainly come as a surprise. As explained in detail by Legendre and Legendre [8], DCA ordination is strictly based on the result of CA, but modified twofold. First, the location of data points is altered, such that any horseshoe-shaped point cloud is stretched (p. 482, [8]) to mimic a multivariate normal distribution. Whereas this is mainly considered a merely cosmetic manipulation, there is also a nonlinear rescaling of axes involved. This second modification apparently successfully reduces nonlinearity regularly found in vegetation data, and our methods measuring agreement of biotic and abiotic patterns signal an improvement. Although many ecologists use DCA as a standard method, superiority in comparison to other methods has never been shown [28]. Apparently, the most popular method is NMDS based on the Bray-Curtis distance, suggested in a paper by Clarke [29] with 7411 citations to date. In the evaluations of this paper, this also ranks among the best (Figure 4a-d).
Despite the completely different situation in the two case studies presented (in the wetland site, soil conditions differentiate vegetation; in the forests, climate-related factors dominate) the methodological conclusions are similar: transformation of species performance scores is frequently more influential than the choice of the ordination method. This is in agreement with findings of van der Maarel [20]. Further on that, all ordinations tend to converge in terms of similarity towards solutions based on presence-absence types of species scores-in the present case, x = x 1 . All methods relying on Euclidean or Manhattan distance, including PCA based on covariance, strongly respond to transformation of species scores, and the best fit with site factors is often achieved within the range of x = x 1 and x = x 0.5 .
The mathematically estimated fit of site factors with ordination pattern is an indication of the ease of visual interpretation. When using trend surface analysis, it is restricted to two dimensions, although the most important in search of patterns. Furthermore, the resulting isoclines regularly confirm nonlinear interactions. NP-MANOVA, on the other hand, allows for the handling of vegetation data in full dimensionality (by analysing a distance matrix) and relates this to one or several site factors [19]. However, it may fail to reveal an interaction in the case of strong nonlinearity of the vegetation pattern, as observed in Figure 4b for CA.
CA and PCA in all case studies proved to be surprisingly efficient in predicting site factors. However, as shown in Figure 5, when applied to larger samples there is an increasing tendency of generating aberrant observations (outliers) due to the intrinsic transformations. Lepš and Šmilauer [28], and also Legendre and Legendre [8] suggest circumventing the problem by suppressing rare species. However, our own tests showed that the problem cannot be solved in the case of even larger samples, posing a limit to the application spectrum of these methods.
The evaluations presented so far concern correlations between multivariate composition of vegetation and environmental factors. However, this is just one operation among others presently highly debated in vegetation ecology, sometimes addressed as a fourth-corner problem [9]. In this, the formal relationships among four data matrices is assessed: a plot described by a species abundance matrix, as used in this paper (1); a species by functional traits matrix (2); a plot by an environmental factors matrix (3); and a species' traits by environmental factors matrix (4) [29,30]. Hence, a strategy analogous to the one presented here could be chosen to evaluate the relationship between ordinations and species traits, or between environmental factors and traits. Kleyer et al. [31], in their comprehensive overview, present alternative methods for relating such data, where finding appropriate functions to overcome nonlinearity is crucial. Although they mention that species scores can either be included as a presence-absence binary code or by taking species abundance as weights, in view of the strong dependency of results from their transformation, all solutions to the fourth-corner problem may require a respective supplementary evaluation. The same may even apply to the search of optimal classification methods, such as that presented by Lengyel and Podani [32].

Conclusions
Ordination methods differ considerably in the power of predicting numeric environmental factors. The results suggest that a routinely performed evaluation is worthwhile, as the outcome is case-dependent. The method and software proposed allows for evaluating the best ordination methods individually in any data set describing vegetation and independently measured environmental factors. Whereas the method plays a crucial role in the efficiency of ordinations, the results confirm the importance of considering the transformation of species scores. In conclusion, even if the method choice follows some tradition or personal preference, the kind of testing proposed is worthwhile as it shows the potential of alternative ordination methods.