Marie Chavent

Jerome Lacaille

Alex Mourer

Madalina Olteanu

Sparse k-means for mixed data via group-sparse clustering

The present manuscript tackles the issue of variable selection for clustering, in high dimensional data described both by numerical and categorical features. First, we build upon the sparse k-means algorithm with lasso penalty, and introduce the group-L1 penalty -already known in regression -in the unsupervised context. Second, we preprocess mixed data and transform categorical features into groups of dummy variables with appropriate scaling, on which one may then apply the group-sparse clustering procedure. The proposed method performs simultaneously clustering and feature selection, and provides meaningful partitions and meaningful features, numerical and categorical, for describing them.

Introduction

Whereas the issue of feature selection through regularization procedures received a great deal of attention in the supervised learning context and resulted in an abundant literature over the last twenty years, it is only much later and relatively recently that it effectively emerged in the unsupervised framework. The first approaches were model-based, these being naturally suited for including lasso (L 1 ) and related penalties, and one may cite [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] for a L 1 -penalized EM procedure (the mixture consists of Gaussian distributions with equal variances) or [START_REF] Ch | Model-based clustering of high-dimensional data: A review[END_REF] for a detailed review on model-based clustering for high-dimensional data. In the more general framework where no assumption is made on the underlying distribution, a sparse k-means algorithm with L 1 penalty was introduced in [START_REF] Witten | A framework for feature selection in clustering[END_REF], and later on extended to feature selection within each cluster and reinforced by consistency results, [START_REF] Arias-Castro | A simple approach to sparse clustering[END_REF] [5] [START_REF] Sun | Regularized k-means clustering of high-dimensional data and its asymptotic consistency[END_REF]. Let us also mention that a generalization of the sparse k-means algorithm to overlapping groups of variables was recently introduced in [START_REF] Huo | Integrative sparse k-means with overlapping group lasso in genomic applications for disease subtype discovery[END_REF] .

That being said, all methods cited above are essentially designed for numerical data, while real data is often made of numerical and categorical features. Some of the authors above touch upon the question of categorical features, by mentioning the possibility of making them numerical using a transformation through dummy variables. However, this processing step is not that immediate, since the Euclidean distance on zero-one vectors is not particularly suited for being mixed with Euclidean distances on numerical variables. Other authors ESANN 2020 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning. Online event, 2-4 October 2020, i6doc.com publ., ISBN 978-2-87587-074-2. Available from http://www.i6doc.com/en/.

implicitly suggest that the proposed algorithms may be written in terms of distances or dissimilarities between input data only, and hence it suffices to use an appropriate distance for categorical features. Nevertheless, the distance-based approaches may rapidly translate into an increased complexity if the size of the data becomes large.

The present manuscript aims at proposing an explicit method for variable selection in a mixed-data framework, using a penalized criterion. Our contribution is two folded: first, we introduce a group-sparse version of the k-means algorithm, which allows one to select among priorly defined groups of numerical variables; second, we use the preprocessing proposed in [START_REF] Chavent | Multivariate analysis of mixed data: The pcamixdata r package[END_REF] for mixed data and write each categorical feature as a group of properly scaled numerical features, on which we apply the group-sparse k-means. Eventually, starting from a data set described by both numerical and categorical features, we simultaneously perform k-means clustering and feature selection on both types of features.

Method and data preprocessing

We start by introducing the group-sparse k-means algorithm, followed by the particular writing of the problem in the case of mixed data.

For the moment, let us suppose that the data X consists of n input vectors valued in R p , and that X j ∈ R n represents the j-th feature (all features are scaled to zero mean and unit variance). Furthermore, the data may be partitioned into K clusters, with n k the number and C k the set of indices of the observations in the k-th cluster. Then, one may define X j,k = 1 n k i∈C k X i,j the average of feature j in cluster k, and X j = 1 n n i=1 X i,j the average of the jth feature over the entire data set.

Since the k-means algorithm aims at finding the partition C 1 , ..., C K maximizing the between-class variance, let us write this as:

B(X, C 1 , ..., C K ) = K k=1 n k n p j=1 (X j,k -X j ) 2 = p j=1 B(X j , C 1 , ..., C K ) , (1)
where B(X j , C 1 , ..., C K ) is the between-class variance associated to the j-th feature. In the following, in order to simplify the notations, we shall use b

j = B(X j , C 1 , ..., C K ), j = 1, ..., p.
If one suspects that features do not equally contribute to the clustering, some being more informative than others, she may address this by introducing a weighted version of the between-class variance, penalized by a regularization term:

B(X, C 1 , ..., C K , w) = p j=1 w j b j -λh(w) = w T b -λh(w) , (2)
where w T = (w 1 , ..., w p ), w ≥ 0, w 2 ≤ 1 is the vector of weights, and b T = (b 1 , ..., b p ) are the feature-wise between-class variances, for a given partition C 1 , ..., C K , and λ is a hyper parameter. In the following, K the number of clusters as well as λ are supposed to be priorly fixed. In the sparse clustering algorithm introduced in [START_REF] Witten | A framework for feature selection in clustering[END_REF], the regularization term is chosen by analogy with the lasso framework, hence h(w) = w 1 . This leads to a sparse representation of the variables, with the weights w j directly translating the contribution of the associated features to the clustering.

Group sparse k-means

We generalize the original sparse clustering by introducing a group regularization framework. Suppose that the p features are furthermore divided into L priorly known groups, such that For group data, let us define a specific L 1 -group penalty, which has been already used in the regression framework [START_REF] Yi | Model selection and estimation in regression with grouped variables[END_REF],

X = X 1 |...|X L , with X ∈ R n×p ,
h(w) = w 1,group = L =1 α w 2 , (3)
where (α ) is a vector of weights applied to the groups of variables. In the literature on group-sparse regression, two common choices appear to emerge, either α = 1, ∀ = 1, ..., L, or α = √ p ∀ = 1, ..., L. The latter is penalizing each group by its size, and we shall use it in the following.

With the previous notations, the optimization problem writes as follows:

max w,C1,...,C K w T b -λ L =1 √ p w 2 , w 2 ≤ 1 , w ≥ 0 (4)
For a fixed number of clusters K and for a fixed hyper parameter λ, the problem in Equation 4 is solved using an iterative algorithm. After having set an initial value for the weight vector w (generally one uses weights with equal norms for the L groups), one is alternating the two following steps until convergence:

1. Keeping w fixed, find C 1 , ..., C K maximizing w T b, which is actually equivalent to performing a usual k-means training on the scaled features Xj = √ w j X j , j = 1, ..., p.

2. Keeping C 1 , ..., C K fixed, find w ≥ 0 maximizing B(X, C 1 , ..., C K , w) and such that w 2 ≤ 1. The solution may be analytically obtained using Lagrange multipliers and the KKT conditions, and may be expressed in terms of group soft-thresholding operators:

w =      b b 2 if λ = 0, 0 if λ > 0 and SG (b, √ pλ) 2 = 0, SG (b, √ pλ) SG (b, √ pλ) 2 if λ > 0 and SG (b, √ pλ) 2 = 0, (5)
where

SG (b, √ pλ) T = S G (b 1 , √ p 1 λ) T , ..., S G (b L , √ p L λ) T ∈ R p , (6)
√ p T = ( √ p 1 , ..., √ p L ), and S G (b , √ p λ) = b b 2 b 2 - √ p λ + ∈ R p , = 1
, ..., L, is the group soft-thresholding operator.

The group soft-thresholding operator will take out from the model all groups which have a norm of the corresponding between-class variance b smaller than the fixed threshold λ normalized by the size of the group, and shrink by the same amount the norms of the remaining groups of variables.

The hyper parameter λ may be tuned using various criteria for assessing the quality of the model, such as the gap statistic as described in [START_REF] Witten | A framework for feature selection in clustering[END_REF], the ratio of explained variance, etc. Usually, one considers a fine grained grid valued between zero and some upper bound chosen as the maximum b 2 when training a simple k-means procedure. Choosing the optimal λ and the optimal number of clusters K is not trivial, but, because of the reduced number of pages, we do not discuss here the different criteria in the literature, nor how they fit in the group-sparse framework.

Sparse clustering for mixed data

Eventually, let us discuss the case where the data is summarized by both numerical and categorical features. If one has n input data X described by d 1 categorical features and d 2 numerical ones, such that each of the categorical features has p j possible values, j = 1, ..., d 1 , she may transform each categorical feature X j into p j dummy variables Xj = ( X1 j , ..., Xpj j ) ∈ {0, 1} n×pj and thus define a natural group structure on the transformed data

Y = [Y 1 |...|Y d1+d2 ],
where Y = X for = 1, ..., d 1 , Y = X for = d 1 + 1, ..., d 1 + d 2 , with respective group sizes p T = (p 1 , ..., p d1 , 1, ..., 1) ∈ R d1+d2 .

When training group-sparse k means on the group structure above, this amounts to performing variable selection in a mixed-data context. Before applying the algorithm described in the previous section, Y must be properly preprocessed. This prior step is described for example in [START_REF] Chavent | Multivariate analysis of mixed data: The pcamixdata r package[END_REF]: numerical variables are scaled to zero mean and unit variance, while the dummy variables are centered and normalized by 1/ n n l,j , where n l,j is the number of input data taking the j-th value of the lth feature, or equivalently the sum over Xj . The scaling applied to the dummy variables actually leads to using a χ 2 distance on the categorical features, while the numerical features, after scaling, are compared through the usual Euclidean distance.

A real-life example

