- Email: [email protected]
Mixture resolution according to the percentage of robusta variety in order to detect adulteration in roasted coffee by near infrared spectroscopy
Analytica Chimica Acta 585 (2007) 266–276
Mixture resolution according to the percentage of robusta variety in order to detect adulteration in roasted coffee by near infrared spectroscopy C. Pizarro ∗ , I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz Department of Chemistry, University of La Rioja, C/Madre de Dios 51, 26006 Logro˜no (La Rioja), Spain Received 7 July 2006; received in revised form 19 October 2006; accepted 21 December 2006 Available online 17 January 2007
Abstract Near infrared spectroscopy (NIRS), combined with multivariate calibration methods, has been used to quantify the robusta variety content of roasted coffee samples, as a means for controlling and avoiding coffee adulteration, which is a very important issue taking into account the great variability of the final sale price depending on coffee varietal origin. In pursuit of this aim, PLS regression and a wavelet-based pre-processing method that we have recently developed called OWAVEC were applied, in order to simultaneously operate two crucial pre-processing steps in multivariate calibration: signal correction and data compression. Several pre-processing methods (mean centering, first derivative and two orthogonal signal correction methods, OSC and DOSC) were additionally applied in order to find calibration models with as best a predictive ability as possible and to evaluate the performance of the OWAVEC method, comparing the respective quality of the different regression models constructed. The calibration model developed after pre-processing derivative spectra by OWAVEC provided high quality results (0.79% RMSEP), the percentage of robusta variety being predicted with a reliability notably better than that associated with the models constructed from raw spectra and also from data corrected by other orthogonal signal correction methods, and showing a higher model simplicity. © 2007 Elsevier B.V. All rights reserved. Keywords: Near infrared spectroscopy; Multivariate calibration; Orthogonal wavelet correction; Wavelet transform; Genetic algorithms; Orthogonal signal correction; Coffee authenticity
1. Introduction Food authentication is a major challenge that has become increasingly important in recent years, due to the drive to guarantee the actual origin of a product and for determining whether it has been adulterated with contaminants or filled out with cheaper ingredients. Food industry, regulatory authorities and consumers are all interested in authenticating raw materials and food products in order to satisfy food quality and safety requirements [1–4]. In particular, assurance of the quality of roasted coffees has attracted widespread attention as a means for controlling and preventing coffee adulteration, and also given the great difference in the final sale price depending on a wide range of factors, including coffee varietal and geographic origin. The two varieties of
∗
Corresponding author. Tel.: +34 941299626; fax: +34 941299621. E-mail address: [email protected] (C. Pizarro).
0003-2670/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.aca.2006.12.057
economic significance in the global coffee trade are Coffea arabica L. (arabica coffee) and Coffea canephora Pierre (robusta coffee) [5]. Both species differ not only in relation to their botanical characteristics and physicochemical composition, but also in terms of commercial value, with arabica coffees commanding market prices 20–25% higher and being considered to be of better quality than robustas because of their superior taste and aroma. However, most commercially available coffees are actually blends of the two species, making it even more difficult to tackle coffee adulteration issues. Therefore, suitable methods are required in order to discriminate between coffee varieties and to detect potential adulterations of high quality coffee beans with poorer and cheaper types, thus ensuring authenticity, quality, safety and efficacy of final products to be commercialised. Extensive research has been carried out in recent years to develop reliable and specific coffee authentication methods [6–24]. Many of the recently developed approaches for determining coffee authenticity have focused mainly on coffee identification and classification on the basis of different types
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
of compositional data thanks to the application of different pattern recognition techniques [6–16]. Despite the relative success achieved by many of these approaches, it is important to consider that many analytical reference methods used to assess the chemical components to be later used as discriminant parameters may be quite expensive, elaborate and/or time-consuming. For this reason, simpler and faster methods, such as those based on spectroscopic techniques which can be easily implemented in routine analysis, have emerged as a very attractive and useful alternative tool for varietal identification purposes [17–24]. In particular, several studies have described the effective application of near infrared spectroscopy (NIRS) to address the problem of coffee authentication [20–24]. In this study we decided to adopt a less explored NIR-based quantitative approach to coffee authentication, namely a combination of NIR spectroscopy and chemometrics to optimize calibration procedures in order to both detect coffee adulteration in blends and the actual adulteration degree. Certain effects that occur inherently in diffuse reflectance NIR spectroscopy for solid samples, such as light scattering and influence of particle size, may be responsible for unexpected perturbations in spectra (baseline shifts, slope changes and non linearity) [25]. As a result, NIR signals contain not only chemical but also physical information about samples and measuring conditions, which may be irrelevant and may mask the significant information in the spectra and deteriorate the quality and predictive ability of the calibration models developed from raw spectra. Therefore, the application of suitable pre-processing methods, aimed at correcting spectral data by minimizing the contribution of physical effects to NIR signals and thus enhancing the relevant chemical information contained, is a crucial step in regression model development and improvement. Additionally, the huge amount of overlapping and redundant information present in NIR spectra often causes the useful information contained to be spread out over multiple points and at regions that are greatly affected by noise and interference. For this reason, NIR spectra are suitable for the application of efficient signal de-noising and compression methods, in order to remove both noise and not significant information from spectra. Consequently, the objective of this study was precisely to examine the feasibility of NIR spectroscopy, in combination with multivariate calibration techniques, to develop improved and reliable regression models which could later be used as a quick and accurate scanning tool for detecting coffee adulterations by quantifying the actual content of robusta variety in roasted coffee samples. We propose the application of a preprocessing method, recently developed by our research group, consisting of orthogonal wavelet correction (OWAVEC) coupled to the calibration step. The OWAVEC method was developed, and subsequently improved [26–28], to simultaneously accomplish two essential needs in multivariate calibration – signal correction and feature selection (data compression) – in an attempt to enhance the reliability and robustness of the finally constructed PLS regression model by combining the versatility that wavelet analysis has shown for signal processing [29–38] and the application of a particular orthogonal signal correction algorithm to remove from spectra information not related to the
267
response variable to be modelled and predicted. In particular, the last version of the OWAVEC method [28], which implemented a genetic algorithm as a feature selection method for extracting the significant corrected wavelet coefficients to be later used in PLS model development has been applied since it enabled a very substantial improvement in method performance. In order to better evaluate the performance of OWAVEC method applied to mixture resolution in roasted coffee samples, results provided by the final OWAVEC-PLS models constructed were compared with those obtained when working on both original data and spectra filtered by other orthogonal signal correction methods (OSC and DOSC). It should be noted that the purpose of this study was not to provide a definitive and immutable NIR calibration model for predicting robusta variety content in coffee samples, but to present a preliminary study to propose a strategy capable of accomplishing this task and for verifying its reliability and effectiveness in order to assess the genuineness of coffee samples. Although the data set used in this work was designed to cover, as best as possible, the great natural variability inherent in commercially available coffees by considering different roasting conditions and degrees, and varied geographical origins, it is quite clear that the dynamic nature of the coffee market and the particular needs and production lines of a given coffee company would require a more exhaustive or specific collection of a suitable set of calibration samples to serve as the basis for developing, using the methodology proposed here, a calibration model for use in the routine quality control of coffee samples. 2. Materials and methods 2.1. OWAVEC method This paper does not aim to describe in detail the fundamentals of all pre-processing methods that were tested to remove NIR spectral variation not related to the response variable and to enhance spectral selectivity, moreover taking into account that some of these methods (first derivative, OSC and DOSC) are quite well known and usually applied in multivariate calibration. Nevertheless, considering the recent introduction of OWAVEC method for signal pre-processing, a brief discussion about its main application features will be provided in order to facilitate a clear understanding of the procedure carried out by OWAVEC. Fig. 1 shows a schematic view of the four main sequential steps into which the OWAVEC algorithm can be broken down. As a first step, OWAVEC applied the discrete wavelet transform (DWT) to decompose each spectrum in the wavelet domain, using a selected wavelet function in this transformation. In DWT, the choice of different wavelet basis functions gives rise to different decompositions of a signal in the resulting wavelet domain. In the present study, and taking into account our experience in wavelet analysis of NIR signals, the biorthogonal Daubechies-4 wavelet was selected, partly because its two vanishing moments provide high compressibility for smooth signals such as NIR data. The DWT requires the length of the analysed signal to be dyadic (2n , where n is an integer). However, if this constraint is not fulfilled several methods can be applied
268
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
Fig. 1. Schematic view of all the steps involved in the OWAVEC method. (1) Centered spectra are padded to the nearest length 2n . Each spectrum is transformed using DWT and Daubechies-4 as wavelet basis function. (2) The wavelet coefficients matrix is deflated of information unrelated to the considered response (orthogonalization). (3) Wavelet coefficients matrix compression by using GAs for subset selection of significant wavelet coefficients. (4) Corrected signals can be reconstructed into the wavelength domain using IWT. The final PLS regression model was developed from the resulting corrected and compressed wavelet coefficients matrix.
in order to extrapolate the signal. In the OWAVEC wavelet decomposition step, all signal extrapolations were carried out by linear padding in order to minimise edge effects. Another key parameter to be optimised when applying DWT for regression purposes is the optimal wavelet decomposition level. Thus, in order to select the most suitable wavelet decomposition to be performed, several decomposition levels were tested, always respecting the maximum allowed wavelet decomposition level associated with the extended length of the analysed signals and the use of Daubechies-4 as wavelet function (in this case, the maximum decomposition level permitted was six). The final choice of the optimal decomposition level was based on several criteria, trying to both maximize the compression rate achieved and the predictive ability of the calibration model developed once the OWAVEC method had been fully executed. Once the wavelet coefficients matrix corresponding to the calibration set had been computed, an orthogonalization algorithm was applied to deflate it of information not related to the response variable (i.e. the robusta content of coffee samples). Further details on the orthogonalization procedure implemented in OWAVEC can be found in reference [26]. Once the orthogonalization algorithm had been applied on the calibration wavelet coefficients matrix, and the spectra in the prediction set had been properly decomposed by wavelet analysis, the computed matrix of weights of the original wavelet coefficients in the orthogonal subspace of these coefficients (W) could be used directly to obtain the corrected wavelet coefficients matrix for test (new) data. A critical parameter to be optimized when performing wavelet coefficients matrix orthogonalization is the value of the tolerance factor used to calculate the generalized inverse of the wavelet coefficients matrix. The use of a generalized inverse to
obtain the matrix of weights W implies the loss of the complete orthogonality constraint as an attempt to only take into account stable directions in wavelet coefficients matrix for computing its orthogonal subspace, which means that some correlation between the information removed from the wavelet coefficients matrix and the response variable is allowed. The adopted degree of deviation from orthogonality will depend precisely on the used value for the tolerance factor, so that for each data set studied the tolerance factor should be properly tuned. Thus, different tolerance factor values were tested in each case, so that maximizing the predictive power of the regression model subsequently developed from the resulting corrected and compressed wavelet matrix, while not showing any noticeable signs of overfitting, was finally selected as the optimal value. Finally, as far as the OWAVEC data compression (feature selection) step was concerned, before searching for significant subsets of corrected wavelet coefficients by GAs, the settings required in the genetic algorithm run need to be selected. The effect of varying the values of some GAs design parameters (number of chromosomes; number of elitist chromosomes; percentage of mutation; mean, maximum and minimum numbers of coefficients to be retained in the same chromosome; maximum number of cycles without fitness function (RMSECV) improvement) was evaluated to favour algorithm convergence and to search for an improved solution. The final choice of the optimal wavelet decomposition parameters, tolerance factor and GAs settings was based on cross-validation results, also seeking to reach as high a compression rate as possible without losing model performance. Moreover, for interpretative and comparative purposes, the spectra corrected after pre-processing by OWAVEC were recon-
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
structed into the original wavelength domain from the respective corrected and compressed wavelet coefficients matrices by inverse wavelet transform (IDWT), for both calibration and prediction sets. This can be very useful for identifying the spectral regions responsible for a successful prediction. Since the only pre-processing procedure performed during OWAVEC application was mean-centering, the wavelet coefficients matrix retained the original variance of input spectra, so after the optional reconstruction step of the filtered signals the global shape of the starting spectra would be maintained. 2.2. Samples Eighty three pure variety roasted coffee samples from varied origins (36 arabica and 47 robusta coffees), which were processed under different roasting conditions, were used for this study, in such a way that they were representative of the variability present in coffee market. In addition, a number of varietal blends were prepared in the laboratory from these pure variety coffees. Taking into account that the preparation of all possible blends resulting from mixing the available roasted coffees from different varieties would have been impractical, we decided to select the samples most representative of both arabica and robusta in order to later combine them to obtain suitable blends covering the overall variability of the whole set of blends potentially derived from the pure coffees considered. The application of PCA on mean centered NIR spectra provided an effective and easy-to-implement tool for selecting representative samples from arabica and robusta varieties. Fig. 2 shows a bidimensional representation of PC1 and PC2 scores accounting for 88.33% of the variance in the roasted coffee NIR spectral data, labelled according to their coffee variety: (1) arabica coffees; (2) robusta coffees. Two sample groups appeared slightly separated by the first bisectrix of the two component axes, suggesting the presence of two different clusters just associated with the two varieties considered. Thus, the centroid of each class was selected as representative of each variety and the two extreme
269
samples within each cluster were selected in order to take into account the variability associated to each group. Next, the three samples selected as representative for each variety were combined to make suitable coffee blends with robusta content in the final blends ranging from 0 to 60% (w/w), in such a way that a total number of 108 blends were consequently generated. This specific range of blending (based on the gradual addition of increasing robusta amounts in final blends up to 60%) was particularly studied in order to cover the actual range of adulteration encountered, which is more specifically located within relatively low and still financially profitable adulteration levels. Fraudulent practice in the elaboration of coffee blends from arabica and robusta pure coffee varieties is largely due to economic reasons, i.e. to reduce costs due to the lower price of robusta coffee as compared to arabica coffee, labelling mixed variety coffees as pure arabica or commercialising blends which actually contain a lower percentage of arabica than that declared. Therefore, the resulting data set used in the present study comprised a total of 191 roasted coffee samples (83 pure varieties and 108 blends). In order to develop and validate the regression models, the entire data set was split into two independent subsets: a calibration set with 133 coffees (comprising 70% of samples and covering the whole variability of the original data set to ensure its representativity); and a test set with the remaining 58 samples. Special care was taken to ensure suitable composition of the external test set, checking that samples of both pure varieties and all different compositional blends were included. 2.3. Apparatus and software NIR spectra were recorded on a NIRSystems 5000 near infrared spectrophotometer (FOSS, NIRSystems) equipped with a reflectance detector and a sample transport module. The instrument was controlled by a compatible PC and Vision 2.22 (FOSS, NIRSystems) was used to acquire the data. The OWAVEC method was entirely programmed, including the data compression step by GAs, using MATLAB® technical computing language, version 6.5 [39]. The steps of the OWAVEC procedure involving wavelet analysis were implemented using the Wavelet Toolbox version 2 (for use with MATLAB® ) [40]. The orthogonal signal correction (OSC) and direct orthogonal signal correction (DOSC) pre-processing methods were applied simply to serve as a comparative reference for evaluating the results provided by OWAVEC. Thus, the OSC method developed by Wise and Gallagher [41] (available in PLS-Toolbox 2.1 for use with MATLAB® [42]), and the MATLAB® code of the DOSC method developed by Westerhuis et al. [43] were used for calculations. For the construction of all the calibration models, PLS regression was performed using the PLS program contained in the V-Parvus package [44]. 2.4. Recording of NIR spectra
Fig. 2. Scores of the 83 roasted coffee samples from arabica (labelled as 1) and robusta (labelled as 2) pure varieties on the first two principal components explaining the variability in the NIR spectral data.
Reflectance spectra were obtained directly from untreated samples. Due care was taken to ensure that the same amount of
270
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
sample was always used to fill up the sample cell. Each spectrum was obtained from 32 scans performed at 2 nm intervals within the wavelength range 1100–2500 nm, with five replicates for each individual sample. The samples were decompacted between recordings. An average spectrum was subsequently obtained from the replicates collected for each coffee sample. 2.5. Data analysis Several data pre-processing methods (mean-centering, first derivative, orthogonal signal correction, OWAVEC or a combination of several of these options) were applied before calibration development in order to find models with as high a predictive power as possible and to compare OWAVEC method performance with that achieved using other pre-treatments. Fig. 3 shows a schematic diagram of the different steps executed for the different pre-processing methods tested. In this way, various PLS regression models were constructed between calibration NIR spectra pre-processed by the different pretreatments used and the robusta variety content of roasted coffee samples (when the OWAVEC method was applied, the corrected and compressed wavelet coefficients matrix was considered for calibration development instead of using NIR spectra in the original wavelength domain). All PLS models were developed from
the column-wise centered data matrix, regardless of whether or not an additional pre-processing method was also applied on spectral variables. The first derivative of spectra was applied in some cases in order to improve spectral differences and to correct baseline effects. In order to avoid noise enhancement as a consequence of derivative, data were smoothed by the Savitzky–Golay method, using a cubic degree polynomial with a window size of 13 points. The optimal complexity of every calibration model was assessed by cross-validation (all the PLS models were built by cross-validation using five deletion groups). The actual predictive ability of all the constructed regression models was also validated by testing their performance on an external test set, in order to control and avoid possible overfitting (particularly when applying orthogonal signal correction methods due to the high risk of overfitting inherent in them). When OSC was applied as a pre-processing method, before the mean-centering step, spectra were transformed into their first derivative spectra, since this additional step substantially improved the quality of the resulting models. When the orthogonal signal correction methods were applied, in order to search for models with high predictive ability and to avoid over-fitted solutions, the suitable orthogonal correction degree to be applied was determined as follows: the number of orthogonal-LVs (OSC) or orthogonal-PCs (DOSC)
Fig. 3. Flow diagram summarising the different pre-processing methods (or combination of methods) applied to the data set.
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
to be removed from spectra was varied from 1 to 6, in such a way that the optimal number of orthogonal components to be subtracted was chosen according to the results obtained in the cross-validation procedure. The root mean square error (RMSE) of the residuals obtained (termed RMSEC in calibration, RMSECV in cross-validation and RMSEP in external prediction) was used to evaluate and compare the respective performance of the different PLS models constructed, since this parameter represented an objective measurement of the resulting predictive ability of the model, which was dimensionally comparable to the studied response. In this particular calibration study, the RMSE values were already obtained in percentages, since the range of robusta variety content in coffee samples varied precisely from 0 to 100%. Likewise, the correlation plots between measured and computed/predicted response values, as well as their corresponding coefficients of determination (R2 ), were also obtained to serve as additional tools to evaluate models accuracy. Nevertheless, the prediction ability of a regression model must be accompanied by an estimate of its uncertainty. Thus, prediction intervals at 95% confidence were computed for each regression model (in both calibration and external prediction) in order to better estimate and compare their respective reliability. 3. Results and discussion 3.1. Regression models for quantifying robusta variety content from roasted coffee NIR spectra Table 1 shows the calibration and prediction errors prompted by the PLS models developed from mean-centered and first derivative NIR spectra to quantify the robusta content of roasted coffee samples. In spite of the high predictive ability achieved when working on mean-centered data (2.98% RMSEP), a large number of latent variables (10 LVs) had to be used to develop the model in order to compensate for the structured noise (baseline and scatter effects) present in the original NIR spectra. A high model complexity could be seen as a limitation for the use of the corresponding calibration model as an accurate tool for detecting coffee adulteration, since it would be desirable to use a small number of components providing a minimal error for a simpler, and probably more robust, calibration model. The transformation of raw signals into their first derivative spectra, prior to calibration model construction, provided a significant decrease in prediction error with respect to the model developed from original spectra (2.04% RMSEP), but at the expense of maintaining high model complexity (8 LVs), which suggested that
271
the contribution to spectra of variation not due to the parameter of interest was still significant. The results corresponding to the calibration models developed from NIR spectra corrected by both OSC and DOSC pre-processing methods (after selection of the optimal number of orthogonal factors to be removed in each case) are also included in Table 1. Both OSC and DOSC methods were applied on NIR spectra taking into account, in order to correct signals, the percentage of robusta variety contained in each sample as response variable. The corrected resulting spectra were then used to develop the respective PLS models. Considering all the models based on signals corrected by both methods, it was observed that the calibration errors and the corresponding model complexities decreased once each orthogonal component was removed from the spectra. Nevertheless, this apparent gradual improvement in model performance was limited for both correction methods, since in the case of OSC the removal of more than two orthogonal latent variables reduced the calibration error at the expense of increasing the prediction error (thus giving an extremely overfitted solution), whereas when applying DOSC a correction beyond five orthogonal PCs did not provide any substantial improvement. The calibration models finally selected as optimal after applying both methods provided similar results with regard to predictive ability (1.36–1.33% RMSEP for OSC and DOSC methods, respectively), and represented a considerable improvement compared to the original model in terms of both model reliability and complexity. However, it should be noted that, despite carefully selecting the number of orthogonalLVs to be removed, OSC led to an overfitted solution since the prediction error almost doubled the corresponding calibration error. Despite the quality improvement prompted by the models developed by applying OSC and DOSC methods, the aim of this study was to search for a suitable pre-processing method capable of both performing efficient signal correction and effective data compression, in such a way that it would be possible to obtain an improved final regression model that would be less complex, more robust and parsimonious and with an increased predictive ability, and which could be applied for the accurate detection of potential coffee adulterations. The major advantage of the OWAVEC method is its ability to compress in few variables (the selected subset of corrected wavelet coefficients) the information responsible for a successful prediction. This will likely not only improve model performance with regard to prediction accuracy, but also make the resulting model more parsimonious, i.e. reduce the dimensionality of the predictors (coefficients) matrix prior to calibration, thus providing a more robust model.
Table 1 Calibration (RMSEC) and prediction (RMSEP) errors and percentages of explained variance provided by the original PLS model (constructed from mean-centered data), the PLS model developed from first derivative spectra and the PLS models built after applying OSC and DOSC methods to NIR spectra
Centering First derivative First derivative + OSC (2 O-LVs removed) DOSC (5 O-PCs removed)
PLS-LVs
%RMSEC
%EXPL.VAR. (cal)
%RMSEP
%EXPL.VAR. (test)
10 8 3 1
2.13 1.67 0.77 1.03
99.67 99.80 99.96 99.92
2.98 2.04 1.36 1.33
99.38 99.71 99.87 99.88
Only the model yielding the best performance for each method is shown.
272
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
Table 2 Calibration (RMSEC) and prediction (RMSEP) errors and percentages of explained variance for PLS models developed from the corrected and compressed coefficients matrix obtained after the application of OWAVEC, when working on both original and first derivative spectra LVs
%RMSEC
%EXPL.VAR. (cal)
%RMSEP
%EXPL.VAR. (test)
LEVEL
TOL
MED/MAX/MIN
CROM/ELIT/CP
NCOF
%ZERO COEFF.
OWAVEC (centered data) 1 1.18 99.90 1 1.17 99.90 1 1.18 99.90 1 1.18 99.90 1 1.18 99.90
1.36 1.34 1.33 1.37 1.36
99.87 99.87 99.88 99.87 99.87
2 2 2 2 2
0.0015 0.0015 0.0015 0.0015 0.0015
8/10/2 4/10/2 4/6/2 5/8/2 5/7/3
20/2/5 20/2/5 20/2/5 20/2/5 20/2/5
4 4 3 4 5
99.44 99.44 99.58 99.44 99.30
First derivative + OWAVEC 2 0.73 99.96 2 0.66 99.97 2 0.68 99.97 2 0.67 99.97 2 0.64 99.97
0.83 0.81 0.79 0.84 0.81
99.95 99.95 99.96 99.95 99.95
2 2 2 2 2
0.0005 0.0005 0.0005 0.0005 0.0005
7/10/2 7/10/5 8/10/5 8/10/3 8/10/4
20/2/5 20/2/5 20/2/5 25/1/5 20/2/5
9 10 10 10 10
98.74 98.60 98.60 98.60 98.60
For each pre-treatment combination tested, the results obtained after performing several runs using different choices of GA parameters are shown. The particular OWAVEC-GA parameters used in each run (including optimal wavelet decomposition level, tolerance factor to compute the generalized inverse, number of chromosomes, number of elitist chromosomes, percentage of mutation, mean, maximum and minimum numbers of coefficients allowed for selection in the same chromosome and maximum number of cycles) are also included. Resulting compression scores (number of retained coefficients and percentage of zero coefficients) are shown.
Fig. 4. Plot of NIR computed (in cross validation) vs. reference robusta contents for calibration samples. Results from the PLS model based on: (a) original data (only mean centered); (b) first derivative spectra; (c) OSC corrected first derivative spectra; (d) DOSC filtered spectra; (e) OWAVEC corrected and compressed wavelet coefficients matrices computed from original data; (d) OWAVEC corrected and compressed wavelet coefficients matrices working on first derivative spectra. Only a model representative of all constructed models in the different GA runs performed when applying OWAVEC method was shown. The 95% confidence interval of the linear fitting was shown (grey dashed line).
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
Due to the stochastic character of the optimisation procedure carried out by OWAVEC-GA, a preliminary explanatory study, using different choices of GA parameters settings, was performed to find an optimal search region, thus to obtain a number of selected coefficients large enough to achieve a good fit but small enough to avoid overfitting. The agreement observed between results obtained after different GA runs for each data pre-treatment tested indicated a good convergence of the genetic algorithm applied. Table 2 summarizes the results obtained with the different calibrations models developed from the respective subsets of OWAVEC corrected and compressed wavelet coefficients selected by genetic algorithms in five different runs of the method performed when working on both original and first derivative spectra. All very similar PLS models developed after the direct application of the OWAVEC-GA method to original data provided high predictive ability (≈1.35% RMSEP) with the minimum possible complexity and no appreciable signs of overfitting (≈1.18% RMSEC), which represented a substan-
273
tial improvement with respect to the models constructed from spectra not processed by orthogonal signal correction methods. Moreover, the impressive compression rates resulting from the feature selection step carried out by GAs (99.30–99.58% zero coefficients), served to demonstrated the enhanced effectiveness of the OWAVEC-GA method with respect to the other orthogonal signal correction methods tested here, since although no substantial improvement in predictive ability was achieved, the final model improved substantially in terms of parsimony (only 3–5 from a total of 713 corrected wavelet coefficients were selected in the different runs) and overfitting was reduced compared to OSC based model. Finally, we tried to further improve the reliability and stability of OWAVEC-based PLS models by taking advantage of the proven ability of first derivative to reduce baseline offset variations in NIR signals and to increase spectral resolution. Thus, Table 2 also includes the results corresponding to the regression models constructed from the respective corrected and compressed wavelet coefficients matrices resulting from
Fig. 5. Plot of NIR predicted (in external validation) vs. reference robusta contents for test samples. Results from the PLS model based on: (a) original data (only mean centered); (b) first derivative spectra; (c) OSC corrected first derivative spectra; (d) DOSC filtered spectra; (e) OWAVEC corrected and compressed wavelet coefficients matrices computed from original data; (d) OWAVEC corrected and compressed wavelet coefficients matrices working on first derivative spectra. Only a model representative of all constructed models in the different GA runs performed when applying OWAVEC method was shown. The 95% confidence interval of the linear fitting was shown (grey dashed line).
274
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
the application of OWAVEC on mean-centered first derivative NIR spectra, with the performance of several GA-searches. All PLS models developed after the coupled application of the first derivative and OWAVEC methods yielded substantially improved calibration and prediction results for assessing the robusta content of roasted coffees not only with respect to the models constructed from original data, but also with regard to the models developed from data processed by both other signal correction methods and to the OWAVEC method when working directly on mean-centered data. In addition to the extremely low prediction errors achieved (0.79–0.84% RMSEP) with low model complexities (only 2 LVs) and lack of overfitting (0.64–0.73% RMSEC), special consideration must be given to the high compression rates obtained thanks to the implementation of GAs as a feature selection method in OWAVEC (only 9–10 from a total of 713 wavelet coefficients were retained; 98.60–98.74% zero coefficients), which allowed us to obtain a more parsimonious calibration model. We could also observe that the value of the tolerance factor used in the OWAVEC orthogonalization step when working on first derivative spectra was lower than for raw spectra. This fact made clear that the
treatment of NIR spectra by first derivative prior to OWAVEC application served to effectively minimize background effects, since the previous correction for unwanted systematic variation achieved by derivative resulted, to a certain a extent, in a reduction of non-stable directions in the corresponding wavelet matrix, thus enabling us to consider a lower correlation between the discarded information and the response variable. In order to gain further insight into the accuracy of the different calibration models developed and to compare their relative goodness, linear regression of the computed/predicted values for the robusta content versus reference values, for both calibration and test sets, was applied (Figs. 4 and 5). The grey dashed lines are the prediction bands at 95% confidence limit, in such a way that the smaller was the prediction interval, the model precise was the model. Thus, of all spectral pre-processing methods which were tested to optimize the calibration model for detecting potential coffee adulterations, the coupled application of first derivative and OWAVEC not only gave rise to the model displaying the lowest complexity, RMSEC and RMSEP, but the highest correlation for both the calibration and prediction sets (R2 = 1) together with the lowest uncertainty in prediction (the most nar-
Fig. 6. Spectral profiles of both calibration and test roasted coffee samples corresponding to (a) Original NIR spectra. (b) First derivative spectra. (c) OSC corrected first derivative spectra. (d) DOSC corrected spectra. (e) OWAVEC corrected reconstructed spectra. (f) OWAVEC corrected reconstructed first derivative spectra.
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
row prediction intervals). Moreover, to provide an even more precise idea of the level of accuracy achieved with the proposed model, the 95% confidence interval for the prediction error might be estimated as ±1.98 RMSEP, i.e. ±1.60% adulteration. This high accuracy achieved in the resolution of arabica/robusta mixtures was even more remarkable considering that the detection limit of NIR spectroscopy is about 0.1–1%. 3.2. Observations on NIR spectra The scatter effects occurring inherently in NIRS can be significant and produce a considerable expansion of the absorbance interval for the individual wavelengths, which can often hinder the extraction of significant information from the raw data. This phenomenon was clearly revealed in the spectral profiles corresponding to original NIR signals shown in Fig. 6a. For this reason, NIR spectra need to be pre-processed in order to remove systematic noise. First derivative pre-treatment did not remove all the spectral variations caused by the scatter effects, or rather, it only reduced them at some spectral regions (Fig. 6b). One of the main objectives that orthogonal signal correction methods attempted to accomplish was the minimization of systematic effects, by removing information not related to the response variable, in order to improve the quality of the resulting calibration models. In fact, as can be observed from the spectra obtained after data filtering by the OSC (Fig. 6c) and DOSC (Fig. 6d) methods, the spectral differences were significantly minimized with respect to the raw signals, proving the relative efficiency of both methods to correct data, trying to preserve only spectral information closely related to the robusta content of coffee samples. In the previous sub-section, we demonstrated the usefulness and effectiveness of the OWAVEC-GA method for correcting and compressing data into a wavelet domain prior to the development of a reliable and robust calibration model, as well as its better performance potential with respect to the other preprocessing methods tested, mainly by expounding numerical results. However, a graphical analysis of the spectral profiles reconstructed from the respective subsets of wavelet coefficients finally selected by OWAVEC can also help to demonstrate the goodness of the calibration strategy proposed here. These reconstructed profiles can be used to extract significant information about which specific spectral regions, representing the underlying features responsible for a successful prediction of the robusta content of roasted coffees, were associated with the corrected coefficients selected for the development of the final regression model. As can be seen, from plots displaying both original (Fig. 6e) and first derivative (Fig. 6f) reconstructed spectra after pre-processing by OWAVEC, the spectral differences between samples were not only minimised but restricted to a very limited number of regions since exclusively information related to robusta content was selected to be kept. 4. Conclusions The results reported in this study show that a method based on the combination of NIR spectroscopy and multivariate cal-
275
ibration can be used to detect and quantify potential coffee adulterations. However, the success of this methodology will depend largely on signal pre-processing methods applied to minimize the spectral variation not due to the parameter of interest but due to variation in experimental or sample conditions. The regression models developed for predicting the robusta variety content of coffee samples, after pre-processing first derivative NIR spectra by OWAVEC, were considerably better compared to the models obtained from raw data (in terms of both reliability and parsimony) and improved results were also obtained with respect to the other orthogonal signal correction methods tested. The OWAVEC-GA selection of a small subset of corrected wavelet coefficients responsible for accurate modelling not only led to an improvement of predictive ability but also to a decrease of the complexity associated with the resulting model, making it more parsimonious and probably more robust, since spectral regions exhibiting a large variance in raw spectra, as well as noise, not related to the robusta content of roasted coffee samples were discarded. Although the work reported here is only a feasibility study and further research should be carried out before its value for authentication purposes can be truly estimated, the high quality calibration models constructed, thanks to the high effectiveness of OWAVEC as a correction and compression tool, suggest that the proposed regression methodology would be a suitable tool for the quantitative detection of robusta contamination of arabica coffees in the authenticity assessment of coffee. Acknowledgements The authors thank the Ministry of Science and Technology (Project No. 2FD1997-0491) for its financial support, as well as Professor Michele Forina for providing us with the last version of Parvus package. References [1] P.R. Ashurst, M.J. Dennis (Eds.), Food Authentication, Blackie Academic & Professional, London, UK, 1996. [2] R.S. Singhal, P.R. Kulkarni, D.V. Rege, Handbook of Indices of Food Quality and Authenticity, Woodhead Publishing Ltd., Cambridge, UK, 1997. [3] P.R. Ashurst, M.J. Dennis (Eds.), Analytical Methods of Food Authentication, Blackie Academic & Professional, London, UK, 1998. [4] M. Lees (Ed.), Food Authenticity and Traceability, Woodhead Publishing Ltd., Cambridge, UK, 2003. [5] R. Viani, A. Illy (Eds.), Espresso Coffee: The Science of Quality, second ed., Academic Press, London, UK, 2004. [6] N. Frega, F. Bocci, G. Lercker, Industrie Alimentari 34 (1995) 705. [7] C.P. Bicchi, M.P. Ombretta, G. Pellegrino, A.C. Vanni, J. Agric. Food Chem. 45 (1997) 4680. [8] M.J. Mart´ın, F. Pablos, A.G. Gonz´alez, Talanta 46 (1998) 1259. [9] F. Carrera, M. Le´on-Camacho, F. Pablos, A.G. Gonz´alez, Anal. Chim. Acta 370 (1998) 131. [10] F. Pablos, A.G. Gonz´alez, M.J. Mart´ın, M.S. Valdenebro, M. Le´onCamacho, Analyst 124 (1999) 999. [11] M.J. Mart´ın, F. Pablos, A.G. Gonz´alez, Food Chem. 66 (1999) 365. [12] M.J. Mart´ın, F. Pablos, A.G. Gonz´alez, M.S. Valdenebro, M. LeonCamacho, Talanta 54 (2001) 291. [13] A.G. Gonz´alez, F. Pablos, M.J. Mart´ın, M. Le´on-Camacho, M.S. Valdenebro, Food Chem. 73 (2001) 93.
276
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
[14] S. Casal, M.R. Alves, E. Mendes, M.B.P.P. Oliveira, M.A. Ferreira, J. Agric. Food Chem. 51 (2003) 6495. [15] C.G. Zambonin, L. Balest, G.E. De Benedetto, F. Palmisano, Talanta 66 (2005) 261. [16] L. Mondello, R. Costa, P.Q. Tranchida, P. Dugo, M. Lo Presti, S. Festa, A. Fazio, G. Dugo, J. Sep. Sci. 28 (2005) 1101. [17] E.K. Kemsley, S. Ruault, R.H. Wilson, Food Chem. 54 (1995) 321. [18] R. Briandet, E.K. Kemsley, R.H. Wilson, J. Agric. Food Chem. 44 (1996) 170. [19] R. Briandet, E.K. Kemsley, R.H. Wilson, J. Sci. Food Agric. 71 (1996) 359. [20] G. Downey, J. Boussion, D. Beauchˆene, J. Near Infrared Spectrosc. 2 (1994) 85. [21] G. Downey, J. Boussion, J. Sci. Food Agric. 71 (1996) 41. [22] G. Downey, B. Spengler, J. Agric. Food Res. 35 (1996) 179. [23] G. Downey, R. Briandet, R.H. Wilson, E.K. Kemsley, J. Agric. Food Chem. 45 (1997) 4357. [24] I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz, C. Pizarro, Anal. Chim. Acta 514 (2004) 57. [25] B.G. Osborne, T. Fearn, P.H. Hindle, Practical NIR Spectroscopy with Applications in Food and Beverage Analysis, second ed., Longman Scientific and Technical, Harlow, UK, 1993. [26] I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz, C. Pizarro, Anal. Chim. Acta 515 (2004) 31. [27] I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz, C. Pizarro, Anal. Chim. Acta 544 (2005) 89. [28] I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz, D. G´omez-C´amara, C. Pizarro, Anal. Chim. Acta 555 (2006) 84. [29] R.R. Coifman, Y. Meyer, M.V. Wickerhauser, Wavelet analysis and signal processing, in: M.B. Ruskai, G. Beylkin, R. Coifman, I. Daubechies, S.
[30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]
Mallat, Y. Meyer, L. Raphael (Eds.), Wavelets and Their Applications, Jones & Bartlett, Boston, USA, 1992, pp. 153–178. B. Walczak, D.L. Massart, Trends Anal. Chem. 16 (1997) 451. B.K. Alsberg, D.B. Kell, J.J. Rowland, M.K. Winson, A.M. Woodward, Analyst 122 (1997) 645. B. Walczak, D.L. Massart, Chemometr. Intell. Lab. Syst. 36 (1997) 81. V.J. Barclay, R.F. Bonner, I.P. Hamilton, Anal. Chem. 69 (1997) 78. A.K.M. Leung, F.T. Chau, J.B. Gao, Chemometr. Intell. Lab. Syst. 43 (1998) 165. J. Trygg, S. Wold, Chemometr. Intell. Lab. Syst. 42 (1998) 209. S. Mallat, A Wavelet Tour of Signal Processing, second ed., Academic Press, London, UK, 1999. K. Jetter, U. Depczynski, K. Molt, A. Niem¨oller, Anal. Chim. Acta 420 (2000) 169. B. Walczak (Ed.), Wavelets in Chemistry, Elsevier Science B.V., The Netherlands, 2000. MATLAB® 6.5, The MathWorks, Natick, USA, 2002. M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi, Wavelet Toolbox for Use with MATLAB® , The MathWorks, Natick, USA, 2000. B.M. Wise, N.B. Gallagher, http://www.eigenvector.com/MATLAB/ OSC.html. B.M. Wise, N.B. Gallagher, PLS Toolbox 2.1, Eigenvector Research Inc., USA, 1998. J.A. Westerhuis, S. de Jong, A.K. Smilde, Chemometr. Intell. Lab. Syst. 56 (2001) 13. M. Forina, S. Lanteri, C. Armanino, C. Cerrato-Oliveros, C. Casolino, VPARVUS 2004, an extendable package of programs for explorative data analysis, classification and regression analysis, Dip. Chimica e Tecnologie Farmaceutiche ed Alimentari, University of Genova.
Mixture resolution according to the percentage of robusta variety in order to detect adulteration in roasted coffee by near infrared spectroscopy C. Pizarro ∗ , I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz Department of Chemistry, University of La Rioja, C/Madre de Dios 51, 26006 Logro˜no (La Rioja), Spain Received 7 July 2006; received in revised form 19 October 2006; accepted 21 December 2006 Available online 17 January 2007
Abstract Near infrared spectroscopy (NIRS), combined with multivariate calibration methods, has been used to quantify the robusta variety content of roasted coffee samples, as a means for controlling and avoiding coffee adulteration, which is a very important issue taking into account the great variability of the final sale price depending on coffee varietal origin. In pursuit of this aim, PLS regression and a wavelet-based pre-processing method that we have recently developed called OWAVEC were applied, in order to simultaneously operate two crucial pre-processing steps in multivariate calibration: signal correction and data compression. Several pre-processing methods (mean centering, first derivative and two orthogonal signal correction methods, OSC and DOSC) were additionally applied in order to find calibration models with as best a predictive ability as possible and to evaluate the performance of the OWAVEC method, comparing the respective quality of the different regression models constructed. The calibration model developed after pre-processing derivative spectra by OWAVEC provided high quality results (0.79% RMSEP), the percentage of robusta variety being predicted with a reliability notably better than that associated with the models constructed from raw spectra and also from data corrected by other orthogonal signal correction methods, and showing a higher model simplicity. © 2007 Elsevier B.V. All rights reserved. Keywords: Near infrared spectroscopy; Multivariate calibration; Orthogonal wavelet correction; Wavelet transform; Genetic algorithms; Orthogonal signal correction; Coffee authenticity
1. Introduction Food authentication is a major challenge that has become increasingly important in recent years, due to the drive to guarantee the actual origin of a product and for determining whether it has been adulterated with contaminants or filled out with cheaper ingredients. Food industry, regulatory authorities and consumers are all interested in authenticating raw materials and food products in order to satisfy food quality and safety requirements [1–4]. In particular, assurance of the quality of roasted coffees has attracted widespread attention as a means for controlling and preventing coffee adulteration, and also given the great difference in the final sale price depending on a wide range of factors, including coffee varietal and geographic origin. The two varieties of
∗
Corresponding author. Tel.: +34 941299626; fax: +34 941299621. E-mail address: [email protected] (C. Pizarro).
0003-2670/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.aca.2006.12.057
economic significance in the global coffee trade are Coffea arabica L. (arabica coffee) and Coffea canephora Pierre (robusta coffee) [5]. Both species differ not only in relation to their botanical characteristics and physicochemical composition, but also in terms of commercial value, with arabica coffees commanding market prices 20–25% higher and being considered to be of better quality than robustas because of their superior taste and aroma. However, most commercially available coffees are actually blends of the two species, making it even more difficult to tackle coffee adulteration issues. Therefore, suitable methods are required in order to discriminate between coffee varieties and to detect potential adulterations of high quality coffee beans with poorer and cheaper types, thus ensuring authenticity, quality, safety and efficacy of final products to be commercialised. Extensive research has been carried out in recent years to develop reliable and specific coffee authentication methods [6–24]. Many of the recently developed approaches for determining coffee authenticity have focused mainly on coffee identification and classification on the basis of different types
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
of compositional data thanks to the application of different pattern recognition techniques [6–16]. Despite the relative success achieved by many of these approaches, it is important to consider that many analytical reference methods used to assess the chemical components to be later used as discriminant parameters may be quite expensive, elaborate and/or time-consuming. For this reason, simpler and faster methods, such as those based on spectroscopic techniques which can be easily implemented in routine analysis, have emerged as a very attractive and useful alternative tool for varietal identification purposes [17–24]. In particular, several studies have described the effective application of near infrared spectroscopy (NIRS) to address the problem of coffee authentication [20–24]. In this study we decided to adopt a less explored NIR-based quantitative approach to coffee authentication, namely a combination of NIR spectroscopy and chemometrics to optimize calibration procedures in order to both detect coffee adulteration in blends and the actual adulteration degree. Certain effects that occur inherently in diffuse reflectance NIR spectroscopy for solid samples, such as light scattering and influence of particle size, may be responsible for unexpected perturbations in spectra (baseline shifts, slope changes and non linearity) [25]. As a result, NIR signals contain not only chemical but also physical information about samples and measuring conditions, which may be irrelevant and may mask the significant information in the spectra and deteriorate the quality and predictive ability of the calibration models developed from raw spectra. Therefore, the application of suitable pre-processing methods, aimed at correcting spectral data by minimizing the contribution of physical effects to NIR signals and thus enhancing the relevant chemical information contained, is a crucial step in regression model development and improvement. Additionally, the huge amount of overlapping and redundant information present in NIR spectra often causes the useful information contained to be spread out over multiple points and at regions that are greatly affected by noise and interference. For this reason, NIR spectra are suitable for the application of efficient signal de-noising and compression methods, in order to remove both noise and not significant information from spectra. Consequently, the objective of this study was precisely to examine the feasibility of NIR spectroscopy, in combination with multivariate calibration techniques, to develop improved and reliable regression models which could later be used as a quick and accurate scanning tool for detecting coffee adulterations by quantifying the actual content of robusta variety in roasted coffee samples. We propose the application of a preprocessing method, recently developed by our research group, consisting of orthogonal wavelet correction (OWAVEC) coupled to the calibration step. The OWAVEC method was developed, and subsequently improved [26–28], to simultaneously accomplish two essential needs in multivariate calibration – signal correction and feature selection (data compression) – in an attempt to enhance the reliability and robustness of the finally constructed PLS regression model by combining the versatility that wavelet analysis has shown for signal processing [29–38] and the application of a particular orthogonal signal correction algorithm to remove from spectra information not related to the
267
response variable to be modelled and predicted. In particular, the last version of the OWAVEC method [28], which implemented a genetic algorithm as a feature selection method for extracting the significant corrected wavelet coefficients to be later used in PLS model development has been applied since it enabled a very substantial improvement in method performance. In order to better evaluate the performance of OWAVEC method applied to mixture resolution in roasted coffee samples, results provided by the final OWAVEC-PLS models constructed were compared with those obtained when working on both original data and spectra filtered by other orthogonal signal correction methods (OSC and DOSC). It should be noted that the purpose of this study was not to provide a definitive and immutable NIR calibration model for predicting robusta variety content in coffee samples, but to present a preliminary study to propose a strategy capable of accomplishing this task and for verifying its reliability and effectiveness in order to assess the genuineness of coffee samples. Although the data set used in this work was designed to cover, as best as possible, the great natural variability inherent in commercially available coffees by considering different roasting conditions and degrees, and varied geographical origins, it is quite clear that the dynamic nature of the coffee market and the particular needs and production lines of a given coffee company would require a more exhaustive or specific collection of a suitable set of calibration samples to serve as the basis for developing, using the methodology proposed here, a calibration model for use in the routine quality control of coffee samples. 2. Materials and methods 2.1. OWAVEC method This paper does not aim to describe in detail the fundamentals of all pre-processing methods that were tested to remove NIR spectral variation not related to the response variable and to enhance spectral selectivity, moreover taking into account that some of these methods (first derivative, OSC and DOSC) are quite well known and usually applied in multivariate calibration. Nevertheless, considering the recent introduction of OWAVEC method for signal pre-processing, a brief discussion about its main application features will be provided in order to facilitate a clear understanding of the procedure carried out by OWAVEC. Fig. 1 shows a schematic view of the four main sequential steps into which the OWAVEC algorithm can be broken down. As a first step, OWAVEC applied the discrete wavelet transform (DWT) to decompose each spectrum in the wavelet domain, using a selected wavelet function in this transformation. In DWT, the choice of different wavelet basis functions gives rise to different decompositions of a signal in the resulting wavelet domain. In the present study, and taking into account our experience in wavelet analysis of NIR signals, the biorthogonal Daubechies-4 wavelet was selected, partly because its two vanishing moments provide high compressibility for smooth signals such as NIR data. The DWT requires the length of the analysed signal to be dyadic (2n , where n is an integer). However, if this constraint is not fulfilled several methods can be applied
268
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
Fig. 1. Schematic view of all the steps involved in the OWAVEC method. (1) Centered spectra are padded to the nearest length 2n . Each spectrum is transformed using DWT and Daubechies-4 as wavelet basis function. (2) The wavelet coefficients matrix is deflated of information unrelated to the considered response (orthogonalization). (3) Wavelet coefficients matrix compression by using GAs for subset selection of significant wavelet coefficients. (4) Corrected signals can be reconstructed into the wavelength domain using IWT. The final PLS regression model was developed from the resulting corrected and compressed wavelet coefficients matrix.
in order to extrapolate the signal. In the OWAVEC wavelet decomposition step, all signal extrapolations were carried out by linear padding in order to minimise edge effects. Another key parameter to be optimised when applying DWT for regression purposes is the optimal wavelet decomposition level. Thus, in order to select the most suitable wavelet decomposition to be performed, several decomposition levels were tested, always respecting the maximum allowed wavelet decomposition level associated with the extended length of the analysed signals and the use of Daubechies-4 as wavelet function (in this case, the maximum decomposition level permitted was six). The final choice of the optimal decomposition level was based on several criteria, trying to both maximize the compression rate achieved and the predictive ability of the calibration model developed once the OWAVEC method had been fully executed. Once the wavelet coefficients matrix corresponding to the calibration set had been computed, an orthogonalization algorithm was applied to deflate it of information not related to the response variable (i.e. the robusta content of coffee samples). Further details on the orthogonalization procedure implemented in OWAVEC can be found in reference [26]. Once the orthogonalization algorithm had been applied on the calibration wavelet coefficients matrix, and the spectra in the prediction set had been properly decomposed by wavelet analysis, the computed matrix of weights of the original wavelet coefficients in the orthogonal subspace of these coefficients (W) could be used directly to obtain the corrected wavelet coefficients matrix for test (new) data. A critical parameter to be optimized when performing wavelet coefficients matrix orthogonalization is the value of the tolerance factor used to calculate the generalized inverse of the wavelet coefficients matrix. The use of a generalized inverse to
obtain the matrix of weights W implies the loss of the complete orthogonality constraint as an attempt to only take into account stable directions in wavelet coefficients matrix for computing its orthogonal subspace, which means that some correlation between the information removed from the wavelet coefficients matrix and the response variable is allowed. The adopted degree of deviation from orthogonality will depend precisely on the used value for the tolerance factor, so that for each data set studied the tolerance factor should be properly tuned. Thus, different tolerance factor values were tested in each case, so that maximizing the predictive power of the regression model subsequently developed from the resulting corrected and compressed wavelet matrix, while not showing any noticeable signs of overfitting, was finally selected as the optimal value. Finally, as far as the OWAVEC data compression (feature selection) step was concerned, before searching for significant subsets of corrected wavelet coefficients by GAs, the settings required in the genetic algorithm run need to be selected. The effect of varying the values of some GAs design parameters (number of chromosomes; number of elitist chromosomes; percentage of mutation; mean, maximum and minimum numbers of coefficients to be retained in the same chromosome; maximum number of cycles without fitness function (RMSECV) improvement) was evaluated to favour algorithm convergence and to search for an improved solution. The final choice of the optimal wavelet decomposition parameters, tolerance factor and GAs settings was based on cross-validation results, also seeking to reach as high a compression rate as possible without losing model performance. Moreover, for interpretative and comparative purposes, the spectra corrected after pre-processing by OWAVEC were recon-
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
structed into the original wavelength domain from the respective corrected and compressed wavelet coefficients matrices by inverse wavelet transform (IDWT), for both calibration and prediction sets. This can be very useful for identifying the spectral regions responsible for a successful prediction. Since the only pre-processing procedure performed during OWAVEC application was mean-centering, the wavelet coefficients matrix retained the original variance of input spectra, so after the optional reconstruction step of the filtered signals the global shape of the starting spectra would be maintained. 2.2. Samples Eighty three pure variety roasted coffee samples from varied origins (36 arabica and 47 robusta coffees), which were processed under different roasting conditions, were used for this study, in such a way that they were representative of the variability present in coffee market. In addition, a number of varietal blends were prepared in the laboratory from these pure variety coffees. Taking into account that the preparation of all possible blends resulting from mixing the available roasted coffees from different varieties would have been impractical, we decided to select the samples most representative of both arabica and robusta in order to later combine them to obtain suitable blends covering the overall variability of the whole set of blends potentially derived from the pure coffees considered. The application of PCA on mean centered NIR spectra provided an effective and easy-to-implement tool for selecting representative samples from arabica and robusta varieties. Fig. 2 shows a bidimensional representation of PC1 and PC2 scores accounting for 88.33% of the variance in the roasted coffee NIR spectral data, labelled according to their coffee variety: (1) arabica coffees; (2) robusta coffees. Two sample groups appeared slightly separated by the first bisectrix of the two component axes, suggesting the presence of two different clusters just associated with the two varieties considered. Thus, the centroid of each class was selected as representative of each variety and the two extreme
269
samples within each cluster were selected in order to take into account the variability associated to each group. Next, the three samples selected as representative for each variety were combined to make suitable coffee blends with robusta content in the final blends ranging from 0 to 60% (w/w), in such a way that a total number of 108 blends were consequently generated. This specific range of blending (based on the gradual addition of increasing robusta amounts in final blends up to 60%) was particularly studied in order to cover the actual range of adulteration encountered, which is more specifically located within relatively low and still financially profitable adulteration levels. Fraudulent practice in the elaboration of coffee blends from arabica and robusta pure coffee varieties is largely due to economic reasons, i.e. to reduce costs due to the lower price of robusta coffee as compared to arabica coffee, labelling mixed variety coffees as pure arabica or commercialising blends which actually contain a lower percentage of arabica than that declared. Therefore, the resulting data set used in the present study comprised a total of 191 roasted coffee samples (83 pure varieties and 108 blends). In order to develop and validate the regression models, the entire data set was split into two independent subsets: a calibration set with 133 coffees (comprising 70% of samples and covering the whole variability of the original data set to ensure its representativity); and a test set with the remaining 58 samples. Special care was taken to ensure suitable composition of the external test set, checking that samples of both pure varieties and all different compositional blends were included. 2.3. Apparatus and software NIR spectra were recorded on a NIRSystems 5000 near infrared spectrophotometer (FOSS, NIRSystems) equipped with a reflectance detector and a sample transport module. The instrument was controlled by a compatible PC and Vision 2.22 (FOSS, NIRSystems) was used to acquire the data. The OWAVEC method was entirely programmed, including the data compression step by GAs, using MATLAB® technical computing language, version 6.5 [39]. The steps of the OWAVEC procedure involving wavelet analysis were implemented using the Wavelet Toolbox version 2 (for use with MATLAB® ) [40]. The orthogonal signal correction (OSC) and direct orthogonal signal correction (DOSC) pre-processing methods were applied simply to serve as a comparative reference for evaluating the results provided by OWAVEC. Thus, the OSC method developed by Wise and Gallagher [41] (available in PLS-Toolbox 2.1 for use with MATLAB® [42]), and the MATLAB® code of the DOSC method developed by Westerhuis et al. [43] were used for calculations. For the construction of all the calibration models, PLS regression was performed using the PLS program contained in the V-Parvus package [44]. 2.4. Recording of NIR spectra
Fig. 2. Scores of the 83 roasted coffee samples from arabica (labelled as 1) and robusta (labelled as 2) pure varieties on the first two principal components explaining the variability in the NIR spectral data.
Reflectance spectra were obtained directly from untreated samples. Due care was taken to ensure that the same amount of
270
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
sample was always used to fill up the sample cell. Each spectrum was obtained from 32 scans performed at 2 nm intervals within the wavelength range 1100–2500 nm, with five replicates for each individual sample. The samples were decompacted between recordings. An average spectrum was subsequently obtained from the replicates collected for each coffee sample. 2.5. Data analysis Several data pre-processing methods (mean-centering, first derivative, orthogonal signal correction, OWAVEC or a combination of several of these options) were applied before calibration development in order to find models with as high a predictive power as possible and to compare OWAVEC method performance with that achieved using other pre-treatments. Fig. 3 shows a schematic diagram of the different steps executed for the different pre-processing methods tested. In this way, various PLS regression models were constructed between calibration NIR spectra pre-processed by the different pretreatments used and the robusta variety content of roasted coffee samples (when the OWAVEC method was applied, the corrected and compressed wavelet coefficients matrix was considered for calibration development instead of using NIR spectra in the original wavelength domain). All PLS models were developed from
the column-wise centered data matrix, regardless of whether or not an additional pre-processing method was also applied on spectral variables. The first derivative of spectra was applied in some cases in order to improve spectral differences and to correct baseline effects. In order to avoid noise enhancement as a consequence of derivative, data were smoothed by the Savitzky–Golay method, using a cubic degree polynomial with a window size of 13 points. The optimal complexity of every calibration model was assessed by cross-validation (all the PLS models were built by cross-validation using five deletion groups). The actual predictive ability of all the constructed regression models was also validated by testing their performance on an external test set, in order to control and avoid possible overfitting (particularly when applying orthogonal signal correction methods due to the high risk of overfitting inherent in them). When OSC was applied as a pre-processing method, before the mean-centering step, spectra were transformed into their first derivative spectra, since this additional step substantially improved the quality of the resulting models. When the orthogonal signal correction methods were applied, in order to search for models with high predictive ability and to avoid over-fitted solutions, the suitable orthogonal correction degree to be applied was determined as follows: the number of orthogonal-LVs (OSC) or orthogonal-PCs (DOSC)
Fig. 3. Flow diagram summarising the different pre-processing methods (or combination of methods) applied to the data set.
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
to be removed from spectra was varied from 1 to 6, in such a way that the optimal number of orthogonal components to be subtracted was chosen according to the results obtained in the cross-validation procedure. The root mean square error (RMSE) of the residuals obtained (termed RMSEC in calibration, RMSECV in cross-validation and RMSEP in external prediction) was used to evaluate and compare the respective performance of the different PLS models constructed, since this parameter represented an objective measurement of the resulting predictive ability of the model, which was dimensionally comparable to the studied response. In this particular calibration study, the RMSE values were already obtained in percentages, since the range of robusta variety content in coffee samples varied precisely from 0 to 100%. Likewise, the correlation plots between measured and computed/predicted response values, as well as their corresponding coefficients of determination (R2 ), were also obtained to serve as additional tools to evaluate models accuracy. Nevertheless, the prediction ability of a regression model must be accompanied by an estimate of its uncertainty. Thus, prediction intervals at 95% confidence were computed for each regression model (in both calibration and external prediction) in order to better estimate and compare their respective reliability. 3. Results and discussion 3.1. Regression models for quantifying robusta variety content from roasted coffee NIR spectra Table 1 shows the calibration and prediction errors prompted by the PLS models developed from mean-centered and first derivative NIR spectra to quantify the robusta content of roasted coffee samples. In spite of the high predictive ability achieved when working on mean-centered data (2.98% RMSEP), a large number of latent variables (10 LVs) had to be used to develop the model in order to compensate for the structured noise (baseline and scatter effects) present in the original NIR spectra. A high model complexity could be seen as a limitation for the use of the corresponding calibration model as an accurate tool for detecting coffee adulteration, since it would be desirable to use a small number of components providing a minimal error for a simpler, and probably more robust, calibration model. The transformation of raw signals into their first derivative spectra, prior to calibration model construction, provided a significant decrease in prediction error with respect to the model developed from original spectra (2.04% RMSEP), but at the expense of maintaining high model complexity (8 LVs), which suggested that
271
the contribution to spectra of variation not due to the parameter of interest was still significant. The results corresponding to the calibration models developed from NIR spectra corrected by both OSC and DOSC pre-processing methods (after selection of the optimal number of orthogonal factors to be removed in each case) are also included in Table 1. Both OSC and DOSC methods were applied on NIR spectra taking into account, in order to correct signals, the percentage of robusta variety contained in each sample as response variable. The corrected resulting spectra were then used to develop the respective PLS models. Considering all the models based on signals corrected by both methods, it was observed that the calibration errors and the corresponding model complexities decreased once each orthogonal component was removed from the spectra. Nevertheless, this apparent gradual improvement in model performance was limited for both correction methods, since in the case of OSC the removal of more than two orthogonal latent variables reduced the calibration error at the expense of increasing the prediction error (thus giving an extremely overfitted solution), whereas when applying DOSC a correction beyond five orthogonal PCs did not provide any substantial improvement. The calibration models finally selected as optimal after applying both methods provided similar results with regard to predictive ability (1.36–1.33% RMSEP for OSC and DOSC methods, respectively), and represented a considerable improvement compared to the original model in terms of both model reliability and complexity. However, it should be noted that, despite carefully selecting the number of orthogonalLVs to be removed, OSC led to an overfitted solution since the prediction error almost doubled the corresponding calibration error. Despite the quality improvement prompted by the models developed by applying OSC and DOSC methods, the aim of this study was to search for a suitable pre-processing method capable of both performing efficient signal correction and effective data compression, in such a way that it would be possible to obtain an improved final regression model that would be less complex, more robust and parsimonious and with an increased predictive ability, and which could be applied for the accurate detection of potential coffee adulterations. The major advantage of the OWAVEC method is its ability to compress in few variables (the selected subset of corrected wavelet coefficients) the information responsible for a successful prediction. This will likely not only improve model performance with regard to prediction accuracy, but also make the resulting model more parsimonious, i.e. reduce the dimensionality of the predictors (coefficients) matrix prior to calibration, thus providing a more robust model.
Table 1 Calibration (RMSEC) and prediction (RMSEP) errors and percentages of explained variance provided by the original PLS model (constructed from mean-centered data), the PLS model developed from first derivative spectra and the PLS models built after applying OSC and DOSC methods to NIR spectra
Centering First derivative First derivative + OSC (2 O-LVs removed) DOSC (5 O-PCs removed)
PLS-LVs
%RMSEC
%EXPL.VAR. (cal)
%RMSEP
%EXPL.VAR. (test)
10 8 3 1
2.13 1.67 0.77 1.03
99.67 99.80 99.96 99.92
2.98 2.04 1.36 1.33
99.38 99.71 99.87 99.88
Only the model yielding the best performance for each method is shown.
272
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
Table 2 Calibration (RMSEC) and prediction (RMSEP) errors and percentages of explained variance for PLS models developed from the corrected and compressed coefficients matrix obtained after the application of OWAVEC, when working on both original and first derivative spectra LVs
%RMSEC
%EXPL.VAR. (cal)
%RMSEP
%EXPL.VAR. (test)
LEVEL
TOL
MED/MAX/MIN
CROM/ELIT/CP
NCOF
%ZERO COEFF.
OWAVEC (centered data) 1 1.18 99.90 1 1.17 99.90 1 1.18 99.90 1 1.18 99.90 1 1.18 99.90
1.36 1.34 1.33 1.37 1.36
99.87 99.87 99.88 99.87 99.87
2 2 2 2 2
0.0015 0.0015 0.0015 0.0015 0.0015
8/10/2 4/10/2 4/6/2 5/8/2 5/7/3
20/2/5 20/2/5 20/2/5 20/2/5 20/2/5
4 4 3 4 5
99.44 99.44 99.58 99.44 99.30
First derivative + OWAVEC 2 0.73 99.96 2 0.66 99.97 2 0.68 99.97 2 0.67 99.97 2 0.64 99.97
0.83 0.81 0.79 0.84 0.81
99.95 99.95 99.96 99.95 99.95
2 2 2 2 2
0.0005 0.0005 0.0005 0.0005 0.0005
7/10/2 7/10/5 8/10/5 8/10/3 8/10/4
20/2/5 20/2/5 20/2/5 25/1/5 20/2/5
9 10 10 10 10
98.74 98.60 98.60 98.60 98.60
For each pre-treatment combination tested, the results obtained after performing several runs using different choices of GA parameters are shown. The particular OWAVEC-GA parameters used in each run (including optimal wavelet decomposition level, tolerance factor to compute the generalized inverse, number of chromosomes, number of elitist chromosomes, percentage of mutation, mean, maximum and minimum numbers of coefficients allowed for selection in the same chromosome and maximum number of cycles) are also included. Resulting compression scores (number of retained coefficients and percentage of zero coefficients) are shown.
Fig. 4. Plot of NIR computed (in cross validation) vs. reference robusta contents for calibration samples. Results from the PLS model based on: (a) original data (only mean centered); (b) first derivative spectra; (c) OSC corrected first derivative spectra; (d) DOSC filtered spectra; (e) OWAVEC corrected and compressed wavelet coefficients matrices computed from original data; (d) OWAVEC corrected and compressed wavelet coefficients matrices working on first derivative spectra. Only a model representative of all constructed models in the different GA runs performed when applying OWAVEC method was shown. The 95% confidence interval of the linear fitting was shown (grey dashed line).
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
Due to the stochastic character of the optimisation procedure carried out by OWAVEC-GA, a preliminary explanatory study, using different choices of GA parameters settings, was performed to find an optimal search region, thus to obtain a number of selected coefficients large enough to achieve a good fit but small enough to avoid overfitting. The agreement observed between results obtained after different GA runs for each data pre-treatment tested indicated a good convergence of the genetic algorithm applied. Table 2 summarizes the results obtained with the different calibrations models developed from the respective subsets of OWAVEC corrected and compressed wavelet coefficients selected by genetic algorithms in five different runs of the method performed when working on both original and first derivative spectra. All very similar PLS models developed after the direct application of the OWAVEC-GA method to original data provided high predictive ability (≈1.35% RMSEP) with the minimum possible complexity and no appreciable signs of overfitting (≈1.18% RMSEC), which represented a substan-
273
tial improvement with respect to the models constructed from spectra not processed by orthogonal signal correction methods. Moreover, the impressive compression rates resulting from the feature selection step carried out by GAs (99.30–99.58% zero coefficients), served to demonstrated the enhanced effectiveness of the OWAVEC-GA method with respect to the other orthogonal signal correction methods tested here, since although no substantial improvement in predictive ability was achieved, the final model improved substantially in terms of parsimony (only 3–5 from a total of 713 corrected wavelet coefficients were selected in the different runs) and overfitting was reduced compared to OSC based model. Finally, we tried to further improve the reliability and stability of OWAVEC-based PLS models by taking advantage of the proven ability of first derivative to reduce baseline offset variations in NIR signals and to increase spectral resolution. Thus, Table 2 also includes the results corresponding to the regression models constructed from the respective corrected and compressed wavelet coefficients matrices resulting from
Fig. 5. Plot of NIR predicted (in external validation) vs. reference robusta contents for test samples. Results from the PLS model based on: (a) original data (only mean centered); (b) first derivative spectra; (c) OSC corrected first derivative spectra; (d) DOSC filtered spectra; (e) OWAVEC corrected and compressed wavelet coefficients matrices computed from original data; (d) OWAVEC corrected and compressed wavelet coefficients matrices working on first derivative spectra. Only a model representative of all constructed models in the different GA runs performed when applying OWAVEC method was shown. The 95% confidence interval of the linear fitting was shown (grey dashed line).
274
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
the application of OWAVEC on mean-centered first derivative NIR spectra, with the performance of several GA-searches. All PLS models developed after the coupled application of the first derivative and OWAVEC methods yielded substantially improved calibration and prediction results for assessing the robusta content of roasted coffees not only with respect to the models constructed from original data, but also with regard to the models developed from data processed by both other signal correction methods and to the OWAVEC method when working directly on mean-centered data. In addition to the extremely low prediction errors achieved (0.79–0.84% RMSEP) with low model complexities (only 2 LVs) and lack of overfitting (0.64–0.73% RMSEC), special consideration must be given to the high compression rates obtained thanks to the implementation of GAs as a feature selection method in OWAVEC (only 9–10 from a total of 713 wavelet coefficients were retained; 98.60–98.74% zero coefficients), which allowed us to obtain a more parsimonious calibration model. We could also observe that the value of the tolerance factor used in the OWAVEC orthogonalization step when working on first derivative spectra was lower than for raw spectra. This fact made clear that the
treatment of NIR spectra by first derivative prior to OWAVEC application served to effectively minimize background effects, since the previous correction for unwanted systematic variation achieved by derivative resulted, to a certain a extent, in a reduction of non-stable directions in the corresponding wavelet matrix, thus enabling us to consider a lower correlation between the discarded information and the response variable. In order to gain further insight into the accuracy of the different calibration models developed and to compare their relative goodness, linear regression of the computed/predicted values for the robusta content versus reference values, for both calibration and test sets, was applied (Figs. 4 and 5). The grey dashed lines are the prediction bands at 95% confidence limit, in such a way that the smaller was the prediction interval, the model precise was the model. Thus, of all spectral pre-processing methods which were tested to optimize the calibration model for detecting potential coffee adulterations, the coupled application of first derivative and OWAVEC not only gave rise to the model displaying the lowest complexity, RMSEC and RMSEP, but the highest correlation for both the calibration and prediction sets (R2 = 1) together with the lowest uncertainty in prediction (the most nar-
Fig. 6. Spectral profiles of both calibration and test roasted coffee samples corresponding to (a) Original NIR spectra. (b) First derivative spectra. (c) OSC corrected first derivative spectra. (d) DOSC corrected spectra. (e) OWAVEC corrected reconstructed spectra. (f) OWAVEC corrected reconstructed first derivative spectra.
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
row prediction intervals). Moreover, to provide an even more precise idea of the level of accuracy achieved with the proposed model, the 95% confidence interval for the prediction error might be estimated as ±1.98 RMSEP, i.e. ±1.60% adulteration. This high accuracy achieved in the resolution of arabica/robusta mixtures was even more remarkable considering that the detection limit of NIR spectroscopy is about 0.1–1%. 3.2. Observations on NIR spectra The scatter effects occurring inherently in NIRS can be significant and produce a considerable expansion of the absorbance interval for the individual wavelengths, which can often hinder the extraction of significant information from the raw data. This phenomenon was clearly revealed in the spectral profiles corresponding to original NIR signals shown in Fig. 6a. For this reason, NIR spectra need to be pre-processed in order to remove systematic noise. First derivative pre-treatment did not remove all the spectral variations caused by the scatter effects, or rather, it only reduced them at some spectral regions (Fig. 6b). One of the main objectives that orthogonal signal correction methods attempted to accomplish was the minimization of systematic effects, by removing information not related to the response variable, in order to improve the quality of the resulting calibration models. In fact, as can be observed from the spectra obtained after data filtering by the OSC (Fig. 6c) and DOSC (Fig. 6d) methods, the spectral differences were significantly minimized with respect to the raw signals, proving the relative efficiency of both methods to correct data, trying to preserve only spectral information closely related to the robusta content of coffee samples. In the previous sub-section, we demonstrated the usefulness and effectiveness of the OWAVEC-GA method for correcting and compressing data into a wavelet domain prior to the development of a reliable and robust calibration model, as well as its better performance potential with respect to the other preprocessing methods tested, mainly by expounding numerical results. However, a graphical analysis of the spectral profiles reconstructed from the respective subsets of wavelet coefficients finally selected by OWAVEC can also help to demonstrate the goodness of the calibration strategy proposed here. These reconstructed profiles can be used to extract significant information about which specific spectral regions, representing the underlying features responsible for a successful prediction of the robusta content of roasted coffees, were associated with the corrected coefficients selected for the development of the final regression model. As can be seen, from plots displaying both original (Fig. 6e) and first derivative (Fig. 6f) reconstructed spectra after pre-processing by OWAVEC, the spectral differences between samples were not only minimised but restricted to a very limited number of regions since exclusively information related to robusta content was selected to be kept. 4. Conclusions The results reported in this study show that a method based on the combination of NIR spectroscopy and multivariate cal-
275
ibration can be used to detect and quantify potential coffee adulterations. However, the success of this methodology will depend largely on signal pre-processing methods applied to minimize the spectral variation not due to the parameter of interest but due to variation in experimental or sample conditions. The regression models developed for predicting the robusta variety content of coffee samples, after pre-processing first derivative NIR spectra by OWAVEC, were considerably better compared to the models obtained from raw data (in terms of both reliability and parsimony) and improved results were also obtained with respect to the other orthogonal signal correction methods tested. The OWAVEC-GA selection of a small subset of corrected wavelet coefficients responsible for accurate modelling not only led to an improvement of predictive ability but also to a decrease of the complexity associated with the resulting model, making it more parsimonious and probably more robust, since spectral regions exhibiting a large variance in raw spectra, as well as noise, not related to the robusta content of roasted coffee samples were discarded. Although the work reported here is only a feasibility study and further research should be carried out before its value for authentication purposes can be truly estimated, the high quality calibration models constructed, thanks to the high effectiveness of OWAVEC as a correction and compression tool, suggest that the proposed regression methodology would be a suitable tool for the quantitative detection of robusta contamination of arabica coffees in the authenticity assessment of coffee. Acknowledgements The authors thank the Ministry of Science and Technology (Project No. 2FD1997-0491) for its financial support, as well as Professor Michele Forina for providing us with the last version of Parvus package. References [1] P.R. Ashurst, M.J. Dennis (Eds.), Food Authentication, Blackie Academic & Professional, London, UK, 1996. [2] R.S. Singhal, P.R. Kulkarni, D.V. Rege, Handbook of Indices of Food Quality and Authenticity, Woodhead Publishing Ltd., Cambridge, UK, 1997. [3] P.R. Ashurst, M.J. Dennis (Eds.), Analytical Methods of Food Authentication, Blackie Academic & Professional, London, UK, 1998. [4] M. Lees (Ed.), Food Authenticity and Traceability, Woodhead Publishing Ltd., Cambridge, UK, 2003. [5] R. Viani, A. Illy (Eds.), Espresso Coffee: The Science of Quality, second ed., Academic Press, London, UK, 2004. [6] N. Frega, F. Bocci, G. Lercker, Industrie Alimentari 34 (1995) 705. [7] C.P. Bicchi, M.P. Ombretta, G. Pellegrino, A.C. Vanni, J. Agric. Food Chem. 45 (1997) 4680. [8] M.J. Mart´ın, F. Pablos, A.G. Gonz´alez, Talanta 46 (1998) 1259. [9] F. Carrera, M. Le´on-Camacho, F. Pablos, A.G. Gonz´alez, Anal. Chim. Acta 370 (1998) 131. [10] F. Pablos, A.G. Gonz´alez, M.J. Mart´ın, M.S. Valdenebro, M. Le´onCamacho, Analyst 124 (1999) 999. [11] M.J. Mart´ın, F. Pablos, A.G. Gonz´alez, Food Chem. 66 (1999) 365. [12] M.J. Mart´ın, F. Pablos, A.G. Gonz´alez, M.S. Valdenebro, M. LeonCamacho, Talanta 54 (2001) 291. [13] A.G. Gonz´alez, F. Pablos, M.J. Mart´ın, M. Le´on-Camacho, M.S. Valdenebro, Food Chem. 73 (2001) 93.
276
C. Pizarro et al. / Analytica Chimica Acta 585 (2007) 266–276
[14] S. Casal, M.R. Alves, E. Mendes, M.B.P.P. Oliveira, M.A. Ferreira, J. Agric. Food Chem. 51 (2003) 6495. [15] C.G. Zambonin, L. Balest, G.E. De Benedetto, F. Palmisano, Talanta 66 (2005) 261. [16] L. Mondello, R. Costa, P.Q. Tranchida, P. Dugo, M. Lo Presti, S. Festa, A. Fazio, G. Dugo, J. Sep. Sci. 28 (2005) 1101. [17] E.K. Kemsley, S. Ruault, R.H. Wilson, Food Chem. 54 (1995) 321. [18] R. Briandet, E.K. Kemsley, R.H. Wilson, J. Agric. Food Chem. 44 (1996) 170. [19] R. Briandet, E.K. Kemsley, R.H. Wilson, J. Sci. Food Agric. 71 (1996) 359. [20] G. Downey, J. Boussion, D. Beauchˆene, J. Near Infrared Spectrosc. 2 (1994) 85. [21] G. Downey, J. Boussion, J. Sci. Food Agric. 71 (1996) 41. [22] G. Downey, B. Spengler, J. Agric. Food Res. 35 (1996) 179. [23] G. Downey, R. Briandet, R.H. Wilson, E.K. Kemsley, J. Agric. Food Chem. 45 (1997) 4357. [24] I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz, C. Pizarro, Anal. Chim. Acta 514 (2004) 57. [25] B.G. Osborne, T. Fearn, P.H. Hindle, Practical NIR Spectroscopy with Applications in Food and Beverage Analysis, second ed., Longman Scientific and Technical, Harlow, UK, 1993. [26] I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz, C. Pizarro, Anal. Chim. Acta 515 (2004) 31. [27] I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz, C. Pizarro, Anal. Chim. Acta 544 (2005) 89. [28] I. Esteban-D´ıez, J.M. Gonz´alez-S´aiz, D. G´omez-C´amara, C. Pizarro, Anal. Chim. Acta 555 (2006) 84. [29] R.R. Coifman, Y. Meyer, M.V. Wickerhauser, Wavelet analysis and signal processing, in: M.B. Ruskai, G. Beylkin, R. Coifman, I. Daubechies, S.
[30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]
Mallat, Y. Meyer, L. Raphael (Eds.), Wavelets and Their Applications, Jones & Bartlett, Boston, USA, 1992, pp. 153–178. B. Walczak, D.L. Massart, Trends Anal. Chem. 16 (1997) 451. B.K. Alsberg, D.B. Kell, J.J. Rowland, M.K. Winson, A.M. Woodward, Analyst 122 (1997) 645. B. Walczak, D.L. Massart, Chemometr. Intell. Lab. Syst. 36 (1997) 81. V.J. Barclay, R.F. Bonner, I.P. Hamilton, Anal. Chem. 69 (1997) 78. A.K.M. Leung, F.T. Chau, J.B. Gao, Chemometr. Intell. Lab. Syst. 43 (1998) 165. J. Trygg, S. Wold, Chemometr. Intell. Lab. Syst. 42 (1998) 209. S. Mallat, A Wavelet Tour of Signal Processing, second ed., Academic Press, London, UK, 1999. K. Jetter, U. Depczynski, K. Molt, A. Niem¨oller, Anal. Chim. Acta 420 (2000) 169. B. Walczak (Ed.), Wavelets in Chemistry, Elsevier Science B.V., The Netherlands, 2000. MATLAB® 6.5, The MathWorks, Natick, USA, 2002. M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi, Wavelet Toolbox for Use with MATLAB® , The MathWorks, Natick, USA, 2000. B.M. Wise, N.B. Gallagher, http://www.eigenvector.com/MATLAB/ OSC.html. B.M. Wise, N.B. Gallagher, PLS Toolbox 2.1, Eigenvector Research Inc., USA, 1998. J.A. Westerhuis, S. de Jong, A.K. Smilde, Chemometr. Intell. Lab. Syst. 56 (2001) 13. M. Forina, S. Lanteri, C. Armanino, C. Cerrato-Oliveros, C. Casolino, VPARVUS 2004, an extendable package of programs for explorative data analysis, classification and regression analysis, Dip. Chimica e Tecnologie Farmaceutiche ed Alimentari, University of Genova.