AMQ 31(1) 59-73 DiDonato et al ProofCopy Available online http://amq.aiqua.it ISSN (print): 2279-7327, ISSN (online): 2279-7335 Alpine and Mediterranean Quaternary, 31 (1), 2018, 59 - 73 PALAEOENVIRONMENTAL RECONSTRUCTIONS THROUGH COMPOSITIONAL DATA ANALYSIS Valentino Di Donato1, Josep Antoni Martín-Fernández 2, Marc Comas-Cufí 2, Joanna Jamka1 1 Dipartimento di Scienze della Terra, dell’Ambiente e delle Risorse, Università degli Studi di Napoli “Federico II”, Naples, Italy. 2 Departament d’Informàtica, Matemàtica Aplicada i Estadística, Universitat de Girona, Girona, Spain. Corresponding author: V. Di Donato ABSTRACT: The modern analogue technique has been accordingly revised with compositional data analysis framework. The method adopts the Aitchison distance, obtained from isometric log-ratio coordinates of relative abundances, as a natural measure of similarity among assemblages. The number of analogues from which obtain the estimates was determined through leave-one- out verification of modern assemblages. Mean distances and local outlier factor are considered to evaluate the quality of palaeoes- timates. The method has been tested on Atlantic Ocean and Mediterranean planktonic foraminiferal assemblages to reconstruct past sea surface temperatures (SST). In comparison with previous planktonic foraminiferal based reconstructions, the Portugal offshore SST record for the last 210 ka shows a higher coherence with other paleoclimate proxies (i.e. the stable isotope and al- kenone records). In comparison with raw data analysis, CoDa-MAT yields lower estimates of the western Mediterranean last gla- cial period SST. Keywords: modern analogue technique, Aitchison distance, sea surface temperatures Supplementary Appendices only online at http://amq.aiqua.it 1. INTRODUCTION Quantitative estimation of past environmental pa- rameters is one of the most challenging engagements of palaeoclimatical and palaeoecological investigations. In palaeoceanographical studies quantitative palaeocli- matic reconstructions can be attained by means of geo- chemical methods, such as the analysis of long chain alkenones (Brassel et al., 1986) and the Mg/Ca ratio in foraminifera (Elderfield & Ganssen, 2000), or more strictly paleontological methods based of statistical analysis of census counts of assemblages. These meth- ods, known as transfer function methods, have been introduced to provide quantitative estimates from counts of fossils, including regression-based methods (Imbrie & Kipp, 1971), modern analogue techniques (MAT) (Hutson, 1979; Pflaumann et al., 1996; Waelbroeck et al., 1998) and artificial neural networks (Malmgren et al., 2001). The starting point of these methods is a modern dataset usually consisting of counts of modern assem- blages and measured environmental parameters. In this paper we focus on the MAT which compares fossil as- semblages with modern ones using a distance measure or a similarity coefficient. The palaeoenvironmental esti- mates are obtained from the environmental parameters measured at the location of the most similar modern assemblages. For each fossil samples the nearest mod- ern ones are found by adopting an appropriate distance (d). Then, the palaeoestimate p̂ for a fossil is repre- sented by the mean of the environmental parameters pi measured at the geographical location of the h modern analogues: where h represents the number of analogues. Following Hutson (1979) and Pfaumann et al. (1996), the mean can be weighted on the inverse of the distance (1/d(i)), so as increase the influence of the closer analogues on the palaeoestimate: ( ) ( ) 1 ( ) 1 (1/ ) ˆ 1 (1/ ) h i ii h i i d p p d = = ⋅ =   The assessment of the goodness of the method is usually carried out through cross-validation in the refer- ence modern dataset. Following the discussion along the works Telford & Birks (2009, 2011) and Guiot & de Vernal (2011a, 2011b), we assume that in MAT the spa- tial structure of data has relatively low effect on the cal- culation of prediction errors. In consequence, in this work we assume no effect of spatial autocorrelation. MAT mostly differ for the distance measures or similarity indexes used, among which the cosine-theta (Hutson, 1979), the scalar product of the normalized https://doi.org/10.26382/AMQ.2018.05 60 Di Donato V. et al. assemblages vectors (Pflaumann et al., 1996), squared chord distance (Waelbroeck et al., 1998) and the Euclid- ean distance on the logarithm of the species relative abundances in permil (Guiot & de Vernal, 2011a). The information in a modern and a fossil assemblage is the relative abundance or percentages of species, i.e., they can be considered as compositional data (CoDa) (Aitchison, 1986). That is, the information contained in a vector of counts x is the same as in k·x, for any real scalar k>0, property known as scale invariance (Aitchison, 1986). That is, any observation x is a mem- ber of an equivalence class (Barceló-Vidal & Martín- Fernández, 2016). This type of data is common in Earth Sciences when the constituents and compounds are described in terms their concentration (e.g., Buccianti et al., 2006). As with some of the cluster analysis tech- niques, despite it is also possible to experiment with different distances, we follow Everitt et al. (2011, p. 69) recommendation “…the choice of measure will be guided largely by the type of variables being used and the intuition of the investigator”. In the particular case of CoDa the Aitchison distance has been proved to be an appropriate measure for the geometry of the sample space (Aitchison et al., 2000). Palarea-Albaladejo et al. (2012) present a summary of the properties, advantages and difficulties of several different measures when they are used for CoDa. The approach proposed on this paper has the aim to develop a MAT in a coherent fashion with the nature of the data, hereafter the CoDa-MAT. In the following section, the basic elements of the log-ratio methodology to CoDa are introduced and the CoDa-MAT is de- scribed. Next the method is tested, and its results illus- trated using different assemblage datasets. Finally, Sec- tion 4 presents some conclusions and final remarks. 2. THE CoDa-MAT CoDa refers to vectors of positive components showing the relative weight of a set of parts in a total. Nowadays, there is a general agreement that applying the standard statistical methods to CoDa may yield mis- leading results (Pawlowsky-Glahn et al., 2015). The log- ratio methodology proposed by Aitchison (1986) repre- sents a powerful set of methods and techniques to apply to CoDa. During last decades, numerous innovative ideas and strategies to CoDa were presented at the four CoDaWork meetings (e.g., Martín-Fernández & Thió- Henestrosa, 2016a) and collected in special publications (e.g., Martín-Fernández & Thió-Henestrosa, 2016b). The approach adopted in this paper follows the principle of working on coordinates (Mateu-Figueras et al. 2011), that is, the standard statistical analysis is conveniently performed after choosing orthonormal log-ratio coordi- nates. In particular, we expressed each D-vector x=(x1, . . . , xD) of percentages of species as a (D-1)-dimensional real vector y=(y1, . . . , yD-1) of isomet- ric log-ratio coordinates (y =ilr(x)) (Egozcue et al., 2003) (see Appendix 1 for definitions and details). To develop the CoDa-MAT, the following points were considered: 1) pre-processing techniques; 2) choice of the distance measure; 3) number of analogues. 2.1. Data pre-processing: subcomposition, amalgamation and zero replace- ment Modern assemblages datasets often include a large number of zeros. As an example, planktonic fo- raminifera coretop datasets include zeros which are partly related to the broadly latitudinal distribution of most species, partly to the fact that there are several rare species whose abundance is often below the detec- tion limit. The zero values present in the data require a pre-processing because the log-ratio methodology needs the data to be strictly positive (Aitchison, 1986). To reduce the number of zero values to be replaced, rarer species can be excluded from the assemblages, by considering a subcomposition of the original assem- blages (Aitchison, 1986). Moreover, it is possible to amalgamate species characterised by similar ecological requirements (Aitchison, 1986). The choice of a limited numbers of species is motivated, as pointed out by Buc- cianti & Esposito (2004), by the necessity to reduce the zero substitutions needed by the log ratio transforma- tions, and to avoid, as much as possible, the generation of extreme clusters of points, corresponding to small values (Tauber, 1999). As pointed out by Kucera et al. (2005) the main problem with abundances of rare spe- cies is the inevitably low signal to noise ratio. Thus, amalgamation should not involve a reduction in the ac- curacy of results. About this, it is also worth recalling the subcompositional dominance property of an appropriate distance (Palarea-Albaladejo et al., 2012): distances between observations obtained from subcompositions, are less to equal to those obtained from full composi- tions. To manage the zeros occurring in the data, the first step of the analysis is the conversion of species vector of counts c, both for the fossil and modern data, into relative abundance compositional vectors x. If it can be assumed that zero values correspond to count zero, indicating that the component is absent by a sample size effect, it is possible to replace them by a suitable small value (Martín-Fernández et al., 2015). To do that, we adopted a mixed Bayesian-multiplicative estimation ap- proach, which is recommended when the compositional data arise from counts (Martín-Fernández et al., 2015). This treatment consists on a Bayesian estimation of the zero percentages combined with a multiplicative read- justment (Martín-Fernández et al., 2003) of the non-zero values (see Appendix 1). After the zero replacement, the ilr-coordinates’ vector y for the fossil and modern data are obtained. Thus, fossil assemblages are represented by its ilr- coordinates, whereas the modern database is consisting of ilr-coordinates together with the environmental pa- rameters measured at each location (at sea surface for planktonic assemblages coretops). 2.2 Distance measure The ilr-coordinates are isometric, i.e., an isometry is established between x and its real vector y. Hence, distances in the sample space of CoDa (the simplex S D) are associated with distances in RD-1. This property has important consequences for our aims, since distances in 61 Palaeoenvironmental reconstructions through CoDA where n is the number of modern samples, P represents the vector of k measured environmental parameters, and P̂ the corresponding vector of estimated values. This approach is equivalent to analyse the mean of the norm of the residual rows (difference between meas- ured and estimated parameters). The best result of MSD is obtained when MSD=0, i.e, all the estimates are equal to the measured values. In consequence, we se- lect as optimal number of analogues the number that produces the minimal MSD value. MSD index is the multivariate generalization of the Mean of Squared Er- rors (MSE), applied when only one parameter is evalu- ated. The square root of MSE is the Root Mean Square Error of Prediction (RMSEP), a quality index commonly calculated as a measure of the predictive abilities of the training set (Wallach & Goffmet, 1989; Birks, 1995). The multivariate approach should be applied when there is a set of parameters defining the palaeoclimate, as is usu- ally done in pollen analysis, from which several variables are estimated (i.e. temperature, seasonal or annual pre- cipitation, potential evapo-traspiration). In such a case, the environmental parameters may be normalized be- fore the computation of MSD to avoid dominance effects induced by the different measurement units adopted. Using the standard Pearson correlation index, we can calculate a measure for k environmental parame- ters. This procedure calculates the Pearson correlation coefficient r between a measured environmental pa- rameter p and its estimates p̂ for all the modern ana- logues. A perfect quality of the estimations will produce a Pearson correlation equal to 1. In a multivariate ap- proach, when one deal with a vector P of measured en- vironmental parameters (k parameters) and the corre- sponding vector of estimation P̂ , an index R is defined by ˆ( ) ( ) 1 (3) k P j P j j R r k = = where P(j) and P̂(j) respectively represents the jth pa- rameter. In this case, the best value of R corresponds to 1. In consequence, the number of analogues that pro- vides the maximum value of R is selected as optimal. It can be noted that the Pearson correlation coefficient r is used as a similarity index between a measured and an estimated parameter. However, the results with this index could be misleading in some scenarios. For exam- ple, when the estimates are closely proportional to the measured values this index will also show a reasonable value (close to one), although the estimations are low- quality estimations because they are systematically bi- ased. By combining R index with MSD index, one com- pletes the information about the quality and avoid these undesirable situations because MSD index measures the goodness of the estimates across the samples 2.4. Application on fossil assemblages The above-described leave-one-out evaluation provides the basis for the application of CoDa-MAT on fossil assemblages, as best performing number of ana- logues for the modern conditions was determined by means of MSD and R indexes. In a second step, it can be evaluated if this optimal number of analogues can be adopted to perform palaeoestimates from fossil assem- blages. As in the RAM method (Waelbroeck et al., 1998), the number of analogues may be evaluated for each assemblage, by looking for “jumps”, i.e. increases in distance larger than a fraction of the last modern ana- logue selected. The number of analogues may be also reduced to obtain “oceanographically coherent” ana- the simplex are translated into ordinary ones in the space of coordinates. Thus, the Aitchison distance (da) between two compositional vectors x, x* ∈ S D (Aitchison, 2000) is equal to the Euclidean distance (de) between their ilr-coordinates vectors y, y* ∈ R D-1, i.e., da(x, x*)=de(y, y*) (Palarea-Albaladejo et al., 2012). The Aitchison distance meets the requirements of scale in- variance, perturbation invariance and subcompositional dominance that are needed to achieve a meaningful statistical analysis of CoDa (Palarea-Albaladejo et al., 2012). An alternative distance measure that also meets these requirements is the Mahalanobis distance on ilr- coordinates (Palarea-Albaladejo et al., 2012). On the other hand, the distances related to angular measures do not have a compositional coherent behaviour (Palarea-Albaladejo et al., 2012). This narrows it down to just Euclidean or Mahalanobis distance on ilr- coordinates the metrics that can be adopted in the CoDa-MAT. 2.3. How many analogues? In theory, a palaeoestimate for a fossil sample may be obtained from a single very close modern analogue. However, the application on a too restricted number of modern analogues increases the risk of a bad estimate due to presence of anomalous samples in the modern dataset (i.e. arising from a decoupling of assemblage from the underlying surface conditions, or to taphonomic processes). In consequence, in the application of MAT is essential to determine the optimal number of ana- logues from which the palaeoestimates are obtained. In the CoDa-MAT this subject is handled by first one evalu- ating which number of analogues performs better in the modern conditions, by means of sensitivity analysis of leave-one-out cross-validation method (ter Braak & Jug- gins, 1993; Barrows & Juggins, 2005). Once this num- ber is determined, it is also adopted for palaeoestimates of all fossil assemblages. The sensitivity analysis to define the optimal num- ber of analogues is based on numerical indices. Indeed, once the environmental estimates are obtained for all the modern samples these are compared with the corre- sponding measured values, by computing two indices: the Pearson correlation coefficient and the mean squared distances (MSD) (e.g., Martín-Fernández et al., 2003). The MSD is a multivariate index of quality of the estimates of k environmental parameters defined as: ( ) ( ) 22 1 ˆ( ), ( ) 2 n e i d i i MSD n ==  P P logues (in the case of marine assemblages such as planktonic foraminifera - see section 3.2) or to achieve coherence of the vegetational context of the analogues (in the case of continental proxies such as pollens). 2.4.1. No-analogue conditions, outliers and atypical- ity index. Probably the most important difficulty of any proxy- based reconstruction is represented by the no-analogue problem, occurring when the palaeoenvironmental con- ditions represented in the fossil assemblages do not have a correspondence in modern environments. More in general, it is important to define quality indexes of the palaeoestimates. The most obvious ones are repre- sented by the average distance from the fossil assem- blages to its nearest modern analogues and by the stan- dard deviation of associated to the estimates. Large mean distances and high standard deviation (often as- sociated with oceanographical incoherence) likely indi- cate bad estimates. Overpeck et al. (1985) indicated threshold values to be adopted with squared chord dis- tance. In the CoDa-MAT the problem of detecting no- analogue conditions is faced by first considering the atypicality index and outlier detection for CoDa sets. In a broad sense, atypicality index (e.g., Aitchison, 1986) represents the probability of a more typical composition (or with a smaller Mahalanobis distance, from the centre of the dataset) than the observed one. In CoDa-MAT, atypicality index is first computed for each modern sam- ple using Mahalanobis distances, to determine critical values of distance for which a sample can be consid- ered a potential outlier. After that, atypicality can be computed for each fossil assemblage to test the hy- pothesis of a significant difference respect to the popu- lation of modern ones. A significant difference may high- light no-analogue conditions. Note that, in CoDa-MAT, this index is based on the Mahalanobis distance applied to ilr real vectors. This strategy is also in the basis of the robust methods to outlier detection for CoDa sets (Filzmoser & Hron, 2008) that are freely available in the R packages mvoutlier and robCompositions (R Develop- ment Core Team, 2011). Robust methods to identify multivariate outliers are based on the robust estimation of the covariance structure (e.g., Peña & Prieto, 2001). In our work, the estimated ilr covariance structure was used to assign a robust Mahalanobis distance to each fossil assemblage indicating how far the sample is from the centre of the modern data cloud. A second approach to outlier detection has been developed by computing the Local Outlier Factor (LOF) (Breunig et al., 2000) of assemblages. This index is based on a k-distance neighbourhood approach and assigns to each object the degree of being an outlier by determining how isolated the object is with respect to the surrounding neighbourhood. Values much larger than 1 are associated to suspect local outliers. As done for the atypicality, the LOF has been determined on ilr real vectors. The computations were replicated with increasing k values, which yielded quite stable results. The following results are based on a k=10. The relation- ship between a fossil sample and modern population can be approximately displayed in few dimensions by means of relative variation biplots (RVB) in which it is also possible to include fossil data as supplementary elements (Daunis-i-Estadella et al., 2011). In this way it is possible to easily visualize the space distribution of fossil samples and their modern counterparts. This may be useful to verify if a fossil sample is located within or at the margin of the surrounding closest modern ana- logues (see section 3.2). Matlab routines to obtain RVB and perform the CoDa-MAT are provided as supplemen- tary materials. 2.4.2. Evaluation of errors of estimates in relation to hyperspheres radius This approach provides for each sample an auxil- iary tool for evaluating the reliability of the obtained esti- mates. For each modern sample, the maximum value of the module of the errors of the estimates is computed within hyperspheres with increasing radii centred on the sample itself. The obtained information is applied to fossil assemblages by considering them as laying at the border of h hyperspheres, each centred on one of the h closest analogues. The highest value of the maximum errors associated to these hyperspheres is recorded as potential estimate error, as a measure of the reliability of the obtained estimates for the fossil assemblages. Con- fident estimates are obtained when a fossil assemblage is included within the radius of modern assemblages for which errors and standard deviation of the estimates are low. On the contrary, high values suggest considering with caution the reconstructed values or, at the limit, to discard these estimates if the potential error exceeds acceptable values. 62 Fig. 1 - Location of modern planktonic foraminifera coretop samples adopted for application of CoDa-MAT. Drawn with Ocean Data View software (Schlitzer, 2018). Di Donato V. et al. 3. TESTING THE METHOD The performance of CoDa-MAT in the modern conditions was tested on planktonic foraminifera with a modern dataset consisting of 1252 Atlantic and Mediter- ranean coretop assemblages determined on the >150 μm size fraction (Prell et al., 1999; Hayes et al., 2004; Kucera et al., 2004) (Figure 1). Additional Mediterra- nean Sea coretop assemblages were added from avail- able literature datasets and unpublished data. The oceanographical data, consisting of mean annual, ca- loric summer and caloric winter SST refer to Antonov et al. (2010) and Locarnini et al. (2010). The SST values at coretops location were computed by means of Ocean Data View 4.7.10 (Schlitzer, 2018). Following Kucera et al. (2005), oceanographical data are related to a depth of 10 m. About the modern dataset, it should be noted that the adopted zero substitution approach can't be directly applied to literature data when the planktonic foraminifera assemblages are reported as percentages without the total of counted specimens. Although not an ideal solution (the obvious alternative being to exclude these samples), in few cases the percentages data in- cluded in the modern dataset were “reconverted” into counts by assigning a total of 300 specimens to the 63 Fig. 2 - Relative variation biplot of planktonic foraminiferal as- semblages included in the modern dataset (19 taxonomical groups). Samples are grouped according to their latitude. Fig. 3 - Sensitivity analysis of leave-one-out verification in rela- tion to the number of analogs adopted for the CoDa-MAT esti- mate of the modern measured seasonal SST. a) MSD index; b) R index; c) mean std of estimates. Black: annual SST; red: sum- mer SST; blue: winter SST. Solid line: ilr-Euclidean (Aitchison) distance; Dashed lines: ilr-Mahalanobis distance. Tab. 1 Leave-one-out estimated MSD index, R index and mean std of estimates (6 analogs) Palaeoenvironmental reconstructions through CoDA samples. The adopted modern planktonic foraminiferal assemblage dataset consists of 26 species. With the approach described in Section 2.1 two datasets were generated, consisting of respectively 19 and 15 taxo- nomical groups (see Appendix 2). The RVB of the planktonic foraminifera assemblages included in the modern dataset is shown in Figure 2. The first two axes accounts for about 63% of total variability (71% if a third axis is added). The pattern of the links connecting col- umn points highlights the close relationship (i.e. short links) among low latitude taxa such as Globigerinoides ruber, Globogerinoides sacculifer and Globorotalia menardii-tumida, all of which are located on the positive side of axis 1. On the contrary the polar Neoblobo- quadrina pachyderma left coiled is characterised by a strong log-ratio variability respect to most of the taxa, except for Turborotalita quinqueloba. The location of the samples points in the RVB highlights the well-known broad latitudinal distribution of assemblages. 3.1. Leave-one-out verification The leave-one-out verification (see par. 2.4.1) was carried out on the 15-parts dataset by computing a sen- sitivity analysis of MSD and R indexes in relation to the number of analogues (h). Both Aitchison and log-ratio Mahalanobis distances were considered. In considera- tion of the high correlation between annual and sea- sonal SST, the MSD and R indexes were only tested with a univariate approach. For this sensitivity analysis a number of analogues h in (3) varying from 1 to 30 was considered to obtain the best value of the indexes. For each value of h (from 1 to 30) the estimates of the envi- ronmental parameter of each modern data sample were calculated using a leave-one-out procedure. As shown in Figure 3, lower MSD and higher R values were ob- tained with 6 to 10 analogues, both for unweighted esti- mates and estimates weighted on the inverse of the Aitchison distance. As expected, the weighted estimates were characterised by lower MSD values and higher R indexes, with peak values for annual SST estimates of 0.993. The annual SST R indexes are higher than sea- sonal ones. Moreover, it can be noted, for both R and MSD indexes, a rather flat response to the increase in the number of analogues (Figure 3). Although compara- ble results are obtained from 6 to 10 analogues, the progressive increase in the standard deviation of esti- mates detected for increasing number of analogues, suggests adopting 6 analogues for the estimates. Figure 4 and Table 1 summarise the relationships between measured seasonal SST and its estimates obtained under these conditions. The maximum mean Aitchison distance between a sample and its closest neighbour- hood is about 4.7. The leave-one-out verification was also tested with the Mahalanobis distance on ilr data. However, in com- parison with Aitchison distance, the MSD and standard deviation values are always higher, and the R indexes lower (Figure 2). Overall, the sensitivity analysis indi- cates that the Aitchison distance is more performing than the log-ratio Mahalanobis distance. We replicated the sensitivity analysis by testing the 19 parts dataset. The results are quite similar (the correlation coefficient between the annual SST obtained with 19 and 15 parts is 0.999) and are not reported here. 3.2. Spatial and geographical relationships. As discussed above, on the whole dataset, the best leave-one out estimates are obtained with 6 to 10 analogues. Further, the leave-one-out verification was extended to include the whole analogues, by also con- sidering the geographical location of the analogues. In general, the increase in the number of analogues corre- sponds to an increase in the standard deviation of the estimates. As expected, it can be noted a correspon- dence between increasing Aitchison and geographical distances (Figure 5). For middle to latitude samples, the sample being tested tends to remain in between the analogues (in terms of both spatial and geographical distribution). With increasing Aitchison distances, far- thest away analogues appear, mostly symmetrical with respect to the Equator. However, the behaviour is differ- ent for high latitude samples, as the sample being tested 64 Fig. 4 - Plots of observed versus estimated SST (annual and seasonal) obtained with 6 analogues. The error bars represent the std of estimates. Di Donato V. et al. 65 Fig. 5 - The figure shows, for 4 modern samples located at different latitudes, a) for progressively far analogues, the coordinates on the first 2 axes of the RVB shown in figure 2 and b) their geographical location. The graphs on the left (c) show the differences in SST between a sample and its analogues (blue dots), the SST error for estimates obtained from an increasing number of analogues (green line), and from moving averages of 6 progressively far analogues (red line). Palaeoenvironmental reconstructions through CoDA get into an off-center position as progressively far sam- ples are added. The behaviour is also different for what concern the error of the estimates. For samples located at the border of the cloud of modern data, such as high latitude samples, the increase in the number of ana- logues corresponds to a progressive increase in the error of the estimates (Figure 5). Instead, for middle and low latitude samples the increase is less marked. More- over, in several cases, the increase in the number of the analogues (and, consequently, in the Aitchison dis- tances) does not imply further decreasing in the quality of estimates. Likely, this behaviour arises from the fact that the SST differences obtained moving away from the hypersphere’s centre in opposite directions tend to com- pensate each other. These results suggest evaluating, in the application of the method to fossil assemblages, the oceanographical coherency of modern analogues and the spatial relationships between a fossil samples and its modern counterparts. More confident estimates are likely obtained when a fossil samples is in between its modern analogues, and when the latter are found into a limited geographical region. Summing up, threshold distances for evaluating suspect no-analogue conditions can be determined from atypicality index. In addition, confident estimates should fulfil the criteria of: a) closeness and oceanographical/ geographical coherence of the analogues, b) reduced standard deviation of the estimates, c) staying in the middle of a fossil samples respect to its modern ana- logues. These results agree with the findings of Gujot & de Vernal (2011a; 2011b), as geographical closeness among modern analogues appear a prerequisite for obtaining confident estimates. 3.3. Application examples As an application example, the CoDa-MAT method was applied to two different planktonic foraminifera re- cords. The first one is a literature dataset, consisting of the very detailed record of planktonic foraminifera as- semblages of the core MD95-2040 (de Abreu et al., 2003; Voelker & de Abreu, 2011), recovered in the At- lantic Ocean off the Iberian margin, the second is that of GNS84-C106 core recovered in the Tyrrhenian sea (Buccheri et al, 2002; Di Donato et al., 2008; 2009). Both datasets are obtained from >150 μm size fractions. As noted by Di Donato et al. (2015), the drawback rep- resented by the excessive loss of small sized species in this size fraction can be circumvented, at least in part, by means of CoDA methods. The fact remains never- theless that the unavailability of large modern datasets based on smaller size fractions imposes the choice of the >150 μm one which is not an optimal solution. The atypicality index computed with log-ratio Mahalanobis distances for each sample of the modern database (see section 2.4.1) allowed us to determine a critical value of 26.12 beyond which samples can be considered as outliers. By means of standard and robust atypicality, 91 and 303 outlier samples in the modern database were detected, respectively (Figure 6). These results appear rather restrictive because large part of modern samples is detected as outliers. Note that, by definition, the Ma- halanobis distance is a measure of global atypicality, Fig. 6 - DDplot of squared and robust squared Mahalanobis distances for a) modern data set b) Core MD95-2040 assem- blages c) Core GNS84-C106. The red lines represent the 97.5% percentile level, associated to xi-square value. 66 Di Donato V. et al. Fi g. 7 - R ec on st ru ct io n of s ea so na l S S T f or t he la st 2 10 k a of f th e Ib er ia n m ar gi n fro m c or e M D 95 -2 04 0 an d co m pa ris on b et w ee n C oD a- M A T an d S IM M A X 28 ( de A br eu e t al ., 20 03 , V oe lk er a nd d e A br eu , 2 01 1) r ec on st ru ct ed S S Ts . a : s um m er S S T. b ) w in te r S S T. c ) gr ay -s ha de d do ts : d is ta nc e of fo ss il as se m bl ag es fr om e ac h of th e 6 cl os es t m od er n an al og ue s. F ul l lin e: m ea n va lu es . d ) LO F va lu es e ) at yp ic al ity in de x: 0 : n ot s ig ni fic an t; 1: s ig ni fic an t f -g ) gr ay -s ha de d do ts : m ax im um m od ul es o f t he e st im at e er ro rs fo r th e cl os es t 6 m od er n an al og ue s, at r ad ii eq ua l t o th e m ea su re d A itc hi so n di st an ce . F ul l l in e: M ea n va lu es . h ) A lk en on e ba se d S S T re co ns tru ct io n (P ai lle r & B ar d, 2 00 2) . j ) G lo bi ge rin a bu llo id es s ta bl e is ot op e re co rd a nd M ar in e Is ot op ic S ta ge s (M IS ) (A br eu e t a l., 2 00 3; S ch ön fe ld e t a l., 2 00 3) . 67 Palaeoenvironmental reconstructions through CoDA that is, it is affected by the existence of groups in the data. In contrast, using the LOF of modern assem- blages, only two samples have values larger than 2, and only 35 samples values higher than 1.4. 3.2.1. Atlantic Ocean The foraminiferal record of MD95-2040 core covers the last 210 ka (de Abreu et al., 2003; Voelker, & de Abreu, 2011). Sea surface temperatures for this interval were formerly reconstructed (de Abreu et al., 2003) from planktonic foraminifera with SIMMAX28 method (Pflaumann et al., 1996). As regards the atypicality of assemblages, for about 19% of the samples, the Maha- lanobis distances are above the xi-square critical value of 26.11 corresponding to the 97.5 percentile (Figure 6). In relation to the 99.5 percentile, 4% of the samples have Mahalanobis distances are above the xi-square critical value of 31.32. The robust outlier detection is more restrictive, as about 40% of samples lies beyond the critical value. As for the LOF, glacial assemblages are characterised by higher values of up to 2, while most interglacial assemblages have LOF values not exceed- ing 1.5. Coherently, the mean distance of closest ana- logues is lower for interglacial periods. However, the expected maximum errors of the estimates are rather stable throughout the core (Figure 7). The comparison with the values reconstructed for summer and winter SST by means of SIMMAX28 (Figure 7) (de Abreu et al., 2003), highlights a same general trend, with decreasing SST around the transition between MIS 7 and MIS 6. From MIS 6 to MIS 4, both methods record colder inter- vals (around 160 ka, 130 ka, 85 ka and 60 ka BP) which are also evident in the alkenone-based SST reconstruc- tions and in the oxygen stable isotopes record. In addi- tion, CoDa-MAT gives evidence to a negative SST peak around 108 ka BP, which is also evident in the al- kenones and δ18O record. SIMMAX28 and CoDa-MAT reconstructions are quite different for the interval between 50 ka and 20 ka BP. Although the reconstructed values for the most prominent SST decreases (at 40 ka, 30 ka and 16 ka BP) are similar, for this time interval CoDa-MAT pro- vides less fluctuating SST values, which appear more Fig. 8 - Geographical (left column) and spatial (right column) relationships among selected foraminiferal assemblages of Core MD95-2040 and the modern coretop assemblages. Numbers indicate the ranking of the closest modern analogues. Fig. 9 - Location of regionalised modern planktonic foraminifera coretop assemblages adopted for application of CoDa-MAT to the Core GNS84-C106. Drawn with Ocean Data View software (Schlitzer 2018). 68 Di Donato V. et al. Fi g. 1 0 - C oD a- M A T re co ns tru ct io n of s ea so na l S S T ob ta in ed fr om G N S 84 -C 10 6 C or e. a ) su m m er b ) w in te r. T he e rr or b ar s in di ca te th e st an da rd d ev ia tio n of e ac h re co ns tru ct ed v al ue . c ) gr ay -s ha de d do ts : d is ta nc e of fo ss il as se m bl ag es fr om e ac h of th e 6 cl os es t m od er n an al og ue s. F ul l l in e: m ea n va lu es . d ) LO F va lu es e ) at yp ic al ity in de x: 0 : n ot s ig ni fic an t; 1: s ig ni fic an t f - g) g ra y- sh ad ed d ot s: m ax im um m od ul es o f th e es tim at e er ro rs f or t he c lo se st 6 m od er n an al og ue s, a t ra di i e qu al t o th e m ea su re d A itc hi so n di st an ce . Fu ll lin e: M ea n va lu es . Th e IN TI - M A TE G re en la nd e ve nt s tra tig ra ph y is re po rte d fro m R as m us se n et a l. (2 01 4) . 69 Palaeoenvironmental reconstructions through CoDA coherent with the alkenone-based SST reconstructions and δ18O stable isotope record. The spatial and geo- graphical relationships between selected fossil assem- blages of the Core MD95-2040 and their closest modern analogues are shown in Figure 8. It can be noted that the modern analogues of some fossil assemblages char- acterised by higher LOF, are located in the South and North Atlantic Ocean, mostly symmetrical with respect to the Equator. 3.2.2. Tyrrhenian sea (Western Mediterranean Sea) The Core GNS84-C106 core covers the last 34 ka (Di Donato et al., 2009). Following Kucera et al. (2005), the application of the CoDa-MAT was carried out on a regionalised modern dataset represented by Mediterra- nean and North Atlantic coretop assemblages (Figure 9). The amplitude of SST changes (Figure 10) is of about 10°C, the lower values being recorded between 15 and 17.5 cal ka BP. It can be noted that Holocene samples younger than 6 cal ka BP are characterised by low standard deviation, low mean distances from the closest modern analogues and low LOF. Early Holocene assemblages are characterised by slightly higher mean distances, albeit the low LOF values. On the contrary, the increasing mean distances obtained for late glacial and last glacial period (LGP) assemblages corresponds to an increase in the standard deviation of the estimates. Moreover, the mean distances reach values of up to 5, which are beyond the maximum value obtained for mod- ern assemblages. As regards the atypicality of fossil assemblages, with a critical value of 21.92, several as- semblages of this core, in particular from the LGP, can be considered as outliers with respect to the modern dataset (Figure 6). Comparable results are obtained by considering the expected errors of the estimates ob- tained by means hypersphere radius evaluation (Figure 10). The “expected errors” are below 1°C for the last 5 ka. The higher mean distances recorded in the early Holocene (between 9 and 5 ka BP) corresponds to a first increase in expected errors, with values ranging from 1° to 2°C. Noteworthy this interval is coeval with the eastern Mediterranean sapropel S1 event, during which the Mediterranean oceanographical asset was quite different from today (e.g. Emeis et al., 2000; Ni Fhlaithearta et al., 2010). Higher values, in the order of 2°C were found for the LGP. The spatial and geographi- cal relationships between selected fossil assemblages and their modern counterparts are shown in Figure 11. It can be noted that when the modern analogues are found at higher distances, there are some oceano- graphic incoherencies. As an example, it can be noted that the closest modern analogues of an LGP assem- blage dated at 21.3 cal ka are partly found in the Atlantic Ocean, partly in the Mediterranean Sea. Moreover, in addition to being characterised by high LOF and signifi- cant atypicality, in the RVB plot this fossil assemblages seems lying off centre with respect to its modern ana- logues. A marked oceanographic incoherence also char- acterises the assemblage dated at 13.5 cal ka. In sum- mary, atypicality of samples and oceanographic inco- herency, suggests caution in the interpretation of the reconstructed SST, for the Late Glacial and the LGP intervals of this core, at least for the >150 μm size frac- tion adopted for the analysis. The results obtained with CoDA-MAT have been compared with those obtained with raw data analysis and squared chord distance adopted as similarity measure (Figure 10). The results appear comparable for the Holocene interval of the core, Fig. 11 - Geographical (left column) and spatial (right column) relationships among selected foraminiferal assemblages of Core GNS84-C106 and the modern coretop assemblages. Numbers indicate the ranking of the closest modern analogues. 70 Di Donato V. et al. however for the Late Glacial and Last Glacial Period, the CoDa-MAT reconstructed SST are much lower than those obtained with raw data MAT, which appear even higher than present for the 19-25 cal ka interval. For this case study, the results are likely influenced by the size fraction commonly adopted for the modern planktonic foraminiferal coretop databases, i.e the excessive loss of small sized specimens of T. quinqueloba in the >150- micron size fraction (Di Donato et al., 2015). In this re- spect, it should be noted that are currently not available core top databases of planktonic foraminifera assem- blages built with smaller size fractions. However, it is likely that the adoption of smaller size fractions would produce lower SST estimates for the LGP intervals of the GNS84-C106 core. 4. CONCLUSIVE REMARKS In general, the leave-one-out cross validation per- formed on CoDa-MAT with modern data yielded correla- tions between estimated and measured modern values in the range of methods such as Imbrie and Kipp trans- fer functions, SIMMAX or RAM if applied to Atlantic Ocean (Waelbroeck et al., 1998). However, as pointed out above, the proposal of a method based on CoDa has not the objective of increase the correlation ob- tained from previous methods, but rather to develop the method coherently with the nature of data. Up to now, despite the advances in the statistical theory, only a relatively few micropaleontological studies adopted an approach coherent with CoDa principles (Buccianti & Esposito, 2004; Di Donato et al., 2009; Sgarrella et al., 2012; Di Donato et al., 2018). Apart from these funda- mental issues, CoDa-MAT also provides tools for evalu- ate the reliability of estimates performed on fossil as- semblages and the occurrence of atypical fossil assem- blages. The application of CoDa-MAT to the examples highlighted a coherent behaviour and, for the Iberian margin, a strong coherence with the stable isotopes record. About the Tyrrhenian sea case study, the cau- tion is suggested regarding LGP SST estimates. How- ever, the CoDA approach seems less subjected to prob- lems related to preparation techniques (i.e. size fraction) than raw data analysis. Although application examples of CoDa-MAT are described in a palaeoceanographic context, this method can be obviously applied to differ- ent proxies, i.e. to reconstruct past atmospheric environ- ments, as usually done in palynological studies. Further studies, involving the application of CoDa-MAT to conti- nental proxies will be carried out to evaluate its perform- ing. ACKNOWLEDGEMENTS We would like to thank Prof. Caterina Morigi and an anonymous reviewer for their kindly and constructive comments, which helped us to improve the manuscript. This work has been supported by the project “CODA- RETOS” (Spanish Ministry of Economy and Competi- tiveness; Ref: MTM2015-65016-C2-1-R). REFERENCES Aitchison J. (1986) - The Statistical Analysis of Compo- sitional Data. Monographs on Statistics and Ap- plied Probability. Chapman & Hall Ltd. (Reprinted 2003 with additional material by The Blackburn Press), London (UK), pp. 416. Aitchison J., Barceló-Vidal C., Martín-Fernández J.A., Pawlowsky-Glahn V. (2000) - Logratio analysis and compositional distance. Mathematical Geol- ogy, 32 (3), 271-275. Antonov J.I., Seidov D., Boyer T.P., Locarnini R.A., Mishonov A.V., Garcia H.E. (2010) - World Ocean Atlas 2009 Volume 2: Salinity. S. Levitus, Ed., NOAA Atlas NESDIS 69, U.S. Government Print- ing Office, Washington, D.C., pp. 184. Barceló-Vidal C, Martín-Fernández J.A. (2016) - The mathematics of compositional analysis. Austrian journal of statistics, 45(4), 57-71. Barrows T.T, Juggins S. (2005) - Sea-surface tempera- tures around the Australian margin and Indian Ocean during the last glacial maximum. Quater- nary Science Reviews, 24, 1017-1047. Birks H.J.B. (1995) - Quantitative palaeoenvironmental reconstructions. In: Maddy D., Brew J.S. (Eds.), Statistical Modelling of Quaternary Science Data. Technical Guide 5. Quaternary Research Associa- tion, Cambridge, 116-254. Brassell S.C., Eglinton G., Marlowe I.T., Pflaumann U., Sarnthein M. (1986) - Molecular stratigraphy: A new tool for climatic assessment. Nature, 320 (6058), 129-133. Breunig M.M., Kriegel H.P., Ng R.T., Sander J. (2000) - LOF: identifying density-based local outliers. In Proceedings of the 2000 ACM SIGMOD interna- tional conference on Management of data (SIGMOD '00). ACM, New York, NY, USA, 93-104. Doi: http://dx.doi.org/10.1145/342009.335388. Buccheri G., Capretto G., Di Donato V., Esposito P., Ferruzza G., Pescatore T., Russo Ermolli E., Senatore, M.R., Sprovieri, M., Bertoldo, M., Carella D., Madonna G. (2002) - A high resolution record of the last deglaciation in the southern Tyrrhenian Sea: environmental and climatic evolution. Marine Geology, 186, 447-470. Buccianti A., Mateu-Figueras G., Pawlowsky-Glahn V. (Eds.) (2006) - Compositional Data Analysis in the Geosciences: From Theory to Practice. In: Special Publications, 264. Geological Society, London. Buccianti A., Esposito P. (2004) - Insights on Late Qua- ternary calcareous nannoplankton assemblages under the theory on statistical analysis for compo- sitional data: an application example. Palaeo- geography, Palaeoclimatology, Palaeoecology, 202, 209-227. Daunis-i-Estadella J., Thió-Henestrosa S., Mateu- Figueras G. (2011 - Including supplementary ele- ments in a compositional biplot. Computers & Geo- sciences, 37, Issue 5, May 2011, 696-701. de Abreu L., Shackleton N.J., Schonfeld,J., Hall M., Chapman M. (2003) - Millennial-scale oceanic climate variability off the Western Iberian margin 71 Palaeoenvironmental reconstructions through CoDA during the last two glacial periods. Marine Geol- ogy, 196 (1-2), 1-20. Di Donato V., Esposito P., Garilli V., Naimo D., Buccheri G., Caffau M., Ciampo G., Greco A., Stanzione D. (2009) - Surface-bottom relationships in the Gulf of Salerno (Tyrrhenian Sea) over the last 34 kyr: Compositional data analysis of palaeontological proxies and geochemical evidence. Geobios, 42, 561-579. Di Donato V., Esposito P., Russo Ermolli E., Scarano A., Cheddadi R. (2008) - Coupled atmospheric and marine palaeoclimatic reconstruction for the last 35 ka in the Sele Plain-Gulf of Salerno area (southern Italy). Quaternary International, 190, 146 -157. Di Donato V., Insinga D.D., Iorio M., Molisso F., Rumolo P., Cardines C. Passaro S. (2018) - The palaeocli- matic and palaeoceanographic history of the Gulf of Taranto (Mediterranean Sea) in the last 15 ky. Global and Planetary Change, 172, 278-297. Doi: 10.1016/j.gloplacha.2018.10.014. Di Donato V., Martin-Fernandez J.A., Daunis-i-Estadella J. Esposito P. (2015) - Size fraction effects on planktonic foraminifera assemblages: a composi- tional contribution to the golden sieve rush. Mathe- matical Geosciences, 47 (4), 455-470. Egozcue J.J., Pawlowsky-Glahn V., Mateu-Figueraz G., Barceló-Vidal C. (2003) - Isometric logratio trans- formations for compositional data analysis. Mathe- matical Geology 35 (3), 279-300. Elderfield H., Ganssen G. (2000) - Past temperature and δ18O of surface ocean waters inferred from foraminiferalMg/Ca ratios. Nature, 405, 442-445. Emeis K.C., Struck U., Schulz H.M., Rosenberg R., Bernasconi S., Erlenkeuser H., Sakamoto T., Martinez-Ruiz F. (2000) - Temperature and salinity variations of Mediterranean Sea surface waters over the last 16,000 years from records of plank- tonic stable oxygen isotopes and alkenone unsatu- ration ratios. Palaeogeography, Palaeoclimatol- ogy, Palaeoecology, 158 (3-4), 259-280. Doi: 10.1016/S0031-0182(00)00053-5. Everitt B.S., Landau S., Leese M., Stahl D. (2011) - Cluster analysis. John Wiley & Sons, Ltd, Chiches- ter, United Kingdom, 5th Edition, pp. 330. Filzmoser P., Hron, K. (2008) - Outlier detection for compositional data using robust methods. Math. Geosciences, 40, 233-248. Guiot J., de Vernal A. (2011a) - Is spatial autocorrelation introducing biases in the apparent accuracy of paleoclimatic reconstructions?. Quaternary Sci- ence Reviews 30, 1965-1972. Guiot J., de Vernal A. (2011b) - QSR Correspondence “Is spatial autocorrelation introducing biases in the apparent accuracy of palaeoclimatic reconstruc- tions?” Reply to Telford and Birks. Quaternary Science Reviews, 30, 3214-3216. Hayes A., Kucera M., Kallel N., Sbaffi L., Rohling E. (2004) - Compilation of planktic foraminifera modern data from the Mediterranean Sea. Pan- gaea. Doi:10.1594/PANGAEA.227305. Hutson W.H. (1979) - The Agulhas Current during the Late Pleistocene: Analysis of modern faunal ana- logues. Science, 207 (1), 64-66. Imbrie J. Kipp N.G. (1971) - A new micropaleontological method for paleoclimatology: Application to a Late Pleistocene Caribbean core. The Late Cenozoic Glacial Ages. New Haven, Yale University Press. 71-181. Kucera M., Weinelt M., Kiefer T., Pflaumann U., Hayes A., Weinelt M., Chen M.T., Mix A.C., Barrows T., Cortijo E., Duprat J., Juggins S., Waelbroeck C. (2004) - Compilation of planktic foraminifera cen- sus data, modern from the Atlantic Ocean. Pan- gaea. Doi:10.1594/PANGAEA.227322. Kucera M., Rosell-Melé A., Schneider R., Waelbroeck C., Weinelt M. (2005) - Multiproxy approach for the reconstruction of the glacial ocean surface (MARGO). Quaternary Science Reviews, 24 (7-9), 813-819. Locarnini R.A., Mishonov A.V., Antonov J.I., Boyer T.P., Garcia H.E. (2010) - World Ocean Atlas 2009, Volume 1: Temperature. S. Levitus, Ed., NOAA Atlas NESDIS 68, U.S. Government Printing Of- fice, Washington, D.C., pp.184. Malmgren B.A., Kucera M., Nyberg J. Waelbroeck C. (2001) - Comparison of statistical and artificial neural network techniques for estimating past sea surface temperatures from planktonic foraminifer census data. Paleoceanography, 16 (5), 520-530. Doi:10.1029/2000PA000562. Martín-Fernández J.A., Barceló-Vidal C., Pawlowsky- Glahn V. (2003) - Dealing with zeros and missing values in compositional datasets using nonpara- metric imputation. Mathematical Geology, 35 (3), 253-278. Martín-Fernández J.A., Hron K., Templ M., Filzmoser P,. Palarea-Albaladejo J. (2015) - Bayesian- multiplicative treatment of count zeros in composi- tional datasets. Statistical Modelling, 15(2), 134- 158 Martín-Fernández J.A., Thió-Henestrosa S. (Eds.) (2016a) - Compositional Data Analysis: CoDa- Work, L’Escala, Spain, June 2015. Springer Pro- ceedings in Mathematics & Statistics, 187. Springer International Publishing, New York (USA). Martín-Fernández J.A., Thió-Henestrosa S. (Guest eds.) (2016b) - Compositional Data Analysis. Austrian Journal of Statistics, 45(4). Mateu-Figueras G., Pawlowsky-Glahn V., Egozcue J.J. (2011) - (Pawlowsky-Glahn V, Buccianti A., Eds). The Principle of Working on Coordinates, in Com- positional Data Analysis: Theory and Applications. John Wiley & Sons, Ltd, Chichester, UK, 29-42 Ni Fhlaithearta S., Reichart G.J., Jorissen F.J., Fontanier C., Rohling E.J., Thomson J., De Lange, G.J. (2010) - Reconstructing the sea floor environ- ment during sapropel formation using trace metals and sediment composition. Paleoceanography, 25, PA4225. Doi:10.1029/2009PA001869, 2010. 72 Di Donato V. et al. Overpeck J.T., Webb T., Prentice I.C. (1985) - Quantita- tive interpretation of fossil pollen spectra: dissimi- larity coefficients and the method of modern ana- logues. Quaternary Research, 23, 87-108. Pailler D., Bard E. (2002) - High frequency palaeocean- ographic changes during the past 140 000 yr re- corded by the organic matter in sediments of the Iberian Margin. Palaeogeography, Palaeoclimatol- ogy, Palaeoecology, 181, 431-452. Palarea-Albaladejo J., Martín-Fernández J.A., Soto J.A. (2012) - Dealing with Distances and Transforma- tions for Fuzzy C-Means Clustering of Composi- tional Data. Journal of Classification 29 (2), 144- 169. Pawlowsky-Glahn V, Egozcue J.J., Tolosana-Delgado R. (2015) - Modeling and analysis of compositional data. John Wiley & Sons, Chichester, pp. 378. Peña D., Prieto, F. (2001) - Multivariate outlier detection and robust covariance matrix estimation. Tech- nometrics, 43 (3), 286-310. Pflaumann U., Duprat J., Pujol C., Labeyrie L.D. (1996) - SIMMAX: a modern analogue technique to de- duce Atlantic sea surface temperatures from planktonic foraminifer in deep-sea sediments. Paleoceanography, 11, 15-36. Prell W., Martin A., Cullen J., Trend M. (1999) - The Brown University Foraminiferal Data Base. IGBP PAGES/World Data Center-A for Paleoclimatology Data Contribution Series # 1999-027. NOAA/ NGDC Paleoclimatology Program, Boulder, CO, USA. Rasmussen S.O., Bigler M., Blockley S.P., Blunier T., Buchardt S.L., Clausen H.B., Cvijanovic I., Dahl-Jensen D., Johnsen S.J., Fischer H., Gkinis V., Guillevic M., Hoek W.Z., Lowe J.J., Pedro J.B., Popp T., Seierstad I.K., Peder Steffensen J., Svensson A.M., Vallelonga P., Vinther B.M., Walker M.J.C., Wheatley J.J., Winstrup M. (2014) - A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event strati- graphy. Quaternary Science Reviews, 106, 14-28. R Development Core Team (2011) - R: A Language and Environment for Statistical Computing. R Founda- tion for Statistical Computing, Vienna, Austria. http://www.r-project.org. Schlitzer R. (2018) - Ocean Data View. http://odv.awi.de. Schönfeld J., Zahn R., de Abreu L. (2003) - Stable iso- tope ratios and foraminiferal abundance of sedi- ment cores from the Western Iberian Margin. Doi:10.1594/PANGAEA.733303. supplement to: Schönfeld J., Zahn R., de Abreu L. (2003) - Surface to deep water response to rapid climate changes at the western Iberian Margin. Global and Planetary Change, 36(4), 237-264. Doi:10.1016/S0921-8181(02)00197-2. Sgarrella F., Di Donato V., Sprovieri R. (2012) - Benthic foraminiferal assemblage turnover during intensifi- cation of the Northern Hemisphere glaciation in the Piacenzian Punta Piccola section (Southern Italy). Palaeogeography, Palaeoclimatology, Palaeoecol- ogy, 333-334, 59-74. Doi: 10.1016/j.palaeo.2012.03.009. Tauber F. (1999) - Spurious clusters in granulometric data caused by logratio transformation. Mathemati- cal Geology, 31, 491-504. Telford R.J., Birks H.J.B. (2009) - Evaluation of transfer functions in spatially structured environments. Quaternary Science Reviews, 28, 1309-1316. Telford R.J., Birks H.J.B. (2011) - Is spatial autocorrela- tion introducing biases in the apparent accuracy of palaeoclimatic reconstructions? Quaternary Sci- ence Reviews, 30, 3210-3213. ter Braak C.l.F., Juggins S. (1993) - Weighted averaging partial least squares regression (WA-PLS): an improved method for reconstructing environmental variables from species assemblages. Hydrobiolo- gia. 2691270, 485-502. Voelker A.H.L., de Abreu L. (2011) - A Review of Abrupt Climate Change Events in the Northeastern Atlan- tic Ocean (Iberian Margin): Latitudinal, Longitudi- nal and Vertical Gradients. In: Rashid H., Polyak L., Mosley-Thompson E. (Eds), Abrupt Climate Change: Mechanisms, Patterns, and Impacts. Geophysical Monograph Series (AGU, Washington D.C.), 193, 15-37. Doi:10.1029/2010GM001021. Waelbroeck C., Labeyrie L., Duplessy J.C., Guiot J. Labracherie M. (1998) - Improving past sea sur- face temperature estimates based on planktonic fossil faunas. Paleoceanography, 13, 272-283. Wallach D., Goffmet B. (1989) - Mean squared error of prediction as a criterion for evaluating and compar- ing system models. Ecological Modelling, 44, 299- 306. Ms. received: July 18, 2018 Final text received: October 9, 2018 73 Palaeoenvironmental reconstructions through CoDA 74