For illustrating the algorithm introduced above, we use the Statlog Heart dataset available in the UCI machine learning data repository [START_REF] Frank | UCI machine learning repository, statlog (heart) data set[END_REF]. The data consists of 270 inputs and is described by six numerical features (age, resting blood pressure ESANN 2020 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning. Online event, 2-4 October 2020, i6doc.com publ., ISBN 978-2-87587-074-2. Available from http://www.i6doc.com/en/.

-rest.b.p, serum cholesterol in mg/dl -ser.c, maximum heart rate achievedmax.h.r.a, ST depression induced by exercice relative to rest -oldp, number of major vessels colored by fluoroscopy -num.m.v ) and seven categorical ones (sex, chest pain type -chest.p.t, fasting blood sugar > 120mg/dl -fast.b.s, resting electrocardiographic results -rest.e.r, exercice induced angina -exe.i.a, the slope of the peak exercice ST segment -slope, thalassemia -thal ). One control variable assesses the presence or the absence of a heart disease and gives an a priori partition in two clusters. We used this knowledge to select the number of clusters K = 2.

The hyper parameter λ was varied on a grid taking values between 0 and the maximum between-class feature variance, extracted from a k-means algorithm trained without weighting. The regularization paths are illustrated in Figure 1. At the optimum value of the hyper parameter λ (selected using a heuristic on the ratio of explained variance), six features only are kept, four numerical, and two categorical. The averages within cluster for the numerical features, and the frequencies within cluster for the categorical ones, are displayed in Table 1 and show that all selected features are significantly different in the clusters. The method thus provides both meaningful clusters, and meaningful features -both numerical and categorical -for describing the clusters. Let us remark here also that the ratio of explained variance, all features included, varies very little when using the partition at the optimal value of λ: only 2% is lost when using the six selected features instead on the initial thirteen. ESANN 2020 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning. Online event, 2-4 October 2020, i6doc.com publ., ISBN 978-2-87587-074-2. Available from http://www.i6doc.com/en/.

  p being the size of group , and p 1 + ... + p L = p. The between-class variance vector b and the weight vector w can be also decomposed as b T = (b 1 , ..., b L ) and w T = (w 1 , ..., w L ), where b = (b p0+...+p -1 +1 , ..., b p0+...+p ) T and w = (w p0+...+p -1 +1 , ..., w p0+...+p ) T , with the notation p 0 = 0.

Fig. 1 :

 1 Fig.1: Feature regularization paths for the heart data. The optimum selected λ is plotted in the vertical dotted black line. The y axis is w 2 , the norm of the optimal weights associated to each feature.

  ESANN 2020 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning. Online event, 2-4 October 2020, i6doc.com publ., ISBN 978-2-87587-074-2. Available from http://www.i6doc.com/en/.

Conclusion and perspectives

We provide a complete procedure for simultaneously clustering mixed data, and selecting the most relevant features for the clustering. This procedure is illustrated on a real dataset and proves to be efficient. Further work is currently done on finding appropriate criteria for selecting the hyper parameter λ and the number of clusters K. An R-package implementing this method is also currently being developed and will be soon made available for the community.


Alex Mourer
email: alex.mourer@safrangroup.comjerome.lacaille@safrangroup.com

Jérôme Lacaille

Madalina Olteanu
email: madalina.olteanu@univ-paris1.fr

Marie Chavent
email: marie.chavent@math.u-bordeaux.fr

Automatic Detection of Rare Observations During Production Tests Using Statistical Models

Engines are verified through production tests before delivering them to customers. During those tests, numerous measures are taken on different parts of the engine, considering multiple physical parameters. Unexpected measures can be observed. For this very reason, it is important to assess if these unusual observations are statistically significant.

However, anomaly detection is a difficult problem in unsupervised learning. The obvious reason is that, unlike supervised classification, there is no ground truth against which we could evaluate results. Therefore, we propose a methodology based on two independent statistical algorithms to double check the results. One approach is the Isolation Forest model which is specific to anomaly detection and able to handle a large number of variables. The goal of the algorithm is to find rare items, events or observations which raise suspicions by differing significantly from the majority of the data and, at the same time, it discriminates non-informative variables to improve estimation. One main issue of Isolation Forest is its lack of interpretability. Within this scope, we extend the Shapley values, interpretation indicators, to the unsupervised context to interpret the model outputs.

The second approach is the Self-Organizing Map (SOM) model which has nice properties for data mining by providing both clustering and visual representation. The performance of the method and its interpretability depend on the chosen subset of variables. In this respect, we first implement a sparseweighted K-means to reduce the input space, allowing the SOM to give an interpretable discretized representation.

INTRODUCTION

As an aircraft engines manufacturer, Safran verifies all individual engines before delivering to the customer during production tests. Those bench test operations generate lots of measures for different parts of the engines, resulting in multiple physical parameters acquisitions. As we may encounter unexpected measures, it is important to detect their causes and relevance. We build statistical methods to reach this goal.

Variations between performances of engines are common. Nevertheless, the production tests that verify essential engine functions before delivering it to an airline company are done in different bench test cells, under different ambient conditions, etc. A thermodynamic model is applied to compensate for context variations but there still exist some second level residuals we may have to compensate to enhance the quality of the measurements. They essentially depend on test bench components like slave cowls, but also sites and suppliers.

Therefore, one of the objectives is to take into account test bench components effects. Furthermore, there is no universally admitted way to evaluate unsupervised anomaly detection algorithms results. Hence, we proposed a new methodology based on two different algorithms.

• Large number of variables (>50) make statistical estimation challenging, especially w.r.t. the small number of engines (591). Therefore, a specific algorithm for anomaly detection, named Isolation Forest Liu et al. (2008), is proposed. Isolation-based methods measures the probability to be isolated and anomalies are those that have the highest probability. In randomly generated binary trees, where instances are recursively partitioned, the trees produce noticeable shorter paths for anomalies. In fact, regions occupied by anomalies are low density regions, which result in a smaller number of partitions (shorter paths). Furthermore anomalies have, by definition, distinct feature-values and thus they are more likely to be separated early in the partitioning process.

• Then, another unsupervised algorithm named, Self-Organizing Map (SOM) is applied [START_REF] Kohonen | Self-organized formation of topologically correct feature maps[END_REF]. It acts as an extension of the k-means algorithm that preserves as much as possible the topological structure of the data. Moreover, SOM has an intrinsic distance between prototypes and their direct neighbours. This latter representation can validate the estimation of Isolation Forest if the results coincide. Finally, SOM gives a discretized representation of the input space, which categorize anomalies. The categorization helps to understand the origin of the problem.

In addition to the rare events detection task, we need to provide explanations of the different models. Moreover, important parameters must be discovered at a local level to figure out flaws in a single engine and at global level to discern origins of unexpected variations and inherent bias.

DATA ANALYSIS 2.1. Structure of the Data

The characteristics of the test bench data used in the analysis are as follows:

• 14 variables are chosen in an expert-manner by domain experts of the performance team of Safran. They are not generally interested in other variables and thus we limit ourselves to this subset of variables. However, in future works we will consider a larger set of variables. • 591 engines are observed.

• 4 stabilized points are considered.

A stabilized point, is a fixed level of performance for which all engines are tested and measurements acquired. They are ordered from the lowest to the highest level of performance. In the database, six are available but we do not consider the two first because they are taken at low engine speeds where there is a lot of variance in measures which makes them difficult to analyze. The measures of interest are listed below. Nine of them are numerical variables:

• FNIN1 : Thrust (FN: performance). • XN12R : LP spool speed (N1: fan speed). • XN25R : HP spool speed (N25: core speed). • WF36 : Fuel flow. • W2AR : Engine corrected air flow. • P3 : HP compressor discharge pressure. • T49C : LP turbine inlet temperature (EGT: exhausting gas temp).

• T3: HP compressor discharge temperature.

• P18QSC: Pressure section 18 normalized by standard conditions.

where HP and LP stand for High pressure and Low pressure respectively. The variables listed above are described in Figure 1. Moreover, four variables are categorical and represents test bench components:

• CELL : Bench.

• CNOZ : Primary nozzle.

• BMSN : Air nozzle.

• COWL : Nacelle.

Figure 1. Simplified diagram of a turbofan engine where the different measured variables are specified.

Bias in Data

An analysis of the data shows that production tests data are biased by test bench components. However, normalizing variables successively by each test bench component is an acceptable method only if test bench components are independent with each others. In this case we found out, using a pairwise ttest, that this type of normalization does not remove the bias in the data. Thus, it is preferred to keep the non-standardized data and to include the bench components to our models, which will be able to handle interactions between variables.

EXPERT KNOWLEDGE TO DEFINE ANOMALIES

Some pre-treatments are needed before detecting unusual case in the data. As said before, it is normal to have variance in production tests data. Those fluctuations do not necessarily represent unexpected behavior. In practice, the behavior of an engine is not defined with measurements taken independently, but it is defined between pairs of measurements. Thus, an engine has a normal level of functioning if it has a "constant ratio" between some defined pairs of variables across the stabilized points considered.

In other words we do not define an anomaly considering the observed values but we construct a new data set from the observed one where each variable is constructed considering the dependence between pairs of measurements.

In the new data set the dependence between the pairs of variables is of importance, and especially the evolution of these dependencies across the different stabilized points. At each stabilized point a physical equation describes the relationship for each pair of variables and reveals the expected behavior of engines.

Figure 2 shows an example of the physical equation (red line) for the pair of variables (FNIN1, W2AR) and the observed values for the set of engines at the fourth stabilized point.

