Evaluating the Predictive Power of Ordination Methods in Ecological Context

Abstract

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-fit 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 figure among the most successful methods, irrespective of data used.

1. 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, as well as below this. The result is an approximately Gaussian-shaped response curve. A second complication is the presence of a high proportion of statistical noise occurring in nature, due to fluctuations of ecosystems, such as random variations of meteorological conditions, heterogeneity of soils, and disturbances arising from human influence. A third problem is the high dimension of resemblance spaces, as vegetation composition is frequently described by hundreds of variables. As a result, the number of ordination methods to choose from is large, and as von Wehrden et al. [3] mentioned, hardly any of these have recently disappeared, while their continued use is mainly driven by tradition, rather than experience or evaluation.

An ordination is a convenient tool used to visualize patterns in a similarity space, such as trends, groups, outliers, or combinations of these. The classical method used to achieve such graphs is principal component analysis (PCA, [6,7,8,9]). PCA centres the original variables (here: the species) and replaces the original variables by linear combinations, maximizing the total variation of observations. This is nothing else here but rotating the point cloud to the end of maximizing insights into the dominant pattern. Reducing the resulting data matrix to display the first two or three axes with maximum variation simplifies visual inspection ([10]). PCA assumes linearity of interactions, whereas in ecology, nonlinearity is expected to occur. Plant ecologists therefore continue to use various alternative ordination methods. Von Wehrden et al. ([3]) presented an overview of these and recent trends in application, concluding that, “Ongoing discussions and the diversity of available methods have led to considerable uncertainty among researchers and students”. They found that the choice is subject to permanent change, fostering the impression that ordination methods, as yet, lack theoretical foundation. Or, as Austin ([4]) concludes, “There seems to be an inconsistency between theory and methodology”. However, an inconsistency only exists under the assumption that there must be a formal (mathematical) relationship between environmental factors and the composition of vegetation. According to the similarity theory [2], environmental and biological patterns tend to converge, but only to a degree limited by natural disturbance and other forces, like competition, causing divergence. In this paper, I hypothesize that the degree of convergence is dependent on the ecosystem investigated. I therefore propose a software tool to probe for the best ordination method in any single ecosystem, rather than relying on traditions or past experience.

As the number of methods increased, papers comparing these became abundant (e.g., [11,12]). Many authors used simulated vegetation data in the hope of better understanding their interaction with artificial environmental gradients. Minchin ([13]) recognized the need for an extended systematic comparison of methods by measuring the performance of ordinations in some way. His evaluation was entirely based on artificial data, which are highly controlled, but hardly able to mimic all the complications occurring in natural situations. Ruokolainen and Salo ([14]), in their comparison of four ordination methods, concluded that simulations were too simplistic to be able to predict the outcome with real data. They therefore used independently-recoded fire intensity as a reference-site factor to determine the predictive power of vegetation ordinations. Accordingly, in the present paper I evaluate as many as 96 alternative ordinations of one and the same data set simultaneously, by measuring performance through predicting site factors in two real-world examples. This addresses the question on whether the optimum choice of method differs with the kind of vegetation analysed (1), as well as the question on whether, within a given system, there are different ordination methods needed to predict different site factors, such as climatic variables vs. soil parameters, for example (2). As will be shown in the Methods section, additional technical issues emerge from the inclusion of a large number of ordinations, like choosing a method to assess the similarity of entire ordinations, or how to relate vegetation ordinations and site factors.

2. Methods

2.1. 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).

2.2. 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 $\left(0\le x\le 100\right)$, and presence-absence on the other $\left(0\le x\le 1\right)$. As shown in

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^{\prime }=x^{2.5}$ is approximating the cover percentage [9,20].

Code	Cover (%)	$\mathbf{Rank},{\mathit{x}}^{1.0}$	${\mathit{x}}^{0.1}$	${\mathit{x}}^{0.25}$	${\mathit{x}}^{0.5}$	${\mathit{x}}^{2.5}$	${\mathit{x}}^{4.0}$
−	0	0	0	0	0	0	0
+	<1	1	1	1	1	1	1
1	5	2	1.07	1.19	1.41	5.65	16
2	17.5	3	1.12	1.31	1.73	15.58	81
3	37.5	4	1.15	1.41	2.00	32.00	256
4	62.5	5	1.17	1.50	2.24	55.90	625
5	87.5	6	1.19	1.57	2.45	88.18	1296

, 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.

2.3. 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:

$$d_M=\left(1-r_M\right)/2$$

(1)

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.

2.4. 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.

2.5. 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.

3. 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.

4. Results

4.1. 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. 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.

	x′ = x2.5	x′ = x2	x′ = x1.5	x′ = x1.25	x′ = x1	x′ = x0.75	x′ = x0.5	x′ = x0.25	x′ = x0.01
PCA.cor	65.3	64.1	61.9	60.6	59.4	58.7	58.7	59.4	61.0
PCA.cov	176.8	142.1	88.9	76.0	67.5	62.9	60.6	61.0	62.7
PCoA.m	80.3	73.5	68.7	66.8	65.8	65.3	65.5	66.3	67.2
PCoA.bc	64.6	56.0	53.8	53.9	54.4	55.1	55.9	56.7	57.5
PCoA.chord	125.0	118.4	79.0	68.4	59.7	55.2	54.1	54.9	56.3
PCoA.can	59.6	58.9	58.2	57.8	57.4	57.1	56.7	56.5	56.3
PCoA.cd	121.1	93.1	77.3	67.7	61.0	58.4	58.5	59.8	61.5
CA	70.4	57.2	48.3	47.2	47.9	49.3	50.9	52.5	53.9
DCA	75.9	57.4	49.2	48.7	49.1	50.2	51.5	52.6	53.6
NMDS.e	189.4	173.7	90.0	80.3	74.1	66.9	67.6	67.5	67.4
NMDS.bc	65.3	60.7	58.8	57.7	57.7	58.8	59.2	59.9	59.9
NMDS.cd	99.1	77.9	84.5	64.4	60.5	60.8	63.2	64.9	65.7

. 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:

• 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.
• 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.
• 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.
• 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.
• Using the sum of squared distances of analogue point locations from procrustean analysis leads to ordinations of ordinations, and hence similar conclusions.
• Correspondence analysis (CA) is different from any other method in many ways, whereas its detrended version, DCA, is not.

4.2. 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.

Surface fitting offers a quick interpretation of the ordination of ordinations, as shown in Figure 2 and Figure 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].

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 (Figure 2 and Figure 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 (Figure 1, Figure 2, Figure 3 and Figure 4). CA and PCA based on correlation show unique patterns, whereas all types of NMDS are clustered, indicating that they are rather similar.

5. 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¹. 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¹ 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].

6. 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.