The line gives one important information, that is the expected relationship between the thrust (FNIN1) and the mass flow rate (W2AR); which means that for a certain value of thrust we expect a certain mass flow rate. However, the equation of the functioning line (red line) is unknown and its estimation can be done with the help of a Robust-Principal component analysis (PCA) as detailed in [START_REF] Candès | Robust principal component analysis[END_REF]. PCA is a statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components. This transformation is defined in such a way that the first principal component has the largest possible variance, and each succeeding component in turn has the highest variance possible under the constraint that it is orthogonal to the preceding components. Therefore, the first axe of the PCA will model the functioning line, and the second axe will represent the space where the new variable will lie. The Robust PCA (RPCA) is an extension of PCA that is less sensible to extreme values, thus the slope of the functioning line will not depend on outliers and let them easier to detect.

Define Anomalies in Production Tests Data

For different stabilized points, we expect that an engine keeps a constant ratio between pairs of variables, i.e. their relation-ships do not change over the stabilized points. Therefore, we consider the engines values projected on the second axis obtained by the RPCA for each stabilized point. Then, we define a normal engine behavior as follow: "An engine has a normal behavior if it has similar projected values across the stabilized points".

Formally, let us consider an engine i ∈ {1, . . . , n}, a pair of variables j ∈ {1, . . . , J} and a stabilized points p ∈ {1, . . . P }. Let y j i,p be the projection of the engine i on the second axis of the RPCA for the pair of variables j at the stabilized point p. Then, the mean value over the stabilized points for an engine is

m j i = 1 P P p=1 y j i,p . (1)
The difference of an engine from its mean value for a stabilized point p is

x j i,p = y j i,p -m j i . ( 2
)
Thus given j, ∀p a new variable X j p = (x j 1,p , . . . , x j n,p ) T is created. The variable x j i,p represents the deviation of the engine i at a specific point p compared to its mean value m j i to the pair j. Thus, small values of x j i,p imply small deviations thus normal behavior, meanwhile large values lead to anomalies.

Nine pairs of measurements are defined in an expert manner are: (FNIN1, W2AR), (XN12R, W2AR), (P18QSC, W2AR) are the thrust, the LP spool speed and the pressure at section 18 given the engine corrected air flow. (FNIN1, WF36), (T49C, WF36) are the thrust and the exhausting gaz temperature given the fuel flow. The core speed given the fan speed, the pressure and the temperature at section 3 (HP) are also considered (XN25R, XN12R), (XN25R,T3), (XN25R, P3). Finally, the thrust function of the exhausting gaz is taken into consideration (FNIN1, T49C).

Each pair of variables is observed over four stabilized points, which give 36 new variables. In addition, the four test bench components, CELL, CNOZ, COWL and BMSN are kept in the set of variables that will be used in the model because the pre-treatment was not able to make the data independent from those ones. Finally, the models in following sections are applied on this set of 40 variables.

Note that, to keep an understanding of the new variables obtained from a pair, they are named as follow: "first variable name ___ second variable name ___ stabilized point ". For example, the variable created from T49C and FNIN1 on the fourth stabilized point will be called T49C_FNIN1_6 (the 4 stabilized point level are listed from 3 to 6).

ANOMALY DETECTION 4.1. Definition of the Method

Engines, in most of the cases, have solid and adequate measures. Thus, checking for unusual values can be seen as a statistical problem of outlier detection. For our purpose, a statistical method that is both efficient and interpretable is required. Density-based techniques are the most competitives and among these, Isolation Forest in [START_REF] Liu | Isolation forest[END_REF] showed best results on various studies [START_REF] Goix | Apprentissage automatique et extrêmes pour la détection d'anomalies[END_REF]. We employ this method on the new data to detect anomalies.

Isolation Forest is similar in principle to Random Forest [START_REF] Breiman | Random forests[END_REF] and is built on the basis of decision trees. It identifies anomalies or outliers rather than profiling normal data points. Isolation Forest isolates observations by randomly selecting a feature and then randomly selecting a split value between the maximum and minimum values of that selected feature. Then, if an observation lies in a high-density region, the probability to isolate it is small because the values of the splits must be very close. On the other hand, if an observation lies in a low-density region, then many values of splits can isolate it, thus it has a higher probability to be isolated by a random split. Random partitioning produces noticeably shorter paths for anomalies. When a forest of random trees collectively produces shorter path lengths for particular samples, they are highly likely to be anomalies. Let h t (x) the path length of x in the tree t and h(x) = T t=1 h t (x) the average value of h(x) over the trees with T is the total number of trees in the forest. The number of splits required to isolate an observation is influenced by the number of samples n in the data. To account for this a normalized anomaly score, relying on a property of Binary Search Trees (BST)) [START_REF] Liu | Isolation forest[END_REF], is defined as

f (x, n) = 2 -h(x) c(n) , (3)
with c(n) defined as

c(n) =      2H(n -1) - 2n -1 n for n > 2, 1 for n = 2, 0 otherwise, ( 4
)
where n is the size of data set and H is the harmonic number.

The value of c(n) above represents the average of h(x) given n, so we can use it to normalise h(x) and get an estimation of the anomaly score for a given instance x. Note that f ∈ [0; 1] with value closer to 1 indicates that the observation is more likely to be an anomaly.

Once the Isolation Forest has been applied to the data, an anomaly score is estimated for each engine. This anomaly score is pointless if it cannot be completely understood by domain experts.

For complex models, such as ensemble methods, deep networks or Isolation Forest, we cannot use the original model as its own best explanation because it is not easily understandable. Instead, we must use a simpler explanation model, which we define as any interpretable approximation of the original model. We would like to have an average explanation of the model as well as explanation of single prediction and this is called as local explainability [START_REF] Guidotti | A survey of methods for explaining black box models[END_REF] which is also known as "post-hoc" explainability. In, Doshi-Velez & Kim (2017) they assert that a useful local explanation should answer the following questions: What were the main factors in the decision? Would changing a certain factor have changed the decision? Why did two similar-looking cases get different decisions, or vice versa? More precisely, we would like to understand how variables contributed to the score of a single engine as well as how variables contributed in average. In addition, understand whether the contribution of a variable have a positive impact or a negative impact on the score, or in other words if a variable helps making an observation more normal or more abnormal. In this aim, Shapley values will be used.

Model Interpretability with Shapley Values

Shapley values have attracted a great deal of attention in recent years in the field of interpretability, which has been originally discussed in game theory [START_REF] Shapley | A value for n-person games[END_REF] Shapley values measure features importance for models in the presence of interaction between variables. This method requires retraining the model on all feature subsets S ⊆ F , where F is the set of all features. As described in [START_REF] Lundberg | A unified approach to interpreting model predictions[END_REF], it assigns an importance value to each feature that represents the effect on the model prediction of including that feature. To compute this effect, a model f S∪j is trained with that feature present, and another model f S is trained with the feature excluded. Then, predictions from the two models are compared on the current input f S∪j (x S∪j ) -f S (x S ), where x S represents the values of the input features in the set S.

Since the effect of withholding a feature depends on other features in the model, the preceding differences are computed for all possible subsets S ⊆ F \j. The Shapley values are then computed and used as feature attributions. They are a weighted average of all possible differences:

φ j = S⊆F \j |S|!(|F | -|S| -1)! |F |! f S∪j (x S∪j ) -f S (x S ) , (5)
At first glance, the above equation seems ridiculously complicated, but it can be easily explained in one sentence: "The contribution of a feature j is the mean difference between a model trained on a subset of variables S with j and a model trained on the same subset S without j, and this is done for all the possible subsets of variables". All possible sets of feature values have to be evaluated with and without the j-th feature to calculate the exact Shapley value. For more than a few features, the exact solution to this problem becomes problematic as the number of possible coalitions exponentially increases as more features are added. In [START_REF] Štrumbelj | Explaining prediction models and individual predictions with feature contributions[END_REF] an approximation with Monte-Carlo sampling is proposed

φj = 1 M M m=1 f (x m +j ) -f (x m -j ) , (6)
where f (x m +j ) is the prediction for x, but with a random number of feature values replaced by feature values from a random data point z, except for the respective value of feature j. The x-vector x m -j is almost identical to x m +j , but the value x m j is also taken from the sampled z. Each of these M new instances is an "artificial object" assembled from two instances. In our case, f (x) is the anomaly score predicted by the Isolation Forest. If φj is positive the value of the feature j increase the anomaly score, and it decreases if φj is negative.

RESULTS ON THE DATA 5.1. Average Shapley Values

The Table 1 gives the average Shapley values by variable on all observations. It provides a nice interpretation of the effect of each variable on the anomaly score. A variable j with a φ j value close to 0 will not affect the output of the model and thus the variable is not important to detect anomalies. A high absolute φ j value points out that the variable plays an important role in model estimates. The sign of the contribution gives an additional information on the effect of the variable. A positive contribution denotes that the variable helped to increase the estimated anomaly score w.r.t. the average anomaly score whereas a negative contribution decreases it. Therefore, if in average a variable has a negative φ j , it means that it does not globally contribute to make an engine significantly different from others.

Single Prediction Explanation

Averaged explanation are useful and give insights but they are not sufficient. When an engine has a high estimated anomaly score, domain experts would like to understand which variables are responsible for this score. As an example, the engine 41 is observed. Shapley values are applied to explain how variables contributed to the score. On Figure 3, the average score of anomaly for engines is 0.40 and the engine 41 has an anomaly score of 0.47, which is significantly larger. This difference is decomposed variable by variable. The x-axis, φ, gives the weight of the contribution. On the y-axis, the engine values for each variable are displayed, ordered by decreasing importance. A value of 0 on the x-axis indicates that the variable does not play any role in the estimation of the score. Note that, most important variables for this engine correspond to the highest deviations x j 41,p obtained by RPCA. However, due to the possible high-order interaction between variables, a complex model was needed to assess a good estimation of the anomaly score. and they can control the validity of the results. Figure 4 shows the densities of variables with the highest (WF36_T49C_5) and the lowest (T3_XN25R_5) φ detected for engine 41. As expected, the value x j 41,5 for j =WF36_T49C, is far from 0 and isolated in a low density area, which means that it has a really different behavior over the four stabilized points. On the other hand, for j =T3_XN25R, x j 41,5 is close to 0 and the engine has a consistent behavior. using the statistical methodology is presented in Figure 6. Isolation Forest helps domain experts to identify few engines and Shapley values help them to focus on some specific measures. Then a complete inspection of the engine can be done before validating the production test.

Domain experts have access to the contribution of the variables

ANOMALY

CATEGORIZATION USING SELF-ORGANIZING MAP 6.1. Definition of the Method A SOM is a type of artificial neural network that is trained using unsupervised learning to produce a low-dimensional (typically two-dimensional) discretized representation of the input space of the training samples, called a map. Each unit of the map corresponds to a prototype vector in the original high-dimensional space, and new data points are projected on the map by finding the closest prototype vector w.r.t. euclidean distance [START_REF] Kohonen | Self-organized formation of topologically correct feature maps[END_REF]; [START_REF] Olteanu | On-line relational and multiple relational som[END_REF]. Self-organizing maps have been used for aircraft engine fleet monitoring in [START_REF] Cottrell | Fault prediction in aircraft engines using self-organizing maps[END_REF]; Côme et al. (2010b,a); [START_REF] Forest | A generic and scalable pipeline for large-scale analytics of continuous aircraft engine data[END_REF] and to classify transient flight phases [START_REF] Faure | Using self-organizing maps for clustering anc labelling aircraft engine data phases[END_REF]. No specific study has been yet conducted on using SOM to validate and categorize anomalies and especially on production tests data.

SOM has both intrinsic distances between clusters and nice two-dimensional visualization, which make it a good candidate. Nonetheless, clusters are still in high-dimensions and methods such as Shapley value are not tractable in this situation. Therefore, in the next section, before modelling a SOM, specific clustering algorithm for variable selection will be used.

6.2. Choose a Subset of Variables with the Help of Group-Sparse Weighted K-means Group-sparse weighted K-means generalizes the sparse weighted K-means algorithm for numerical variables in [START_REF] Witten | A framework for feature selection in clustering[END_REF], by using the group regularization framework. Suppose that the numerical matrix of data X is described by p features that are divided into L priorly known distinct groups, such that X = X 1 |...|X L , with X ∈ R n×p , p being the size of group , and p 1 + ... + p L = p.

In presence of group data, we would like to discriminate groups of variables X by using a specific L 1 -group penalty, which has been already used in the regression framework [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. This allows us to select variables by group, forcing the model to select or discriminate the entire group. As described in [START_REF] Chavent | Sparse k-means for mixed data via group-sparse clustering[END_REF], the between-class variance of each variable is multiplied by a weight and a parameter λ penalizes the weights. The latter discriminates the groups of variables with the lowest between-class variance. There is a clustering solution (groups weights and clusters) for each fixed λ. The regularization path (clustering solution given lambda) is computed at a grid of values for the penalty factor λ, covering the entire range, from a model with all the groups included to a model with only one group. The optimization procedure is quite straightforward. The algorithm is optimized in an iterative fashion: first the K-means algorithm is performed on the weighted space of features, then the partition is held fixed and the weights are updated. This iterative procedure is continued until a (local) minimum is reached.

In this context, groups are clearly formed by the variables over the stabilized points. For example, T3_XN25R_3, T3_-XN25R_4, T3_XN25R_5 and T3_XN25R_6 belong to the same group. Hence, there is 9 groups of variables, and we would like to know which are the most discriminative for clustering. Moreover, Table 1 shows that test bench components were not significant to detect unusual behavior. Shapley values attribute them in average a negative contribution, which implies that they are not useful to model unexpected behavior. Thus, these variables are not considered in the analysis. Furthermore, the number of clusters is set to five and was found with the Silhouette method [START_REF] Rousseeuw | Silhouettes: A graphical aid to the interpretation and validation of cluster analysis[END_REF].

In Figure 7 we provide the path of groups' weights against the λ sequence. Weights of groups, on the y-axis, are represented for each λ. The most important groups are 7, 4 and 5 which are respectively the groups P3_XN25R, T49C_FNIN1 and WF36_FNIN1.

Figure 7. Groups weights for each value of λ. The vertical line represents the selected clustering solution where three groups of variables have non-zero weights.

In Figure 8 we see a big gap in terms of explained variance before λ = 0.6 If one give a closer look, the analysis points out that a subset of variables (3 groups) will give similar clustering, in terms of explained variance, to the one with all the variables included (see Figure 7 and Figure 8). We choose the clustering obtained for value of λ represented by the vertical line. This value of λ allows to select 3 groups with high explained variance. Higher value of λ lead to have a more similar solution to the one after the gap in Figure 8 in terms of weights, which seems to be a bad clustering solution. Therefore, the chosen value of λ seems to be a good tradeoff between interpretability and performance. The weights obtained by groups are w T49C_FNIN1 = 0.67, w WF36_FNIN1 = 0.62 and w P3_XN25R = 0.42. This subspace of 12 variables will be used to represent the data with the help of the SOM algorithm.

Figure 8. The ratio of between-sum of squares over the total sum of squares (BSS/TSS: explained variance) for each value of λ. The vertical line represents the selected clustering solution which has an explained variance that is close to the one with the full set of variables.

Data Representation with SOM

The SOM allows us to have access to several different visualizations. First, we plot the distances between prototypes (Figure 9), which gives a representation of the grid where the colors represents the mean distance to the neighbor prototypes. The color scale goes from blue to purple, where purple indicates a large distance. In Figure 10 the repartition of the engines on the map is provided. Finally the Figure 11 is the SOM map where colors indicate the mean anomaly score of engines estimated with the Isolation Forest by clusters. The color scale goes from yellow to red, where red indicates a higher anomaly score. It is interesting to note that the engines with high anomaly score are distributed on the border of the map. The representation is very similar to the one given by the smooth distances between prototypes on Figure 9 which shows that the two methods agree. In addition, the map allows a categorization of the anomalies. In Figure 11, Super-clusters are identified thanks to the mean anomaly score of the prototypes and their proximities on the SOM map. Two map edges are thus identified as super-clusters. The comparison of these clusters of anomalies with the rest of the population will allow us to understand the discriminative variables.

The ANOVA method is applied to test significant difference between clusters (Figure 12). ANOVA provides a statistical test of whether two or more population means are equal. The results show that, for cluster 2, the group of variables P3_-XN25R seems to be important. On the other hand, for cluster 1, the group of variables WF36_FNIN1 may explain their unusual behavior. The ANOVA shows that the set of explaining variables has been reduced to only one group. For the sake of readability, we show the ANOVA test for only two variables at the stabilized point 3, but for the other stabilized points results are similar since they are all linked by construction. Moreover, the other groups of variables are not significantly different over the three clusters.

CONCLUSIONS

In this work, statistical methods in the context of rare event detection in production tests data demonstrated high degree of efficiency and interpretability in either local (one specific engine) or global level (groups of engines). We propose a multi-scale model, giving a hierarchy in the information allowing the experts to better understand flaws on a particular engine but also allowing them to detect more general problems. We apply two different methods: i) Isolation Forest to estimate anomaliy score and Shapley values to interpretate them; ii) SOM on a subset of variables obtained by group-sparse weighted K-means. Both methods provide similar results: the engines detected as anomalies with the Isolation Forest coincide to the engines that have the largest distances estimated with SOM.

Moreover, an other contribution of this work is the use of SOM to validate anomalies detection methods. On the contrary of Isolation Forest, SOM provides visualizations and categorizations of the anomalies which gives additional information to better understand and verify the estimated anomalies.

Explainability in unsupervised learning is a new field that needs to be explored, and this work is a step forward in this direction. Some further investigations are needed in both theoretical and applied domains. In future works, we plan to explore those points and also we will develop a method to transform the anomaly score in a binary score which will help to domain experts in their decisions.

Figure 2 .

 2 Figure 2. Representation of the functioning line for the pair of variables (FNIN1, W2AR) on the stabilized point four. The value of each engine is represented by the colors in the orthogonal space estimated with the RPCA.

  and recently applied to statistics Štrumbelj & Kononenko (2014); Owen & Prieur (2017); Iooss & Prieur (2017); Lundberg & Lee (2017). Moreover, in the context of anomaly detection, few results have already been reported Antwarg et al. (2019); Giurgiu & Schumann (2019); Takeishi (2019); Takeishi & Kawahara (2020). These results have confirmed the usefullness of the Shapley value for anomaly interpretation. As described in their works, we will adopt general techniques in defining and computing the Shapley values derived from supervised learning.

Figure 3 .

 3 Figure 3. Shapley plot of the engine 41. 0.40 is the average anomaly score and 0.47 is the predicted score of the engine 41. Feature value with positive φ increase the score from 0.40 to 0.47, and negative value of φ decrease the anomaly score of the engine.

Figure 4 .

 4 Figure 4. The density of the variable WF36_T49C_5. Engine 41 is isolated in a low density area which explains why this variable have a high contribution to increase the anomaly score.

Figure 6 .

 6 Figure 6. Diagram of the process of production tests validation using Isolation Forest and Shapley values. The statistical methodology helps to highlight specific engines and measures where further analyses are needed.

Figure 9 .

 9 Figure 9. Smooth distances between prototypes. The background colors indicate the distances between neighboring prototypes where pink corresponds to larges distances.

Figure 10 .

 10 Figure 10. Repartition of the engines on the map.

Figure 11 .

 11 Figure 11. SOM map of engines where the colors represents the anomaly score estimated with Isolation Forest averaged by clusters. Color scale goes from yellow to red, where red corresponds to a higher score.

Figure 12 .

 12 Figure 12. Boxplot of ANOVA test comparing two different area of the map Figure 9. For the two variables considered, only one cluster (the pink one) is truly different from the overall population of engines, implying that the two clusters of anomalies can be explained by different subsets of variables.

Table 1 .

 1 Average Shapley value by variables. Positive φ j values indicate that the variable tends to increase the anomaly score while negative one indicates the opposite.

	variable	φ
	WF36_FNIN1_4	1.7839
	WF36_T49C_4	1.4149
	P3_XN25R_3	1.3275
	T49C_FNIN1_6	1.0992
	XN12R_XN25R_5	0.8475
	WF36_T49C_3	0.7794
	WF36_FNIN1_3	0.7782
	P3_XN25R_6	0.7207
	W2AR_XN12R_3	0.5513
	WF36_T49C_6	0.4814
	W2AR_XN12R_4	0.4696
	XN12R_XN25R_6	0.4378
	XN12R_XN25R_4	0.3720
	WF36_T49C_5	0.2933
	P3_XN25R_5	0.2284
	W2AR_XN12R_6	0.1843
	WF36_FNIN1_6	0.1367
	CNOZ	0.1318
	T3_XN25R_4	0.0514
	CELL	0.0487
	T49C_FNIN1_4	0.0212
	WF36_FNIN1_5	-0.0029
	W2AR_XN12R_5	-0.0201
	XN12R_XN25R_3	-0.0446
	P3_XN25R_4	-0.0780
	W2AR_FNIN1_6	-0.1006
	T49C_FNIN1_5	-0.1676
	W2AR_FNIN1_4	-0.2927
	COWL	-0.4981
	W2AR_P18QSC_5 -0.5035
	BMSN	-0.5225
	W2AR_P18QSC_3 -0.6202
	W2AR_P18QSC_4 -0.6486
	W2AR_FNIN1_3	-0.7632
	W2AR_FNIN1_5	-0.7816
	W2AR_P18QSC_6 -0.7929
	T3_XN25R_5	-0.7946
	T3_XN25R_6	-0.8460
	T49C_FNIN1_3	-1.0517
	T3_XN25R_3	-1.1104

ACKNOWLEDGEMENT

This work was supported by the French National Agency for Research and Technology (ANRT) and Safran Aircraft Engines (Safran group).


Francisco Rony

Gomes Barroso

Vaneicia Dos

Santos Gomes

Carlos Eduardo Carvalho

Marie Pierre Ledru

Charly Favier

Francisca Soares

Laurent Bremond

Francisco Rony Gomes Barroso

Phytoliths from soil surfaces and water reservoirs of the Brazilian semi-arid Caatinga

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Phytoliths, pollen, and charcoals are commonly used proxies for paleoenvironmental and paleoecological studies in the paleobotany field. Phytoliths are very interesting proxies for local records in dry environments because they resist oxidizing conditions, unlike pollen, which is preserved only under anoxic conditions, being more useful for regional analysis [START_REF] Barboni | Phytoliths as palaeoenvironmental indicators, west side middle awash walley, Ethiopia[END_REF][START_REF] Coe | Reconstituições da vegetação e inferências de paleoclimas através da utilização dos indicadores fitólitos e isótopos de carbonoexemplos de estudos no Brasil[END_REF][START_REF] Gosh | Phytolith assemblages of grasses from the Sunderbans, India and their implications for the reconstruction of deltaic environments[END_REF]IPCC, 2011;[START_REF] Piperno | The status of phytolith analyses in the American tropics[END_REF]. The production and record of charcoals depend on the vegetation flammability and can feature the regional or local vegetation, depending on the charcoal type, shape, and size [START_REF] Scott | Charcoal: taphonomy and significance in geology, botany and archaeology[END_REF]. For the use of these proxies, paleoenvironmental studies require special environments with active sedimentation conditions and rare erosive events, as they preserve contiguous layers of sediment, for example, large lakes or peat bogs [START_REF] Neumann | Ambientes lacustres[END_REF][START_REF] Suguio | Recent progress in Quaternary geology of Brazil[END_REF].

In the state of Ceará, the presence of a perennial natural lake is not common. Fortunately, after a series of droughts (from 1877 to the present date) in Ceará, federal-state measures were adopted to construct reservoirs [START_REF] Campos | Secas e políticas públicas no semiárido: ideias, pensadores e períodos[END_REF]. The sediments deposited and accumulated over the years in these reservoirs can contain proxies, such as phytoliths, pollen, and charcoals, as in recent sediments of natural lakes [START_REF] Aleman | Tracking land-cover changes with sedimentary charcoal in the Afrotropics[END_REF]. Therefore, the application of paleoecological methods to study the recent sediments of artificial reservoirs can be promising to assess recent vegetation changes resulting from climate and land-use changes. In this sense, studies of soil and reservoir surface phytoliths can provide data on how different vegetation physiognomies are spatially represented and how the reservoir sediment composition varies in tens or hundreds of years. Thus, it will be possible to understand the effect of degradation and deforestation through phytolith composition in chronosequences of few decades in two reservoirs of the Brazilian semi-arid.

In the present study, we used phytoliths as a proxy of tree cover and herbaceous composition. These microscopic opal silica particles, formed in the cells of different plant organs, are released by leaf litter deposition on the soil surface or after the plant's death [START_REF] Piperno | The Silica Bodies of Tropical American Grasses: Morphology, Taxonomy, and Implications for Grass Systematics and Fossil Phytolith Identification[END_REF]. Phytoliths are preserved both in soil and soaked marshes for hundreds and even thousands of years. Although phytoliths are diagnoses of specific genera or families, some have no taxonomic significance or, in other cases, different families can produce the same type of phytolith (redundancy) [START_REF] Piperno | Phytolith Analysis: an Archaeological and Geological Perspective[END_REF][START_REF] Rovner | Potential of opal phytolith for use in paleoecological reconstruction[END_REF]. In some plants, especially in the Poaceae family, there is a diversity of forms produced by the same species (multiplicity).

Poaceae can produce short silica grass cells (GSSC) [START_REF] Mulholland | A morphological classification of grass silica-bodies[END_REF]. The subfamily Panicoideae produces bilobate, polylobate, and cross morphotypes, which are grasses adapted to hot and humid climates. The saddle type generally occurs in the subfamily Chloridoideae, C4 grasses adapted to dry climates [START_REF] Twiss | Morphological classification of grass phytolith[END_REF]Mulholland, 1989;[START_REF] Fredlund | Modern phytolith assemblages from the north American great plains[END_REF][START_REF] Kondo | Opal Phytoliths of New Zealand[END_REF]. The presence of woody dicotyledonous may be associated with the production of globular granulate, sclereid, and faceted block types [START_REF] Scurfield | Silica in woody stems[END_REF][START_REF] Kondo | Opal Phytoliths of New Zealand[END_REF][START_REF] Alexandre | Phytoliths: indicators of grassland dynamics during the late Holocene in intertropical Africa[END_REF][START_REF] Runge | The opal phytolith inventory of soils in central Africa-quantities, shapes, classification, and spectra[END_REF][START_REF] Bremond | A phytolith index as a proxy of tree cover density in tropical areas: calibration with Leaf Area Index along a forest-savanna transect in southeastern Cameroon[END_REF][START_REF] Neumann | Woods of the Sahara and the Sahel[END_REF]. Cyperaceae can produce papillae, and Arecaceae may originate globular echinate [START_REF] Kondo | Opal Phytoliths of New Zealand[END_REF]. In general, Caatinga plants produce a large number of phytoliths. Some families, such as Arecaceae, Bromeliaceae, Chrysobalanaceae, Malvaceae, Anacardiaceae, and Erythroxylaceae, are major producers [START_REF] Coe | Caracterização de fitólitos de plantas e assembleias modernas de solo da Caatinga como referência para reconstituições paleoambientais[END_REF].

Studies on phytoliths in Brazil are still incipient [START_REF] Luz | Estágio Atual do Conhecimento sobre Fitólitos no Brasil[END_REF], but there is relevant information about different Brazilian biomes. In peat bogs of a mixed rainforest in the state of Paraná, southern Brazil, [START_REF] Rasbold | Reconstrução paleoambiental de um depósito sedimentar por análises multiproxy, turvo, Estado do Paraná, Brasil[END_REF] recognized drier climatic conditions 14,500 years ago, with a gradual increase in humidity starting 6000 years ago. In a restinga in Maricá, in the state of Rio de Janeiro, southeastern Brazil, [START_REF] Santos | Opal phytolith and isotopic studies of "Restinga" communities of Maricá[END_REF] found little production of phytoliths in soil samples from six different plant communities, except for the families Poaceae, Cyperaceae, and Arecaceae, which occurred in large quantities and diversity. In the Caatinga biome, in the state of Ceará, northeastern Brazil, [START_REF] Ricardo | Reference collection of plant phytoliths from the Caatinga biome, Northeast Brazil[END_REF] found large phytolith assemblages extracted from plants, in which the most prevalent were tracheids, followed by polyhedral, globular granulate, elongate, and trichomes (lanceolate hair cells). [START_REF] Dias | Reconstituição paleoambiental de dunas vegetadas na Caatinga, em Aracati, Ceará, através de biomineralizações de sílica[END_REF] have studied phytoliths from soil profile, also in State of Ceará, with a predominance of bulliform, acicular and globular types and the results showed the tree density increased in the period between 2300 and 700 years cal BP.

Studies carried out in the Caatinga biome in Ceará [START_REF] Coe | Caracterização de fitólitos de plantas e assembleias modernas de solo da Caatinga como referência para reconstituições paleoambientais[END_REF][START_REF] Dias | Reconstituição paleoambiental de dunas vegetadas na Caatinga, em Aracati, Ceará, através de biomineralizações de sílica[END_REF][START_REF] Ricardo | Reference collection of plant phytoliths from the Caatinga biome, Northeast Brazil[END_REF] did not present paleoecology data with phytoliths from cores collected in reservoirs. Our study reconstructs paleoecological and paleoenvironmental data and shows possible evidence of anthropogenic action in vegetation cover. It is essential to understand the temporal reconstruction of the Caatinga biome that, since the arrival of the Portuguese, has been threatened by agricultural activities, soil erosion, ecosystem degradation, and alteration in species relative abundance [START_REF] Pompeu-Sobrinho | Povoamento Do Nordeste Brasileiro[END_REF]. These anthropogenic pressures changed the vegetation cover, transforming closed tree vegetations into more open physiognomies, with rich herbaceous pastures suitable to livestock farming [START_REF] Pompeu-Sobrinho | Povoamento Do Nordeste Brasileiro[END_REF]. Our goals are a) to study the original flora of the Caatinga, recorded through phytolith assemblages in three different vegetation covers, b) to retrieve the land-use history from phytoliths present in sediments of water reservoirs, and c) to reconstitute the vegetation cover dynamics.

Study area

The Caatinga Phytogeographical Domain (CPD) is a broad Brazilian natural region located in the northeastern corner of the South American continent [START_REF] Ab'sáber | Os domínios de natureza no Brasil: potencialidades paisagísticas[END_REF]IBGE, 2004). The Caatinga sensu stricto is the main vegetation type of the CPD, comprising approximately 800,000 km2. It is also widespread in the state of Ceará (Fig. 1A). According to [START_REF] Veloso | Classificação da vegetação brasileira, adaptada a um sistema universal[END_REF], "Caatinga" is a local name, and this vegetation can be called savanna.

The Caatinga vegetation comprises a deciduous, thorny, and succulent biome with open and grey aspect during the dry season, when almost all trees and shrubs shed their leaves. The Caatinga is frequently mentioned in the literature as widely variable in composition as well asin physiognomy, in response to climate and topography heterogeneity, soil types, and human influence [START_REF] Andrade-Lima | The Caatingas dominium[END_REF][START_REF] Cabrera | Biogeografia de America Latina[END_REF][START_REF] Sampaio | Overview of the Brazilian Caatinga[END_REF][START_REF] Fernandes | Fitogeografia Brasileira[END_REF][START_REF] Moro | A phytogeographical metaanalyses of the semiarid Caatinga Domain in Brazil[END_REF].

In Ceará, the predominant climatic characteristics are high average air temperatures during most of the year [START_REF] Souza | Condições geo-ambientais do semi-árido brasileiro[END_REF]. Rains are considerably irregular throughout the year, and the rainfall concentrates in the first semester [START_REF] Zanella | As características climáticas e os recursos hídricos do Ceará[END_REF]. Cities have rainfall values ranging between 201.80 and 1700 mm, but usually within 400.01 and 600 mm [START_REF]Palaeoclimate proxy indicators[END_REF]. The Caatinga soils are rocky and shallow, in general, with different characteristics, even within short distances (Ab'Sáber, 1974;[START_REF] Sampaio | Overview of the Brazilian Caatinga[END_REF].

The Caatinga comprises several physiognomies. In Ceará, the vegetation cover can be divided into eleven different types (Fundação, 1994). The Phytoecological Units map shows the area of occurrence of the different vegetation types found in Ceará (Fig. 1). In general, the original vegetation was altered as a result of a disorderly and predatory human action [START_REF] Antongiovanni | Chronic anthropogenic disturbance on Caatinga dry forest fragments[END_REF].

In a forest physiognomy, the treetops touch each other during the rainy season, the favorable season for plants in the semi-arid region. Schinopsis brasiliensis Engler (Anacardiaceae) and Astronium urundeuva (M. Allemão) Engl. (Anacardiaceae) are some of the species present in the Caatinga Forest. The degradation of the Caatinga Forest results in the appearance of the Shrubby Caatinga characterized by the smaller size and twisted and whitish stems. The greater or lesser density of individual members of the community results in Dense Shrubby Caatinga or Open Shrubby Caatinga (Fundação, 1994). According to [START_REF] Andrade-Lima | The Caatingas dominium[END_REF], the Open Shrubby Caatinga occupies small, dispersed areas in the Caatinga, with predominantly small species, including some herbs in open spaces. The Dense Shrubby Caatinga is the vegetation type that occupies the largest area nowadays. The most frequent species in Shrubby Caatingas are Mimosa tenuiflora (Willd.) Poir, Cenostigma bractosum (Tul.) Gagnon & G.P. Lewis, Mimosa caesalpinifolia Benth (Fabaceae), Croton sonderianus Müll. Arg. (Euphorbiaceae), and Cereus jamacaru DC (Cactacea) (IPLANCE, 1995). Those three vegetations comprise the "Crystalline Caatinga", a unique phytoecological unit with a cohesive flora [START_REF] Moro | Vegetação, unidades fitoecológicas e diversidade paisagística do estado do Ceará[END_REF].

There are larger forests along semi-arid river courses (riparian forests) than in the surrounding vegetation. The riparian forest and "Carnaubais" cover the interior and follow the Acaraú River. Copernicia prunifera (Mill.) H.E. Moore (Arecaceae) stands out in riparian forests. The typical soils of these areas are sandy and can reach great depths, allowing the development of riparian forests associated with Copernicia prunifera (Carnaúba Palm), forming "Carnaubais". The riparian forest associated with "Carnaubais" occurs in semi-arid areas from the interior of Ceará to sub-humid coastal areas associated with riverbeds, lagoons, or periodically flooded regions [START_REF] Moro | Vegetação, unidades fitoecológicas e diversidade paisagística do estado do Ceará[END_REF].

Material and methods

Study sites

We analyzed phytoliths from three vegetation types in the Crystalline Caatinga: open shrubby, dense shrubby, and forest vegetation. There were two different analyses: 12 soil samples and 2 cores (in different places). We took 12 soil samples from areas with different vegetation cover, in eight cities of Ceará (Table 1). We did not include riparian forests in our study, but they surrounded the Open Shrubby Caatinga cover. We have collected sediments in two reservoirs for core analysis in the Open Shrubby Caatinga. The Acaraú-Mirim and Araras dam reservoirs, in Massapê and Varjota, respectively, are part of the Acaraú River Basin in northern Ceará. The Araras reservoir is near the Sierra Ibiapada sedimentary plateau at the east, whereas the Acaraú-Mirim reservoir is next to high granitic mountains known as Sierra Meruoca. The criteria for choosing the 12 points from soil surfaces in different vegetation covers was the minimum count of 200 phytolith morphotypes per sample. For the chosen reservoirs, the criteria was the period of construction, being the both oldest in the State of Ceará.

Sampling and sediment dating

During the sediment sampling for phytolith extraction, three field expeditions were carried out in the state of Ceará as part of the activities of the project "Evolution of biodiversity loss in areas under degradation processes" (CNPq process # 400,890/2014-3), initiated in 2014, based on a partnership between the Federal University of Ceará and the University of Montpellier, France. The sediments were collected with a Kajak-Brinkhurst sampler (KB) in the middle of the reservoir (2 replicates per reservoir). The KB sampling method allows collecting a few tens of centimeters of unaltered sediments at the water-sediment interface, representing the last decades. Two cores were collected in two older dams in Ceará: Aracaru-Mirim (1907) and Açude Araras (1958). Samples were collected contiguously in 1 cm deep layers for phytolith analyses. A chronological framework was established to date sediments using short-lived radionuclides measurements. Age-depth models are fully described in [START_REF] Ledru | When archives are missing, deciphering the effects of public policies and climate variability on the brazilian semi-arid region using sediment core studies[END_REF]. The base of lake sediment dated to 1963 ± 8.5 yr AD, in the Acaraú-Mirim reservoir, and 1959.5 ± 11.5 yr AD (facies 1 and 2) in the Araras reservoir. The topmost sample is described in soil analysis.

Phytolith extraction and counting

We extracted phytoliths from 3 g sediment samples using the following protocol [START_REF] Aleman | Tracking land-cover changes with sedimentary charcoal in the Afrotropics[END_REF]: (1) deflocculation of the sediment with NaPO3, (2) dissolution of carbonates using HCl, (3) oxidation of organic matter with H2O2, (4) removal of clay by gravity, (5) drying of the residue with ethanol, and (6) densimetric separation of phytoliths in a dense solution of ZnBr2 (d > 2.3). We mounted the recovered fraction on microscope slides in immersion oil 518 (Carl Zeiss™ Immersol™) for 3D observation and counted phytoliths at a 600X magnification. We deposited the slides in the Laboratory of Phytogeography of the Federal University of Ceará (UFC). We counted a minimum of 200 Grass Silica Short Cells (GSSC). We made phytolytic diagrams using the PSIMPOLL program [START_REF] Bennett | Documentation for Psimpoll 4.25 and Pscomb 1.03: C Programs for Plotting Pollen Diagrams and Analysing Pollen Data[END_REF].

Phytolith classification

We described phytoliths according to their three-dimensional shape and classified them following the international phytolith nomenclature code [START_REF] Madella | International code for phytolith nomenclature 1.0[END_REF]. We calculated three phytolith indices for each sample. 1. Aridity index (Iph) indicates the dominance of shortgrass or tall-grass savannas based on the proportion of Chloridoideae (saddle type) to the sum of Chloridoideae and Panicoideae (saddle, cross, and bilobate) phytoliths [START_REF] Bremond | Definition of grassland biomes from phytoliths in West Africa[END_REF]. 2. D/P index (ratio of woody dicotyledon to Poaceae phytoliths) was also calculated as the ratio of globular granulate and decorated phytoliths to the GSSC [START_REF] Bremond | Definition of grassland biomes from phytoliths in West Africa[END_REF]. 3. Palm Tree Cover index (Pa/P), which is the ratio between the phytolith type characteristic of Arecaceae (Globular echinate) and the sum of phytolith types of Poaceae (Bulliform cuneiform + Bulliform parallelepipedal + Short cells + Acicular) [START_REF] Coe | Caracterização de fitólitos de plantas e assembleias modernas de solo da Caatinga como referência para reconstituições paleoambientais[END_REF].

Results

Phytoliths of soil surfaces

We identified 27 types of phytoliths (Figs. 2 and3). The saddle was the most frequent (39%), followed by other morphotypes, such as the globular granulate (15%), bilobate (14%), globular echinate (11%), and papillae (3%) (Fig. 3A). Sample 12 showed the highest number of phytoliths and sample 9, the lowest. They both came from the Caatinga Forest.

In the Open Shrubby Caatinga, samples 1, 4, 5, and 6 showed a high frequency of saddle, followed by bilobate types. Samples 2 and 3 showed mostly globular echinate, saddle, and bilobate types. In the Dense Shrubby Caatinga, sample 7 showed mainly bilobates, followed by saddle and papillae types, whereas sample 8 had a predominance of saddle and similar values of bilobate, globular echinate, and globular granulate types. Finally, in the Caatinga Forest, sample 12 was the only one with a high frequency of globular granulate (57%). Samples 9, 10, and 11 showed mainly saddle, followed by bilobate types.

We calculated phytolith indices for each vegetation cover (Fig. 3A). The mean D/P index was 0.2; the highest value (2.7) occurred in sample 12 (Caatinga Forest) and the lowest (0.005) in sample 8 (Open Shrubby Caatinga). The value was null for samples 4 (OSC), 9 (CF), and 11 (FC). The mean Iph index (%) was 72%; the highest value (90%) occurred in sample 4 (OSC) and the We identified 27 phytolith types in the Acaraú Mirim core (Fig. 3B). The predominant types were globular echinate (41%), followed by bilobate (14%), saddle (13%), and globular granulate (13%). The level with the highest number of phytoliths corresponded to the year 2000, and the lowest to 1970. There was not enough phytolith for counting corresponding to the year 2009.

Phytoliths from reservoir sediments

The globular echinate type was dominant in nine out of ten samples, only in the sample of the year 1996, globular granulate showed the highest number. The presence of bilobate and saddle types has changed along the core. Between 1965 and 1978, there was a continuous decrease in globular echinate and an increase in globular granulate phytoliths. In 1991, globular echinate phytoliths were highest in the core, and we found the lowest level of saddle and globular granulate types.

The mean D/P index was 0.4, with the highest value (1.0) found in 1996, and the lowest value (0.2) in 1991, 1983, 1970, and 1965. The mean Iph (%) index was 50%, with the highest value (57%) found in the year 2000, and the lowest (42%), in 1983.

We identified 17 types in the Araras reservoir (Fig. 3C), which varied along the sedimentary sequence. We prepared twelve samples for this core, but only three showed enough material for analysis, including the years 1962, 1966, and 1970. The predominant morphology was globular echinate (43%), followed by saddle (28%), bilobate (13%), and globular granulate (4%). The year 1962 showed the highest number of bilobate and saddle phytoliths and the lowest presence of globular echinate. In 1966 and 1970, there was a greater production of granular echinate and a lower number of bilobate and saddle types. The D/P index was 0.004 in 1962 and 1966 and 0.2 in 1970. The Iph (%) index was 70%, 73%, and 78% in 1962, 1966, and 1970, respectively. The Pa/P index was 0.58, 0.95, and 1.15 for 1962, 1966, and 1970, respectively.

Discussion

Phytolith from soils

In the present study, the saddle morphotype (Poaceae family) was very frequent in soil samples. Our data on soil phytolith assemblages differed from the analyzed phytolith of the most common plants in the Caatinga (other families), in which the prevalent types were tracheids, followed by polyhedral, globular granulate, elongate, and trichomes (lanceolate hair cells) [START_REF] Ricardo | Reference collection of plant phytoliths from the Caatinga biome, Northeast Brazil[END_REF]. A study on five soils sampled in Caatinga described globular granulate and globular echinate as the two more abundant morphotypes in the total assemblage [START_REF] Coe | Caracterização de fitólitos de plantas e assembleias modernas de solo da Caatinga como referência para reconstituições paleoambientais[END_REF]. [START_REF] Dias | Reconstituição paleoambiental de dunas vegetadas na Caatinga, em Aracati, Ceará, através de biomineralizações de sílica[END_REF] have found phytoliths well preserved, collected in dunes in the state of Ceará and they have inferred a relatively wetter climate, with a higher tree density than the current semi-arid, about 2200 years ago, for the Caatinga.

Although the saddle was the most frequent morphotype in the total assemblages, there was a different assemblage for each vegetation cover studied. The saddle was the main component in four (1, 4, 5, and 6) out of six samples from the Open Shrubby Caatinga.

This morphotype is linked to dry and warm conditions and occurs predominantly in the Chloridoideae subfamily (C4 grasses). However, the globular echinate is dominant in the assemblages of samples 2 and 3. Palms (Arecaceae family) and some bromeliads (Bromeliacea family) produce this morphotype. In the surroundings of sample 2, Attalea speciosa Mart. ex Spreng. (Arecaceae) is the predominant palm species of secondary ecological successions [START_REF] Souza | Os enclaves úmidos e sub-úmidos do semi-arido do nordeste brasileiro[END_REF], whereas Copernicia prunifera have been present in the riparian forest with "Carnauba" palm, surrounding the area of samples 2 and 3 (Fundação, 1994;[START_REF] Mesquita | Análise Socioambiental Do Médio Curso Do Rio Acaraú -CE[END_REF] (Fig. 1). Therefore, the riparian forest with Carnauba can contribute to the influx of globular echinate in the soil of the Open Shrubby Caatinga. In the Dense Shrubby Caatinga, the bilobate morphotype was the most frequent in sample 7, and it is mainly produced by the Panicoideae subfamily (C4 grasses), with records associated with warm and humid environments. The presence of papillae morphotypes is characteristic of the occurrence of Cyperaceae [START_REF] Piperno | Phytoliths: a Comprehensive Guide for Archaeologists and Palaeoecologists[END_REF]. Cyperaceae is a speciose, cosmopolitan family often found, but not exclusively, associated with poorly drained plant formations, such as swamps, lakes, riverbanks, and other water bodies [START_REF] Judd | Sistemática vegetal: um enfoque filogenético[END_REF]. Besides showing a high saddle percentage, sample 8 had similar values for bilobate, globular echinate, and globular granulate phytoliths, the second most representative types in this sample. In the state of Paraná, the gradual increase of the bilobate morphotype was observed during the transition from arid to more humid conditions [START_REF] Rasbold | Reconstrução paleoambiental de um depósito sedimentar por análises multiproxy, turvo, Estado do Paraná, Brasil[END_REF], and the rainfall rate in this site is 600.01-800 mm.

In the Caatinga Forest, only sample 12 showed a high number of globular granulate, a typical diagnostic morphotype for woody plants. We identified globular granulate producing plants in the surroundings of sample 12, but the production is different in each species, for example, it can be rare in Aspidorperma pyrifolium Mart. & Zucc. (Apocynaceae), occasional in Cordia oncocalyx Allemão (Boraginaceae), and frequent in Handroanthus impetiginosus (Mart. ex DC.) Mattos (Bignoniaceae), Anadenanthera colubrina (Vell.) Brenan (Fabaceae), and Mimosa tenuiflora [START_REF] Ricardo | Reference collection of plant phytoliths from the Caatinga biome, Northeast Brazil[END_REF]. This particular site is approximately 1 km away from the Aiuaba Ecological Station (6°35′55.48″S; 40° 7′16.77″O), a preserved area with relevant natural characteristics, created and protected by the Federal Government. Samples 9, 10, and 11 showed mainly saddle, followed by bilobate types. They were produced by Poaceae and could be related to degradation or anthropic pressure in more closed environments, when grasses, instead of woody plants, arise [START_REF] Ledru | When archives are missing, deciphering the effects of public policies and climate variability on the brazilian semi-arid region using sediment core studies[END_REF].

Therefore, the discrepant production of phytolith types among plants from the Caatinga can lead to inaccurate interpretations of assemblages by super or under-representing some of them. Studies on Caatinga plants [START_REF] Coe | Caracterização de fitólitos de plantas e assembleias modernas de solo da Caatinga como referência para reconstituições paleoambientais[END_REF][START_REF] Ricardo | Reference collection of plant phytoliths from the Caatinga biome, Northeast Brazil[END_REF] showed that in some families, such as Fabaceae, the phytolith production could vary according to the woody species; e.g., Libidibia ferrea (Mart. ex Tul.) L.P.Queiroz showed a very low production of globular granulate, Anadenanthera colubrina showed very high silicification of globular granulate, and Cenostigma bracteosum showed no silicification structures. So, in the Caatinga Forest, the type and abundance of woody species can affect the production of the globular granulate and its concentration in soil surfaces. Also the high degradation level in Caatinga Forest can affect the record of the phytolith type, because the replacement of the original vegetation by grasses can produce Grass Silica Short Cells.

Some elongate phytoliths, such as elongate psilate, have no taxonomic value, as it is produced in the epidermis of all grasses [START_REF] Bremond | Calibration des fonctions de transfert entre assemblages phytolithiques, structure des vegetations et variables bioclimatiques actuelles, pour l'integration de la ynamique des biomes herbaces dans les modeles de vegetation Ph[END_REF]. On the other hand, cuneiform bulliform phytoliths appear in the literature as water stress indicators [START_REF] Bremond | A phytolith index as a proxy of tree cover density in tropical areas: calibration with Leaf Area Index along a forest-savanna transect in southeastern Cameroon[END_REF].

Phytolith from the cores

The analyses in the two cores from the Open Shrubby Caatinga (Acaraú Mirim and Araras reservoirs) revealed a predominance of the globular echinate morphotype. This morphotype can be associated with the riparian forest in the edge of the reservoirs bearing a high number of Copernicia prunifera species (Arecaceae). These results were similar to phytolith assemblages from the cores in Lake Gbali, Central Africa [START_REF] Aleman | Influence of the local environment on lacustrine sedimentary phytolith records[END_REF]. The globular echinate morphotype has dominated (40%), though Arecaceae plants constituted only 5% of the riparian forest and were not part of the savanna matrix [START_REF] Aleman | Influence of the local environment on lacustrine sedimentary phytolith records[END_REF]. Palm trees produce abundant phytoliths and are expected to be over-represented in assemblages [START_REF] Pearsall | Paleoethnobotany: A Handbook of Procedures[END_REF][START_REF] Piperno | Phytoliths: a Comprehensive Guide for Archaeologists and Palaeoecologists[END_REF].

In the Acaraú Mirim core, we found a high frequency of globular echinate and a low frequency of globular granulate at the first record in 1965, which could be related to the open environment with the presence of the Arecaceae family and few woody plants. That scenario changed gradually until 1978. In 1983, the phytoliths revealed again evidence of an open environment, which could be related to the severe drought of 1982-83. According to pollen studies at the same core, the vegetation became more open between 1980 and 1988 [START_REF] Ledru | When archives are missing, deciphering the effects of public policies and climate variability on the brazilian semi-arid region using sediment core studies[END_REF]. The highest number of globular echinate occurred in 1991, whereas globular granulate peaked in 1996. A satellite image acquired in 1991 [START_REF] Ledru | When archives are missing, deciphering the effects of public policies and climate variability on the brazilian semi-arid region using sediment core studies[END_REF] and the significant rise in tree pollen taxa between 1990 and 1997 at the same studied core [START_REF] Ledru | When archives are missing, deciphering the effects of public policies and climate variability on the brazilian semi-arid region using sediment core studies[END_REF] indicated the return of a natural vegetation cover and a reduction in the extent of bare soil. That recovery was associated with the abandonment of cultivated land due to law changes (National Institute for Colonization and Agrarian Reform INCRA Law 97,844 July 19, 1989). Hence, the pollen analyses corroborated our phytolith record for the same reservoir.

In the Araras reservoir, the samples were considered sterile between 1993 and 1972 because they did not show the minimum GSSC number (200 morphotypes). Hence, the information about most layers was unreadable. The sediment size of samples may be favorable to the accumulation of phytoliths (finer particles, such as silt and clay), or may, when coarse (sandy), allow their percolation to deeper layers [START_REF] Coe | Caracterização de fitólitos de plantas e assembleias modernas de solo da Caatinga como referência para reconstituições paleoambientais[END_REF]. However, the main sediment size in this core is fine-grained, from top to bottom (0-60 cm) [START_REF] Ledru | When archives are missing, deciphering the effects of public policies and climate variability on the brazilian semi-arid region using sediment core studies[END_REF]. In the readable layers, from 1962 until 1970, there was an increasing frequency of globular granulate and globular echinate, and a decreasing frequency of GSSC. According to pollen studies in the Araras core [START_REF] Ledru | When archives are missing, deciphering the effects of public policies and climate variability on the brazilian semi-arid region using sediment core studies[END_REF], the main changes in the vegetation between 1963 and 1971 were reforestation with the expansion of Caatinga tree species, mainly Mimosa caesalpiniifolia Benth., M. tenuiflora and Astronium urundeuva, and an abrupt decrease in herb taxa.

So, pollen studies [START_REF] Ledru | When archives are missing, deciphering the effects of public policies and climate variability on the brazilian semi-arid region using sediment core studies[END_REF], in this same core, confirm the phytolith assemblage in Araras and Acaraú-Mirim reservoirs. In general, the pollen and phytolith data sets give a remarkably complementary picture of environmental change, despite the disparate range of taxa identifiable by each technique and the different taphonomy of each microfossil assemblage [START_REF] Kealhofer | A combined pollen and phytolith record for fourteen thousand years of vegetation change in northeastern Thailand[END_REF].

Final considerations

The present study analyzed the phytolith assemblage from a semi-arid region in Ceará, northeastern Brazil, in three different vegetation covers from the Crystalline Caatinga: Open Shrubby Caatinga, Dense Shrubby Caatinga, and Caatinga Forest. The phytolith morphotypes showed good conservation, so we could qualify and quantify them and the presence and variation of the Poaceae (subfamily Panicoideae and Chloridoideae), Arecaceae and Cyperaceae families. In soil samples, the saddle was the most common morphotype. It is related to warm and dry conditions and dominated most samples. However, we noticed exceptions. In the Open Shrubby Caatinga, there was a dominance of the globular echinate morphotype, produced by palms, which can be related to the riparian forest with Carnauba surrounding the collection area. In the Dense Shrubby Caatinga, the dominance of the bilobate morphotype could be associated with warm and humid environments. Finally, in the Caatinga Forest, only one sample, collected in an area with natural vegetation protected by Brazilian laws, showed the predominance of the globular granulate morphotype, which is related to woody plants. Therefore, there is evidence that the original flora is well recorded through phytolith assemblages in Ceará, corresponding to actual vegetation. In core samples from the two reservoirs in the Open Shrubby Caatinga, we could record the main changes in the vegetation structure, based on Grass Silica Short Cells and woody plant morphotype records. The variation in open or more closed vegetation cover could be related to natural phenomena (e.g., droughts) or anthropic pressure (e.g., abandonment of cultivated land due to a new law). The phytoliths present in the sediments of water reservoirs could be related to land use.

Fundings

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Figure1

  Figure1. Study area: A, distribution of the Caatinga in northeastern Brazil (Source: IBGE Instituto Brasileiro de Geografia e Estatística, 2004); B, Location of samples from the waterreservoir (Acaraú Mirim {2} and Araras {3} cores) and soil surfaces (1-12), in three different vegetation covers (Source: Figueiredo, 1997).

Figure 3 .

 3 Figure 3. Phytolith morphotypes found on the soil surface and sedimentary cores of reservoirs. (A) GSSC Morphotypes: a-k, bilobate; l, saddle; m, polylobate; n, trapeziform; o-p, rondel bilobate; q-r, cross. (B) a-g cuneiform bulliform. (C) Other morphotypes: a-b, tracheid; c, elongate crenate; d, elongate crenate; e, rodshaped? f, elongate dendriform; g, elongate crenate? h, elongate tuberculate; i, elongate sinuate; j, elongate sinuate; k, globular granulate; l, elongate lacunose; m-p, papillae; q, globular echinate; (D), diatom. The scale bar is 20 μm.

Figure 3 .

 3 Figure 3. Phytolith diagram. A, soil surface samples; B, Samples from the Acaraú Mirim core (Massapê City); C, Samples from the Araras core (Varjota City). GSSC: Grass Silica Short Cell. The dot indicates values for percentages ≤5%. Acronyms: CF, Caatinga Forest; DSC, Dense Shrubby Caatinga; and OSC, Open Shrubby Caatinga.

  This work was funded by: the Brazilian Coordination for the Improvement of Higher Education Personnel (CAPES), finance code 001 (Brazil), that granted the first author a post-doctoral scholarship (PNPD/CAPES Program) and financial support; "Sciences without boarder" (Portuguese acronym CSF) project "Evolution of Biodiversity loss of areas in desertification processes" (CNPq 400,890-2014-3/CNPq/CAPES/FAPs/#16/2014-3) and covered a visiting foreign professor fellowship for MPL at the Federal University of Ceara (UFC) in the frame of the CSF project (2014-2018); Ceará State Foundation for Scientific and Technological Development (FUNCAP); the Brazilian National Council for Scientific and Technological Development (CNPq) granted F.S. de Araújo a research productivity fellowship.

Acknowledgments

We thank the University of Montpellier, France; Sandrine Canal for the preparation of phytoliths; Olga Aquino for her help in sample preparation; Arlete Soares for providing us with the necessary equipment for phytolith analysis; Fernanda Kelly Gomes da Silva for her help with map data; Vivian Jeske-Pieruschka and Julia Sfair for helping us with data analysis.