Chemometric Strategies for Spectroscopy-Based Food Authentication

Featured Application

This review will offer a global overview of the chemometric approaches most commonly used in the field of spectroscopy-based food analysis and authentication. Three different scenarios will be surveyed: data exploration, calibration and classification. Basic and simple descriptions of the main multivariate techniques exploited in such a domain along with a comprehensive outline of their most recent and interesting applications will be provided.

Abstract

In the last decades, spectroscopic techniques have played an increasingly crucial role in analytical chemistry, due to the numerous advantages they offer. Several of these techniques (e.g., Near-InfraRed—NIR—or Fourier Transform InfraRed—FT-IR—spectroscopy) are considered particularly valuable because, by means of suitable equipment, they enable a fast and non-destructive sample characterization. This aspect, together with the possibility of easily developing devices for on- and in-line applications, has recently favored the diffusion of such approaches especially in the context of foodstuff quality control. Nevertheless, the complex nature of the signal yielded by spectroscopy instrumentation (regardless of the spectral range investigated) inevitably calls for the use of multivariate chemometric strategies for its accurate assessment and interpretation. This review aims at providing a comprehensive overview of some of the chemometric tools most commonly exploited for spectroscopy-based foodstuff analysis and authentication. More in detail, three different scenarios will be surveyed here: data exploration, calibration and classification. The main methodologies suited to addressing each one of these different tasks will be outlined and examples illustrating their use will be provided alongside their description.

1. Introduction

In recent years, consumers’ attention towards the quality of foodstuff has become increasingly lively. The awareness that there exists a close link between health and diet has spread dramatically, leading a growing number of people to develop conscious eating habits. For this reason, and thanks to the so-called food revolution, aliments (especially high value-added food products) are commonly subjected to strict quality controls. These tests are of paramount importance, especially for attesting to some peculiar features (connected, for instance, to their geographical origin, their specific manufacturing process, and/or the know-how of their producers) that might translate into the recognition of distinctive labelling designations like the Protected Designation of Origin (PDO) or the Protected Geographical Indication (PGI). In light of this, a plethora of analytical methodologies aimed at foodstuff authentication and traceability have recently been developed, and a wide variety of applications of such methodologies have been reported in the scientific literature. More specifically, given the notable market value of this type of product, much effort has been put into the design of non-destructive characterization tools for quality control, generally based on the principles of light–matter interaction and spectroscopy. However, unfortunately, the signal profiles yielded by spectroscopy instrumentation (regardless of the spectral range investigated) are often particularly complex and thus, their assessment and interpretation are usually not straightforward. In this context, multivariate statistical approaches (also known as chemometric approaches) have lately played—and currently still play—a crucial role: they have proven, in fact, to be extremely powerful when large amounts of spectral data (characterized by a considerable degree of intercorrelation among the recorded spectral variables) need to be coped with, and when useful and meaningful information is to be extracted from these data for disparate purposes.

Given, therefore, the sheer relevance of such approaches in the field of spectroscopy-based foodstuff analysis and authentication, this review aims at providing a comprehensive overview of those most commonly exploited in real-world case studies. Specifically, three scenarios will be here surveyed: exploration, calibration (necessary to carry out quantification) and classification (for, e.g., adulteration and/or fraud detection). The first one, intended to pursue direct insights into the data with the possibility of revealing hidden/underlying structures and relationships between samples and/or variables, and making extensive use of plots and graphs to highlight similarities, differences, trends, clusters and/or correlations, relies on so-called unsupervised techniques, which do not require any input other than the (spectroscopic) data themselves [1]. On the other hand, calibration and classification problems call for the use of predictive models that exploit the spectroscopic information to predict one or more properties of the objects under study; in order to be reliable, such models have to be constructed by supervised approaches, i.e., techniques which actively take advantage not only of the experimental (spectroscopic) data but also of, e.g., reference values of the aforementioned properties or sample labelling [2]. The main chemometric methodologies suited to addressing each one of these different tasks will be outlined, and examples illustrating their use will be provided alongside their description.

2. Data Exploration

Prior to conducting a proper food authentication study, users might be interested in carrying out preliminary investigations to assess, for instance, whether a particular spectroscopic technique is capable of discerning products of distinct geographical origins or whether an instrumental response is sensitive enough to detect the presence of compounds or adulterants of interest in the specific samples at hand. In this regard, chemometric exploratory tools based on the reduction of the dimensionality of the original data collected are of crucial importance, as they dramatically ease the visualization of these data, allowing possible differences among specimens (related to one or multiple of the aforementioned aspects, like the presence of undesired/unexpected substances) to be spotted in a straightforward way. Given, in fact, the complex multidimensional and multivariate nature of spectroscopic measurements, such differences are very often undiscernible when simply plotting raw spectral profiles since they normally lie not in few scattered regions of the concerned wavelength domain, but in the correlation structure intrinsic to the recorded spectral variables or descriptors (this is the case especially for low-selectivity approaches as Near-InfraRed—NIR—spectroscopy). For this reason, one of the most commonly exploited tools in similar scenarios to facilitate and improve the identification of diverse groups of food samples is Principal Component Analysis (PCA) [3,4]. PCA makes it possible to compress, describe and interpret large sets of multidimensional data. Its basic principle can be summarized as follows: let $X$ be a centered $N\times J$ matrix with $J$ denoting the number of variables (e.g., $J$ wavelengths of light scanned in a spectroscopy experiment) registered, for instance, at $N$ time instants or for $N$ different objects. As specified, when $J$ is very large, the useful and meaningful information in $X$ is usually intercorrelated among various of these variables over the whole set of recordings. Then, for a chosen degree of acceptable accuracy, it is possible to reduce the $J$-dimensional space of the original descriptors to an $A_{\mathrm{PCA}}$-dimensional subspace in which data mostly vary and onto which all the $N$ samples under study can be projected and represented as new points. Mathematically speaking, PCA is based on the bilinear structure model in Equation (1):

$$X=T{P}^{\mathrm{T}}+E_{\mathrm{PCA}}$$

(1)

where $P$ ($J\times A_{\mathrm{PCA}}$) is an array of so-called loadings, which determine the $A_{\mathrm{PCA}}$ basis vectors (principal components or factors) of the PCA subspace, $T$ ($N\times A_{\mathrm{PCA}}$) defines the projection coordinates or scores of all the $N$ rows of $X$ on this lower-dimensional space and $E_{\mathrm{PCA}}$ ($N\times J$) denotes the matrix of unmodelled residuals, i.e., the portion of $X$ not explained at the chosen rank, $A_{\mathrm{PCA}}$—in PCA, the optimal number of principal components to extract can be estimated by a wide range of approaches that, broadly speaking, can be classified into three distinct categories [5]: ad hoc rules (like Kaiser's eigenvalue-greater-than-1 rule [6], Velicer’s minimum average partial rule [7], and Cattell’s scree test [8]), statistical tests (like Bartlett’s Chi-square test [9] and Tracy-Widom’s statistics-based test [10]) and computational criteria (like cross-validation [11,12,13] and permutation testing [14,15,16,17,18,19,20,21].

The PCA solution may be formulated in many equivalent ways and attained by a variety of algorithms, among which the most widespread and popular is certainly Singular Value Decomposition (SVD [22]). SVD decomposes $X$ as:

$$X=UD{V}^{\mathrm{T}}+E$$

(2)

with the columns of $U$ ($N\times A_{\mathrm{PCA}}$) and $V$ ($J\times A_{\mathrm{PCA}}$) being the first $A_{\mathrm{PCA}}$ left and right singular vectors of $X$, respectively, and $D$ ($A_{\mathrm{PCA}}\times A_{\mathrm{PCA}}$) a square diagonal array whose diagonal elements correspond to its first $A_{\mathrm{PCA}}$ non-zero singular values. Therefore, it holds that $T=UD$ and $P=V$.

PCA shows the following property:

$$P^{\mathrm{T}}P=I$$

(3)

where $I$ is an identity matrix of dimensions $A_{\mathrm{PCA}}\times A_{\mathrm{PCA}}$, which translates into the fact that the dimensions of its subspace are orthogonal and, thus, can be inspected assuming that the information they individually capture is mutually unrelated.

The assessment and interpretation of a PCA model is generally performed by examining both scores and loadings profiles (i.e., the columns of $T$ and $P$, respectively). The former yield insights into existing relations among analyzed objects, while the latter enable the identification of correlation patterns within ensembles of multiple spectral variables. The basic idea, here, is to simultaneously analyze these profiles not only to possibly pinpoint the presence of distinctive clusters of samples, but also to recognize sets of spectral features shared within each one of these clusters and responsible for their differentiation. Nevertheless, as a consequence of the orthogonality imposed to the principal components, PCA scores and loadings seldom provide physico-chemically meaningful information, since they do not capture single phenomena or events affecting the registered measurements, but rather heterogeneous effects resulting from combinations of such phenomena or events. To overcome this limitation, alternative approaches like Independent Component Analysis (ICA [23,24]) can be resorted to. Unlike PCA, ICA regards the rows of the data matrix $X$ as a collection of observed signals that are mixtures of a certain number of common unknown source signals or independent components and decomposes it as:

$$X=A{R}^{\mathrm{T}}+E_{\mathrm{ICA}}$$

(4)

with $A$ ($N\times A_{\mathrm{ICA}}$) being the so-called proportion matrix, carrying the mixing coefficients associated with these independent components contained in every column of $R$ ($J\times A_{\mathrm{ICA}}$), and $E_{\mathrm{ICA}}$ ($N\times J$) denoting the ICA residual array.

As for PCA, the interpretation of an ICA model rests on the simultaneous inspection of the retrieved independent components and their respective proportion profiles. Independent components are analogous to PCA loadings, but, differently from them, are not necessarily orthogonal, and therefore, can be easily related to individual phenomena impacting the recorded spectral variables (e.g., Rayleigh or Raman scattering) [25,26,27,28]. Subsequently, ICA proportions can give a quantitative idea of how they influence the raw measured spectra.

There exist several algorithms by which the ICA solution can be attained. Examples are Infomax [29], FastICA [30] and JADE (Joint Approximate Diagonalization of Eigenmatrices) [26,31]. Their basic principle, though, is common: source signals are estimated by maximizing their statistical independence. In essence, they can be looked at as the outcome of a rotation of PCA loadings that enhances their aforementioned interpretability properties.

PCA and ICA have been both extensively used in the field of food characterization (see

Table 1. Application examples of exploratory, regression and classification analysis of spectroscopic data in the field of food science. PDO, PCA, ICA, PLS, LDA, PLSDA, LW-PLSDA, K-PLSDA and QDA stand for Protected Designation of Origin, Principal Component Analysis, Independent Component Analysis, Partial Least Squares, Linear Discriminant Analysis, Partial Least Squares Discriminant Analysis, Locally Weighted Partial Least Squares Discriminant Analysis, Kernel-Partial Least Squares Discriminant Analysis and Quadratic Discriminant Analysis, respectively.

Aim	Method	Reference
Exploration of apple varieties	PCA	[32]
Beer storage monitoring	PCA	[33]
Adulteration detection of Camellia oils	PCA	[34]
Chemical characterization of Mediterranean olive oils	ICA	[35]
Chemical characterization of honey samples	ICA	[36]
Chemical characterization of soft drinks	ICA	[36]
Determination of lard content in chocolate samples	PLS	[92]
Quantification of turmeric adulteration in egg-pasta	PLS	[93]
Straw wine quality parameter prediction	PLS	[94]
Sensory characterization of Trentingrana cheese	PLS	[95]
Egg content quantification in dried egg-pasta	Local PLS	[76]
Characterization of PDO Chianti Classico olive oil	LDA	[123]
Determination of the geographical origin of pistachios	PLSDA	[124]
Classification of rice varieties	LW-PLSDA/K-PLSDA/Artificial neural networks	[117,125]
Classification of honey samples	PLSDA	[126]
Technological classification of egg white powders	PLSDA	[127]
Discrimination of distillates	PLSDA	[128]
Classification of tomato genotypes	LDA/PLSDA/Support vector machines/	[129]
Olive fruit classification	QDA	[130]
Insect infestation detection in stored rice	PLSDA	[131]
Characterization of Italian craft beers	PLSDA	[132]

); by way of illustration, PCA has been coupled to (i) high-resolution ¹H Nuclear Magnetic Resonance (NMR) spectroscopy for the discrimination of three apple varieties (Bramley, Russet and Spartan), making it possible to spot significant variations in their malic acid and sucrose content [32]; (ii) synchronous fluorescence spectroscopy for monitoring changes in the flavin composition of beer during storage either under light exposure or in darkness [33]; and (iii) Attenuated Total Reflectance Mid-InfraRed (ATR-MIR) spectroscopy and diffuse reflectance NIR spectroscopy in a feasibility study on the detection of soybean oil adulteration in Camellia oil samples [34]. Conversely, due to its particular methodological basis, ICA has mainly been applied to extract the individual spectral contributions and proportion profiles of certain classes of analytes constituting heterogeneous food matrices, e.g., primary oxidation and polyphenolic compounds, tocopherol, carotenoids and chlorophylls in extra virgin olive oils [35], or fructose, sucrose and glucose in soft drinks [36].

Sometimes, in the context of unsupervised data analysis, one may be interested in having information about the different degrees of similarity/dissimilarity among samples or, more specifically, on whether the distribution of individuals enables the identification of groups or clusters. Clusters correspond to sets of observations which are more similar to one another than they are with respect to all the remaining objects. In this regard, the search for possible clusters of samples is the objective of the so-called clustering techniques, whose basic idea is that similarity among individuals can be inversely related to a distance measure in a multivariate space [37,38]. Based on this concept, clustering techniques can operate non-hierarchically (partitioning methods) or hierarchically. The former [39], encompassing also the popular $k$-means algorithm [40], distribute the samples among a pre-defined number of groups (that can be either known a priori or inferred empirically from the data under study); a sample is basically placed in the cluster with the closest barycenter. They are relatively straightforward, but they suffer from the limitation of needing the number of groups, $k$, to be defined beforehand.

Conversely, as the name suggests, hierarchical approaches proceed iteratively through successive agglomerations of smaller clusters into larger ones [41]—the same approach can be operated top down by progressively dividing larger clusters into smaller ones, but it is more rarely used. Every object is initially assumed to constitute a separate group of observations (singleton) and, afterwards, at each iteration, the two most similar groups are joined into a single ensemble. The procedure continues until a certain stopping criterion is met or when all the samples are gathered into a unique cluster. Hierarchical approaches provide more detailed insights into the relations between individuals and/or groups of individuals, but they have the disadvantage that, due to the nature of the agglomeration process, the identification of the most plausible number of clusters may not be easy. The interested reader can find more details on clustering techniques in [37,38,42].

3. Calibration

The analytical characterization of a foodstuff, especially with the aim of its authentication, may imply the quantitative determination of one or more constituents of the samples at hand or of one or more of their global properties, such as the sensory scores for some specific attributes, the iodine value or the dry matter content. Moreover, one would ideally want such a determination to be rapid (in some cases, almost instantaneous, to allow controls to be conducted on- or in-line), non-invasive or at least non-destructive, not needing the specimens to be subjected to clean-ups or pre-treatments, and, possibly, not encompassing the use of auxiliary reagents or organic solvents. All these requirements point in the direction of spectroscopic techniques, which, depending on the applications concerned, possess most or all of the aforementioned characteristics. However, the possibility of using spectroscopy to quantify the properties of a food product is conditional, on the one hand, to the existence of a relation between the measured signal and the response to be estimated and, on the other hand, to the postulation of a reasonable mathematical formulation to express (or at least approximate) this relation. In other words, quantifying the value of a property, $y$, based on the measurement of a spectroscopic signal, $x$, rests on inferring a functional relation, $f$, that connects them so that:

$$y=f\left(x\right)$$

(5)

Unfortunately, the exact form of the function $f\left(x\right)$ linking a measured spectrum to the property to be determined is unknown and cannot be retrieved from first principles. A reasonable guess of $f\left(x\right)$ should be therefore empirically drawn, through a procedure which is called calibration [43,44,45].

Calibration makes use of known data to approximate the functional relation in Equation (5). In particular, it requires that, for a sufficient number of samples (which are usually gathered in a so-called training set), both the recorded spectra and the values of the property(-ies) of interest are available. For example, if NIR spectroscopy is used for the quantification of the protein content in wheat, in order to proceed with the calibration phase, one should not only collect the NIR spectra of the training samples, but also determine their protein content by, e.g., reference methods such as Kjeldahl [46] or Dumas [47]. Afterwards, a mathematical expression of the function $f\left(x\right)$ is explicitly assumed: this function will depend not only on $x$, but also on the values of some adjustable parameters (coefficients), whose estimation is the core issue in calibration. In many cases, $f\left(x\right)$ can be assumed to be linear; this means that, if $J$ spectral variables are recorded (i.e., if $x$ is a $J$-dimensional vector of light intensity/absorption measurements, with elements $\left[\begin{array}{cc}{x}_1& x_2\end{array},\begin{array}{ccc}{x}_3& \cdots & x_J\end{array}\right]$), Equation (5) can be written as:

$$y=\widehat{y}+e=b_1{x}_1+b_2{x}_2+b_3{x}_3+\cdots +b_J{x}_J+e$$

(6)

where $\widehat{y}$ is the approximation of $y$ provided by the linear function defined by the coefficients $b_1,b_2,b_3,\cdots ,b_J$, and $e$ is the residual, i.e., the difference between the true value of $y$ and $\widehat{y}$. In calibration, the available ($x,y$) pairs constituting the training set of data are used to calculate the coefficients in Equation (6) through regression, according to some user-defined criterion of optimality. Since the aim of calibration is to use the measured signal $x$ as the basis to obtain the best prediction of $y$, in the majority of cases such a criterion involves the minimization of the residuals (least squares). In particular, if the training set contains $N$ samples for which both $x$ and $y$ are known, for each one of these samples an equation analogous to Equation (6) can be written:

$$y_n={\widehat{y}}_n+e_n=b_1{x}_{n1}+b_2{x}_{n2}+b_3{x}_{n3}+\cdots +b_J{x}_{nJ}+e_n$$

(7)

where the subscript $n$ indicates that the equality holds for the $n$-th training object—here, it should be noticed that, since the same functional relation is valid for all the samples, the coefficients $b_1,b_2,b_3,\cdots ,b_J$ are identical for all the individuals. By gathering all the training set recordings in the spectral matrix $X$ (whose rows contain the profiles of all the individuals) and the corresponding true values of the property of interest in the column vector $y$, Equation (7) can be globally expressed for all the training samples as:

$$y=\widehat{y}+e=Xb+e$$

(8)

where the column vector $b$ carries the regression coefficients ($b=\left[\begin{array}{cc}{b}_1& b_2\end{array},\begin{array}{ccc}{b}_3& \cdots & b_J\end{array}\right]$), the approximated values of the response are the elements of $\widehat{y}\left(=Xb\right)$, and the residuals constitute the vector $e$. As anticipated, the estimation of the optimal value of the coefficients $b$ is usually carried out by the so-called least squares approach:

$$\underset{b}{\mathrm{argmin}}{{\displaystyle \sum }}_{n=1}^N{e}_n^2=\underset{b}{\mathrm{argmin}}{e}^{\mathrm{T}}e=\underset{b}{\mathrm{argmin}}{\left(y-Xb\right)}^{\mathrm{T}}\left(y-Xb\right)$$

(9)

where the operator $\mathrm{argmin}$ indicates that the optimal set of values is the one minimizing the overall approximation error across all the training objects, i.e., ${{\displaystyle \sum }}_{n=1}^N{e}_n^2=e^{\mathrm{T}}e={\left(y-Xb\right)}^{\mathrm{T}}\left(y-Xb\right)$. Resolving Equation (9) makes it possible to estimate the regression coefficients as:

$$b_{\mathrm{MLR}}={\left(X^{\mathrm{T}}X\right)}^{-1}{X}^{\mathrm{T}}y$$

(10)

Since the calibration approach, which postulates a linear regression model and uses the least squares criterion for the calculation of these regression coefficients, is also generally referred to as Multiple Linear Regression (MLR [48,49]), the subscript $\mathrm{MLR}$ was explicitly added to Equation (10).

Once the values of the coefficients have been estimated using the training samples, it is possible to predict the value of the response for any new (unknown) specimen, ${\widehat{y}}_{\mathrm{new}}$, based on its spectral profile, $x_{\mathrm{new}}$, as:

$${\widehat{y}}_{\mathrm{new}}=x_{\mathrm{new}}{b}_{\mathrm{MLR}}$$

(11)

Sometimes, the same spectroscopic profile can be used to estimate the values of multiple properties; for instance, one might want to determine, for the same individual, the protein, lipid, starch and moisture content, based on its measured NIR spectrum. Mathematically, this would result in setting up a calibration equation analogous to Equation (8) for each of the $L$ responses to be predicted:

$$y_l={\widehat{y}}_l+e_l=X{b}_l+e_{l}l=1,\dots ,L$$

(12)

In the case of MLR, it is possible to demonstrate that there is no difference (in terms of estimated values of the regression coefficients $b_l$ and, subsequently, of the predicted responses ${\widehat{y}}_l$) between building separate models for all the single properties of interest and constructing a unique global calibration [50] that can be expressed as:

$$Y=\widehat{Y}+E=XB+E$$

(13)

where each of the $L$ columns of the matrices $Y$, $\widehat{Y}$, $B$ and $E$ contains the values of the corresponding elements in Equation (12) for a particular response:

$$\begin{gathered} Y=\left[\begin{array}{ccccc}{y}_1& \cdots & y_l& \cdots & y_L\end{array}\right] \\ \widehat{Y}=\left[\begin{array}{ccccc}{\widehat{y}}_1& \cdots & {\widehat{y}}_l& \cdots & {\widehat{y}}_L\end{array}\right] \\ B=\left[\begin{array}{ccccc}{b}_1& \cdots & b_l& \cdots & b_L\end{array}\right] \\ E=\left[\begin{array}{ccccc}{e}_1& \cdots & e_l& \cdots & e_L\end{array}\right] \end{gathered}$$

(14)

MLR has the advantages of a clear loss-function, which directly generalizes the least squares criterion to the multivariate case, and of not requiring the tuning of any metaparameter (that is to say, once the $X$-$Y$ training pairs are identified, the solution is univocal) [45]. However, due to the mathematical structure of Equation (10), this method cannot be applied in situations where the matrix $X$ is ill-conditioned, i.e., when the number of training samples is lower than the number of recorded variables (wavelengths) and/or when such variables are highly correlated [51]. Unfortunately, both these conditions are commonly met by spectroscopic data, and thus, MLR is seldom applied if spectral profiles are concerned, unless some variable selection strategy is adopted [52,53,54,55].

On the other hand, a solution to the problem of having to deal with ill-conditioned matrices of descriptors is provided by the use of bilinear models based on the extraction of principal components or latent variables [56]. Indeed, as already outlined in Section 2—where the use of principal components to summarize the relevant information encoded in the data was discussed—an effective dimensionality reduction can be achieved by identifying a reduced number of highly informative orthogonal directions in the multivariate space of the original variables and projecting the data onto the subspace spanned by them, obtaining a set of new coordinates or scores, $T$. Using $T$ instead of $X$ in a calibration framework makes it possible to overcome all the issues related to ill-conditioning, since the number of relevant components is usually lower than the number of samples and because, additionally, they are completely uncorrelated due to their orthogonality. In light of these aspects, various bilinear regression techniques have been proposed in the literature, the most commonly used of which are Principal Component Regression (PCR [48,57,58]) and Partial Least Squares (PLS [59,60,61,62,63,64,65,66]) regression.

PCR, as the name suggests, is a two-step method which is based on the sequential use of PCA and MLR, with the PCA scores, $T$, exploited as predictors:

$$Y={\widehat{Y}}_{\mathrm{PCR}}+E=T\mathcal{B}+E$$

(15)

Here, $\mathcal{B}$ is the matrix of regression coefficients relating the scores to the predicted responses through the MLR model. By resorting to the mathematical relationship defining the PCA projection ($T=XP$, see also Equation (1)), it is possible to express the linear model in Equation (15) directly as a function of $X$ and $Y$, as:

$${\widehat{Y}}_{\mathrm{PCR}}=T\mathcal{B}=XP\mathcal{B}=X{B}_{\mathrm{PCR}}$$

(16)

where the matrix of PCR coefficients, $B_{\mathrm{PCR}},$ is given by:

$$B_{\mathrm{PCR}}=P\mathcal{B}$$

(17)

With respect to the MLR solution reported in Equation (10), PCR is said to be biased, as not all the information originally present in $X$ is used to estimate the parameters of the regression model [67,68]; the calculation of $B_{\mathrm{PCR}}$, in fact, relies solely on the variability captured by the PCA scores in $T$, and a different model will be returned depending on how many principal components are extracted from the data matrix $X$. Accordingly, the number of principal components to retain in a PCR model is a metaparameter which has to be optimized; normally, models with different numbers of principal components are tested and the one resulting in the lowest prediction error (usually computed by cross-validation or other resampling procedures) is selected [69,70,71] (see also Section 7.2).

PCR can be considered a highly performant statistical approach as it enables the possibility of calculating regression models for ill-conditioned data matrices. Nevertheless, even in similar contingencies it may not necessarily constitute the best option to choose: indeed, the two steps PCR combines (i.e., data compression by PCA and regression by MLR) have different objectives, and it is not always verified that the directions of maximum explained variance in $X$ are also those corresponding to the maximum correlation with the response(s) [72]. Conversely, PLS tries to overcome this issue, by actively using the information in $Y$ already at the data compression stage, so that the scores extracted from $X$ are relevant for describing at the same time the variance in the descriptors and in the properties of interest.

PLS rests on the extraction of two sets of scores, one from the independent and one from the dependent data block, having maximum covariance with one another:

$$\underset{r_a,q_a}{\mathrm{argmax}}{u}_a^{\mathrm{T}}{t}_{a,\mathrm{PLS}}=\underset{r_a,q_a}{\mathrm{argmax}}{q}_a^{\mathrm{T}}{Y}^{\mathrm{T}}X{r}_a$$

(18)

where $t_{a,\mathrm{PLS}}$ and $u_a$ are the $X$- and $Y$-scores along the $a$-th component, respectively, and $r_a$ and $q_a$ represent the vectors of coefficients allowing to retrieve such scores from the corresponding original matrices:

$$\begin{array}{c}{t}_{a,\mathrm{PLS}}=X{r}_a\\ u_a=Y{q}_a\end{array}$$

(19)

The functional relation (usually a linear dependence), also known as inner relationship, between the two data blocks is then implemented at the scores level; component by component, the $Y$-scores are approximated by the $X$-scores through a univariate linear regression model:

$${\widehat{u}}_a=t_{a,\mathrm{PLS}}{c}_a$$

(20)

with $c_a$ being a proportionality constant (inner regression coefficient). Since more than a single component is usually needed to summarize the relevant information in $X$ and $Y$, Equations (19) and (20) can be rewritten as:

$$\begin{array}{c}{T}_{\mathrm{PLS}}=XR\\ U=YQ⟹\widehat{Y}=U{Q}^T\\ \widehat{U}=T_{PLS}C\end{array}$$

(21)

where the columns of the scores matrices $T_{\mathrm{PLS}}$ and $U$ and those of the coefficient matrices $R$ ($X$-weights) and $Q$ ($Y$-loadings) are associated with the individual retrieved components (see also Equation (19)), and $C$ is a diagonal array whose non-zero elements are the coefficients $c_a$. The overall regression model can then be expressed as:

$$\widehat{Y}=U{Q}^T=T_{\mathrm{PLS}}C{Q}^T=XRC{Q}^T=X{B}_{\mathrm{PLS}}$$

(22)

Equation (22) shows that, despite the regression model be calculated at the level of the scores, due to the linearity of the projection, it is possible to express it directly in terms of the original variables (similarly to PCR, see also Equation (16)). At the same time, since only a part of the information in the original data is used for the regression, and, in particular, for the estimation of the regression coefficients:

$$B_{\mathrm{PLS}}=RC{Q}^T$$

(23)

PLS constitutes another case of so-called biased regression [67,68], requiring the determination of the optimal number of components to be retained in the model (which, once again, is carried out by selecting the value leading to the minimum prediction error, usually in cross-validation).

A fundamental aspect of PLS regression is that, contrarily to MLR and, therefore, to PCR, when the problem involves the calibration of multiple responses, individual models yield different outcomes compared to a single global model. Since PLS is also based on the extraction of latent variables from the $Y$-block, in fact, the calculation of a single model for calibrating all the responses at the same time entails their sharing part of their systematic variability, i.e., that they will exhibit a certain degree of intercorrelation that is not due to chance or subsampling. Only in this case is it recommended to build a single model for predicting all the responses; if this condition is not met, then it is suggested that each response be calibrated separately [73].

Altogether, MLR, PCR and PLS represent the large majority of regression methods which are used in scenarios demanding calibration based on spectroscopic data. However, all of them return linear models, while circumstances where the complexity of the problem at hand may call for non-linear functional relations might also occur. Given the nature of spectral data, a relatively popular way of implementing a non-linear functional relation—if necessary—is to use locally linear models [74,75,76]. The concept behind local regression is that a problem which is globally non-linear may be approximated by a combination of locally linear regressions, each one covering a relatively small portion of the variable space. Other options are non-linear extensions of the PLS algorithm [77,78,79,80,81,82], support vector machines encompassing non-linear kernel transformations [83,84,85,86], and artificial neural networks [87,88,89,90,91] (which, nonetheless, have more stringent requirements regarding the number of training objects needed for constructing a reliable model).

Throughout the years, regression analysis has played and still plays a crucial role in the domain of food characterization. In particular, PLS, together with its non-linear variants, has been widely used for addressing a large number of issues in such a context (see Table 1), from the estimation of lard content in chocolate specimens by ATR-Fourier Transform-IR (ATR-FT-IR) spectroscopy [92] to the quantification of egg content and turmeric adulteration in egg-pasta by NIR spectroscopy [76,93], the estimation of quality parameters of straw wine by FT-MIR spectroscopy (alcohol and sugar content and acidity) [94], and the analytical profiling of sensory attributes of Trentingrana cheese by Proton Transfer Reaction–Mass Spectrometry (PTR-MS) [95]. This last study, in particular, was conducted as an attempt to relate specific flavor and odor features to characteristic compounds detected and identified by the PTR-MS instrumental platform (the perception of boiled milk flavor, for instance, was found to be connected to the presence of methanetiol in the samples at hand).

4. Classification

Apart from exploration and regression, another task that users and practitioners commonly perform while conducting food authentication studies is classification. Classifying a sample or object implies predicting one or more of its discrete properties based on the information collected during its characterization [96]. An example would be the determination of the geographical origin of a certain product from its spectral profile recorded within a certain wavelength range (UltraViolet-Visible—UV-Vis—NIR, MIR, etc.). More specifically, classification would aim at assigning such a product to one category or class constituted by objects sharing similar features (in this case, objects sharing the same geographical origin). From a slightly different perspective, considering an individual sample observation as a vector (e.g., a spectrum) corresponding to a point in the multivariate space of the experimental variables (e.g., the wavelength channels), classification approaches can also be regarded as tools for the identification of boundaries in this space separating the various categories at hand: the aforementioned sample is therefore predicted as belonging to a particular class when the respective point falls within its associated boundaries.

Given the extreme relevance of similar strategies in food authentication scenarios and the importance they can have for the resolution of real-world issues, it does not come as a surprise that the applications of existing classification methodologies along with the proposal and development of novel ones have dramatically increased in the last decades [97,98]. Although these methodologies may differ in complexity, requirements and assumptions, two broad groups of techniques can be distinguished: discriminant and modelling techniques [99].

4.1. Discriminant Techniques

Discriminant approaches are probably the most classical solutions for dealing with classification-related problems [97]. As their name suggests, they make it possible to highlight differences between samples belonging to distinct classes. Their basic principle is simple: they divide the multivariate space of the registered variables into several subregions equal to the number of categories of objects considered (say $Z$) and they assign each one of such objects to a certain class if the point corresponding to its measurement vector falls within the boundaries associated with that class. Subsequently, owing to these characteristics:

• discriminant classification models need to be calibrated on training sets composed by specimens belonging to all the categories under study;
• every analyzed sample is always assigned to one and only one of these categories;
• samples coming from other classes (not considered in the study) will always be (erroneously) recognized as members of one and only one of these categories.

For all these reasons, discriminant techniques can be particularly useful, e.g., when the amount of expected classes is limited and/or if it can be reasonably supposed that all the samples to be assessed are drawn uniquely from them.

The mother of all discriminant classification strategies is undoubtedly Linear Discriminant Analysis (LDA) [100]. In the simplest case encompassing two categories of objects, LDA estimates a direction of maximum separation between classes (i.e., a so-called canonical variate) as:

$$w^{\mathrm{T}}={\left(m_1-m_2\right)}^{\mathrm{T}}{S}^{-1}$$

(24)

with $m_1$ and $m_2$ (both of dimensions $J\times 1$) being the vectors of variable means (centroids) for class #1 and class #2, respectively, and $S$ ($J\times J$) an estimate of the class covariance matrix (which, here, is assumed to be the same for both categories (this assumption is relaxed in the most direct extension of LDA known as Quadratic Discriminant Analysis (QDA) [101,102])). $m_1$, $m_2$ and $S$ are computed from training samples. A class delimiter, $w_0$, is then calculated as:

$$w_0=w^{\mathrm{T}}\frac{\left(m_1+m_2\right)}{2}$$

(25)

and the classification rule established such that if the projection coordinate of each observation vector, say $x_n$ ($J\times 1$), along $w$ ($t_n=w^{\mathrm{T}}{x}_n$) is found to be lower than $w_0$, the corresponding object is predicted as belonging to class #1 and vice versa. The rationale behind LDA can also be interpreted by assuming that data within each category are Gaussianly distributed so as to calculate the so-called posterior probability of $x_n$ belonging to the z-th class, $p(z|x_n)$, as:

$$p\left(z|x_n\right)=\frac{p_0\left(z\right)}{{\left(2\pi \right)}^{\frac{J}{2}}\left|S\right|}{e}^{-\frac{1}{2}{\left(x_n-m_z\right)}^{\mathrm{T}}{S}^{-1}\left(x_n-m_z\right)}$$

(26)

with $p_0\left(z\right)$ equal to the probability of observing an individual from the same category before carrying out any measurement (prior probability). Accordingly, the final classification rule can be reformulated by stating that a sample should be assigned to the class it has the highest posterior probability of proceeding from.

Even though LDA can be in principle extended to scenarios involving a larger number of categories, it requires the inversion of the covariance matrix $S$, which can be singular or nearly singular in the presence of highly collinear (inter-correlated) descriptors (as for spectral profiles) [51]. To overcome this limitation, one can perform a proper selection of these descriptors (via, e.g., step-wise algorithms) [103], apply LDA to the scores resulting from a preliminary PCA modelling of the collected data in what is called PCA-LDA [101,104,105] or resort to regularized versions of LDA, like in Regularized Discriminant Analysis (RDA) [106,107,108,109,110,111]. While, on the one hand, such approaches are capable of stabilizing the LDA solution against collinearity, on the other hand, some of them might suffer from similar drawbacks as seen before for PCR [104,105,112].

Alternatively, one can think of classification as a regression problem in which the class belonging of a sample (i.e., the dependent variable) is to be estimated from the set of (independent) variables returned by a particular analytical platform (e.g., a spectrum). Similarly to what was outlined in Section 3, such an estimation can be carried out by an extension of PLS named Partial Least Squares Discriminant Analysis (PLSDA) [112,113]. PLSDA regresses the entire data matrix $X$ on a dummy binary-coded response array, say again $Y$, made up of a set of piled $Z$-dimensional row vectors, constructed so that, if their corresponding objects/samples are members of the $z$-th class, they have a 1-value in their $z$-th entry and 0-values in all the other ones. For instance, in a 2-class scenario, samples belonging to the first category will be described by the dependent vector [1 0], while samples belonging to the second one by the vector [0 1]. Notice that the model structure in Equation (22) applies also in this case and, thus, whenever new objects/samples become available, their projection coordinates onto the PLS latent variables, as well as their $Y$-predicted values, can be retrieved as shown in Section 2. The class assignation can be finally carried out based on an LDA model constructed either on PLS scores of training objects (this approach is also known as PLS-LDA [112,114,115]) or on their $Y$-predicted values, according to a highest-prediction rule or to a higher-than-a-threshold prediction rule [116]. As for PLS, more complex variants of PLSDA are available for coping with strong non-linear relationships existing between $X$ and $Y$ (e.g., Locally Weighted PLSDA—LW-PLSDA [117]—Kernel-PLSDA—K-PLSDA [118,119,120]). Support vector machines and artificial neural networks can also be used for the resolution of classification-related problems [83,91,121,122].

For the particular type of problems they make it possible to tackle, discriminant methods have had a long history in the field of food quality assessment and control (see Table 1). LDA, QDA, linear and non-linear PLSDA, support vector machines and artificial neural networks have been broadly used to discriminate foodstuff of different origins [123,124], different varieties [117,125,126], different technological properties [127], different purity degree [128], exhibiting different genotypes (i.e., transgenic/non-transgenic) [129] or affected by different diseases [130,131]. PLSDA has been applied, for example, to NIR spectra of pistachios for their geographical authentication making it possible to highlight that seeds from Turkey and USA show a typical absorption behavior between 8000 and 9000 cm⁻¹ (second overtone of methylenic stretching vibrations) [124] and to UV-Vis spectra of ale beers brewed by distinct producers pinpointing minimal variations in the tone and intensity of their color which could indirectly influence customer perception [132].

4.2. Modelling Techniques

Contrarily to discriminant techniques, modelling approaches are capable of capturing similarities among samples belonging to the same category [97,133]. They basically define a multivariate boundary for each considered class, which delimits a specific region of the multidimensional space of the original descriptors where objects proceeding from it are likely to be found. If an analyzed sample falls within this region, it is assigned to the respective class, otherwise it is considered to be a class outlier and rejected as such [134,135].

A fundamental aspect that needs to be taken into account here is that categories are treated and handled separately and that an individual model is constructed for every one of them. Subsequently, (i) modelling methodologies can easily be applied in scenarios in which a unique class of interest exists (a common contingency in food traceability or authentication studies), and (ii) new objects can be predicted as members of one, none or multiple classes also in the light of the fact that the various class spaces do not necessarily have to cover completely the whole original variable space [136].

Given their characteristic nature, class modelling strategies can all be regarded as outlier detection methods. Nonetheless, mainly depending on the outlyingness criterion adopted, such strategies can exhibit different advantages and disadvantages and be more or less suitable to be applied in certain situations rather than others. Due to the flexibility that the selection of this outlyingness criterion guarantees, a plethora of class modelling approaches have been proposed throughout the last 45 years. Among those appearing in the chemometric literature (see

Table 2. Main class modelling methods and application examples in the field of spectroscopy-based food analysis. PGI and PDO stand for Protected Geographical Indication and Protected Designation of Origin, respectively.

Methods
Name	Aim	Reference
UNEQual class spaces (UNEQ)	Exploratory	[137,138]
Soft Independent Modelling of Class Analogy (SIMCA)	Regression	[139,140]
Non-parametric class modelling	Regression	[151]
Neural networks-based class modelling	Regression	[152,153]
Partial Least Squares Density Modeling (PLSDM)		[154]
Potential function (kernel density) method (POTFUN)	Regression	[168,169,170]
Pattern Recognition by Independent Multicategory Analysis (PRIMA)	Regression	[171]
Multivariate Range Modeling (MRM)	Regression	[172]
Support Vector Domain Description (SVDD)		[173]
Applications
Aim	Method	Reference
Traceability of rice varieties	UNEQ/SIMCA/Neural networks-based class modelling	[152]
Authentication of wine samples	UNEQ/SIMCA/MRM/Neural networks-based class modelling	[153,155,156]
Authentication of beer samples	UNEQ/SIMCA/POTFUN	[132,157,158]
Traceability of extra virgin olive oils	SIMCA/Non-parametric class modelling	[151,159,160]
Authentication of olive seeds	UNEQ/SIMCA/PLSDM	[154,161]
Traceability of coffee	SIMCA/POTFUN	[150,162]
Traceability of walnuts	SIMCA	[163]
Authentication of Avola almonds	SIMCA	[164]
Authentication of PGI Gragnano pasta	SIMCA	[165]
Authentication of Italian PDO hazelnut	SIMCA	[166]
Authentication of Vallerano chestnut	SIMCA	[167]
Plant ripening monitoring	PRIMA	[174]
Wheat straw fermentation monitoring	SVDD	[175]

), probably the most commonly exploited since their development are UNEQual class spaces (UNEQ [137,138]) and Soft Independent Modelling of Class Analogy (SIMCA [139,140]).

UNEQ was introduced by Derde and Massart in 1986 and constitutes the modelling version of QDA. Briefly, UNEQ defines the class model by the centroid of the concerned category, while the class space is represented by the multidimensional ellipsoid (hyperellipsoid) corresponding to a user-defined confidence level, i.e., to a certain probability of finding samples of that category within its boundary. Mathematically speaking, in its initial formulation UNEQ calculates for each investigated object ($x_n-J\times 1$) its squared Mahalanobis distance from the $z$-th class centroid ($m_z-J\times 1$) as:

$$d_n^2={\left(x_n-m_z\right)}^{\mathrm{T}}{S}_z^{-1}\left(x_n-m_z\right)$$

(27)

where $S_z$ ($J\times J$) denotes an estimate of the $z$-th class covariance matrix. Once again, $m_z$ and $S_z$ are computed from training objects. If $d_n^2$ is found to be larger than a critical distance value, the respective sample is rejected by the class model as an outlier. Such a critical value can be retrieved as detailed in [137,138,141]. Variants of both the distance statistic and the way its corresponding threshold is determined have been proposed over the years; an alternative UNEQ framework based on the principles of the Hotelling’s T² statistic was, for instance, developed by Forina et al. in 1995 [142].

Despite its relative simplicity, as for MLR and LDA, UNEQ encompasses the inversion of a covariance matrix, which is seldom attainable in case the registered variables show a high degree of intercorrelation [51]—as a solution, UNEQ can be applied to the scores resulting from a preliminary PCA modelling of the data at hand. In addition, UNEQ requires such variables to follow a multivariate Gaussian distribution, an assumption that is rarely met when distinct sources of variability are present in the data at hand (e.g., in the case of the traceability of a designated food product, these sources of variability may be associated with different cultivars, harvesting years, producers, etc.) [133].

On the other hand, in SIMCA, the data associated with the z-th category of samples ($X_z-N_z\times J$, with $N_z$ equal to the number of training objects belonging to the $z$-th class) are modelled by PCA as:

$$X_z=T_z{P}_z^{\mathrm{T}}+E_z$$

(28)

Every object to be classified ($x_n-J\times 1$) is afterwards projected onto the class principal component subspace as:

$$t_{n,z}^{\mathrm{T}}=x_n^{\mathrm{T}}{P}_z$$

(29)

and its distance from such a subspace (of dimensionality $A_z$) is calculated as:

$$s_n^2=\frac{e_{n,z}^{\mathrm{T}}{e}_{n,z}}{\left(J-A_z\right)}$$

(30)

with $e_{n,z}^{\mathrm{T}}=x_n^{\mathrm{T}}-t_{n,z}^{\mathrm{T}}{P}_z^{\mathrm{T}}$.

This distance value is compared with the average distance to the model subspace of the $N_z$ training samples (computed as $s_z^2={{\displaystyle \sum }}_{n_{z=1}}^{N_z}{{\displaystyle \sum }}_{j=1}^J\frac{e_{n_z,j,z}^2}{\left(J-A_z\right)\left(N_z-A_z-1\right)}$, where $e_{n_z,j,z}^2$ results from the projection of the $n_z$-th training object onto the PCA class model subspace) by means of an $F$ test with ($J-A_z$) and ($N_z-A_z-1$) degrees of freedom. If the null-hypothesis of the test is rejected, the corresponding sample is labelled as an outlier and rejected by the class model.

This first implementation of SIMCA was proposed by Wold et al. in 1976 and was almost immediately amended in order to additionally account for the sample distance within the principal component subspace and not only from the principal component subspace [140]. In addition, together with this alternative formulation, many others have been reported in literature generally based on diverse ways of defining the class space and/or the classification rule (see, e.g., [143]).

Finally, it is important to note that the performance of class modelling approaches mainly depends on the settings of several metaparameters such as the confidence level at which the class space is constructed and/or the number of latent variables extracted from the original data. Different approaches for the optimization of these metaparameters exist, whose efficiency and robustness can vary according to the specific scenario faced (e.g., high class overlap, presence of outliers in the training set, etc.) [144,145,146,147,148,149].

For their particular nature and for the main implication resulting from their use (i.e., each category of objects is handled independently and separately), modelling classification strategies are nowadays regarded as ad hoc solutions for coping with food authentication problems (see Table 2). In particular SIMCA has lately been widely resorted to for the traceability of high value-added foodstuffs like wine, beer, extra virgin olive oil, olive seeds, walnuts, coffee and rice [132,150,151,152,153,154,155,156,157,158,159,160,161,162,163], and has been found to guarantee extremely satisfactory performances when PDO and/or PGI products are to be discerned from lower-quality ones. In this regard, several application studies involving the combination of NIR spectroscopy and SIMCA for the characterization of Avola almonds [164], Gragnano pasta [165], Italian hazelnuts [166] and Vallerano chestnuts [167] have been recently reported in literature.

5. Data Fusion

Although all the illustrated methods certainly represent suitable tools to achieve the purposes they are conceived for, sometimes, the complexity of an analytical problem requires to be regarded from several different perspectives to be comprehensively embraced. As an example, many present-day food chemistry issues cannot be easily addressed unless distinct instrumental techniques (e.g., Gas Chromatography–Mass Spectrometry—GC-MS—and MIR spectroscopy) are combined for a more extensive characterization of the samples at hand. Similar strategies are commonly defined multi-block and the chemometric approach used for the concerted analysis of multi-platform data is known as Data Fusion (DF).

Broadly speaking, DF enables the simultaneous extraction of meaningful and useful information from diverse analytical sources. In recent years, DF has been found to provide more exhaustive descriptions of studied systems compared to when single datasets are assessed separately. For this reason, in the last decades, numerous DF methodologies have been proposed in the literature with the three tasks surveyed in the previous sections of this article in mind: data exploration, regression and classification. Nonetheless, despite the ultimate purposes behind their development, such methodologies are usually classified according to the level at which the data fusion is implemented. In the so-called low-level DF, for instance, raw datasets are simply concatenated and analyzed at the same time as a unique global ensemble. Mid-level DF extracts specific features (e.g., principal components or latent variables) from each considered dataset which are subsequently gathered and investigated concurrently. Finally, in high level DF, the outcomes of individual models constructed for every concerned dataset are pooled in a unique solution. Another common distinction among DF strategies is based on the nature of their underlying algorithms (sequential—if joint information is retrieved iteratively—or parallel—if joint information is retrieved collectively) [176].

Due to the large number of existing DF approaches, the following subsections will offer only a basic and brief description of those most commonly used in the field of food authentication. In

Table 3. Additional data fusion methods and application examples in the field of spectroscopy-based food analysis. MB-PLS, SO-PLS, SO-CovSel and OnPLS stand for Multi-Block-Partial Least Squares, Sequential and Orthogonalized-Partial Least Squares, Sequential and Orthogonalized Covariance Selection and n-block Orthogonal Projections to Latent Structures, respectively.

Methods
Name of the Method	Aim	Reference
Hierarchical PLS (H-PLS)	Regression	[178]
Joint and Individual Variation Explained (JIVE)	Exploratory	[213]
Multiblock PLS serial extension	Regression	[214]
Network-Induced Supervised Learning (NI-SL)	Regression	[215]
Parallel Orthogonalized Partial Least Squares (PO-PLS)	Regression	[216]
Multiblock Redundancy Analysis	Regression	[217]
OnPLS	Regression	[218]
Applications
Aim	Multi-Block method	Reference
Prediction of bread sensory properties	MB-PLS	[190]
Prediction of wine ageing time	MB-PLS/H-PLS/NI-SL/SO-PLS	[192]
Quantification of protein and moisture in soybean flour	MB-PLS	[193]
Discrimination of lemon essential oils	MB-PLS	[194]
Determination of the geographical origin of wine	MB-PLS	[196]
Authentication of spirits	SO-PLS/SO-CovSel	[128,202]
Determination of the geographical origin of saffron	SO-PLS/SO-CovSel	[205]
Path modelling	SO-PLS	[219]
Path modelling	SO-PLS	[220]

, a detailed list of applications of DF in this domain is also reported.

5.1. Multi-Block Data Exploration

One of the first exploratory multi-block methods dates back to 1987. It was proposed by Wold and Hellberg and it is called Consensus PCA (CPCA) [177]. CPCA is the natural extension of PCA for the analysis of multiple datasets. Basically, the algorithm searches for directions of maximum variance common to all the data blocks at hand. CPCA encompasses four different steps: (i) one measured variable is selected as super-score; (ii) block loadings and scores are estimated from such a super-score; (iii) block scores are gathered in an updated super-score block which is used to calculate a set of super-weights; and (iv) super-weights are normalized and a new super-score is computed. The procedure is iterated until convergence.

A few years later, in 1996, Wold himself, together with Kettaneh and Tjessem, developed a novel version of multi-block PCA, called Hierarchical PCA (HPCA), sharing similar algorithmic features with CPCA [178]. Their main difference, though, relates to the fact that in HPCA the super-scores are normalized instead of the super-weights. It was observed that both CPCA and HPCA exhibit severe convergence issues, eventually overcome through the modifications proposed by Westerhuis and Coenegracht [179].

Later on, in 2001, a new sequential exploratory approach (Generalized PCA—GPCA [180]) appeared in literature. GPCA constructs a PCA model on the matrix resulting from the concatenation of the datasets under study and extracts normalized super-scores. These super-scores are resorted to for the retrieval of block scores and loadings. The outcomes of a comprehensive comparison of different variants of CPCA, HPCA and GPCA are reported in [181].

Almost concomitantly, in 2000, Qannari et al. implemented Common Components and Specific Weights Analysis (CCSWA) [182,183]. CCSWA seeks the shared information underlying all the investigated datasets. Specifically, $A_{\mathrm{com}}$ orthogonal common components are extracted from a combination of the cross-product matrices corresponding to the $I$ data blocks coped with ($X_i{X}_i^{\mathrm{T}}\forall i\in \left[1,I\right]$), which is weighted according to their so-called salience, i.e., the portion of variance they explain in each one of them. Over the years, different formulations of CCSWA have been designed. CCSWA has nowadays become very well-known with the name ComDim.

More recently, many alternative multi-block techniques (thoroughly reviewed in [184,185]), based on the extraction of common factors, but at the same time also aimed at the recognition of the distinctive information carried by each of the individual considered datasets, have been proposed (e.g., DIStinctive and COmmon Simultaneous Component Analysis—DISCO-SCA [186]—or 2-block Orthogonal Projections to Latent Structures—O2PLS [187]).

5.2. Multi-Block Regression and Classification

The chemometric literature is rich in multi-block methods for regression and classification. One of the most commonly used is undoubtedly Multi-Block-PLS (MB-PLS). MB-PLS constitutes the direct extension of PLS to the data fusion field. Different variants of its algorithm have been proposed all over the years [179], the most recent being developed by Qin et al. in 2001 [188]. The basic idea behind this latest formulation is that multiple predictor blocks can be concatenated and modelled simultaneously by PLS. Although extremely intuitive and yielding accurate predictions, this strategy does not provide information about the single datasets. Conversely, the solution implemented by Westerhuis and Coenegracht can be exploited for the interpretation of block scores and loadings [179]. In this regard, it is important to stress that MB-PLS requires the datasets at hand to be scaled separately (by, e.g., their individual sum-of-squares or their Frobenius’ norm), to prevent structures with greater variance from preferentially driving the model and blurring the information carried by the others [179].

Due to its remarkable flexibility, MB-PLS is nowadays probably the most widely applied multi-block regression approach. It is also often used in combination with LDA for addressing classification problems—MB-PLS is here known as MB-PLS-LDA or MB-PLSDA. So far, many applications of MB-PLS have been reported in the field of food analysis (see Table 3); it has been resorted to (i) for the investigation of sensory parameters of different nature and of their relationships with technological properties of cheese and bread samples [189,190], (ii) for the prediction of meat spoilage time, wine ageing time and crude protein and moisture content in soybean flour by MIR and NIR spectroscopy [191,192,193], (iii) for the discrimination of botanical varieties of extra virgin olive oil, lemon essential oils and wines of different geographical origin by MIR, NIR and Raman spectroscopy [194,195,196], and (iv) for distinguishing added-value from low-quality products [132].

Despite all the aforementioned advantages MB-PLS offers, usually redundant information and/or the presence of categorical regressors may jeopardize the predictive capability of MB-PLS. In these situations, techniques like the recent Sequential and Orthogonalized-PLS (SO-PLS [197,198]) could represent feasible options to overcome such limitations. SO-PLS enables the sequential modelling of each regressor set. The keystone of this methodology is a preliminary orthogonalization step performed to filter out redundancies among blocks. Briefly, taking into consideration the case where a generic response $Y$ is to be estimated from two predictor matrices, say $X_1$ and $X_2$, the SO-PLS algorithmic scheme encompasses the following four steps (readers are addressed to [140,141] for a more detailed description of the SO-PLS algorithm):

• $Y$ is regressed onto $X_1$ by PLS;
• $X_2$ is orthogonalized with respect to the $X_1$-scores calculated in step 1, yielding the array $X_{2,\mathrm{orth}}$. This ensures the common information shared by $X_1$ and $X_2$ is removed from the latter;
• The PLS residuals resulting from step 1 are regressed onto $X_{2,\mathrm{orth}}$ by PLS;
• The SO-PLS predictive model is expressed by combining the outcomes of steps 1 and 3 as:

  $$Y=X_1{B}_1+X_2{B}_2+E_{\mathrm{SO}-\mathrm{PLS}}$$

  (31)

  where $B_1$ and $B_2$ contain the regression coefficients for $X_1$ and $X_2$, respectively, and $E_{\mathrm{SO}-\mathrm{PLS}}$ carry the final $Y$-residuals.

As a consequence of the deflation of the common information shared by the different datasets under study, SO-PLS provides useful insights into the unique and distinctive sources of variation within each one of them. Analogously to MB-PLS, SO-PLS has already been broadly employed in food science (see Table 3), also coupled to LDA for tackling classification tasks [199]; it has been used (i) to correlate sensory perception and chemical composition descriptors [200], (ii) for the quantification of particular compounds in dietary supplements [201], (iii) for the authentication of spirits (i.e., Italian grappa—Italian grape marc spirit) by MIR and NIR spectroscopy [202], and (iv) for the determination of the geographical origin of foodstuff (e.g., red garlic, semolina and saffron) by multiple spectroscopic modalities (UV-Vis, NIR, MIR and ATR-FT-IR) [203,204,205]. Furthermore, the robustness and versatility of SO-PLS has recently led to the design of hybrid techniques conceived for very disparate purposes: Sequential Multi-Block PLS (SMB-PLS [206]), merging features of both MB-PLS and SO-PLS; Sequential and Orthogonalized-N-PLS (SO-N-PLS [207]), for the multi-block analysis of multi-way structures; Sequential Preprocessing through ORThogonalization (SPORT [208]), for the pretreatment of the SO-PLS predictor blocks; Sequential and Orthogonalized Covariance Selection (SO-CovSel [209]), for the low-level DF of preselected variables; and all their discriminant extensions. Yet, one of the main drawbacks associated with SO-PLS regards the fact that a deep and complete inspection of the factors shared by the regressor matrices cannot be easily carried out.

In similar scenarios, P-ComDim, the predictive version of CCSWA, constitutes the designated approach [210,211,212]. The algorithmic procedures behind P-ComDim and CCSWA are strictly resemblant with the difference that, in the former, for each regressor block, the cross-product matrices $X_i{X}_i^{\mathrm{T}}$ are replaced by the core matrices $X_i{X}_i^{\mathrm{T}}Y{Y}^{\mathrm{T}}$. Common components (together with their salience values) are iteratively retrieved as in ComDim, but, here, also $Y$-scores are computed. Additionally, P-ComDim has been rather often applied for tackling complex food authentication problems, not only due to its generally high prediction accuracy, but also because it may open wide horizons for the identification of commonalities underlying regressors proceeding from very diverse sources.

6. Other Approaches

Data exploration and predictive modelling represent the large majority of chemometric applications in the context of spectroscopy-based assessment of the quality of foodstuff. Nevertheless, there also exist other strategies which may provide relevant insights into the products or raw materials under investigation when the aim is their characterization and/or authentication. In this section, these techniques will be briefly illustrated.

6.1. Curve Resolution

In the previous sections, it was shown how a spectral data matrix $X$ can be decomposed according to a bilinear model and how this decomposition can be regarded as its projection onto a subspace of latent (abstract) directions. However, under normal experimental conditions, it is well known that spectroscopic data can also be approximated by a bilinear model whose elements are directly interpretable in chemical terms. Indeed, by assuming that $X$ contains the spectral intensities or absorbance values for the various analyzed samples at different wavelengths, and, without any loss of generality, that the optical path is unitary, Beer-Lambert-Bouger law can be expressed as:

$$X=C{S}^{\mathrm{T}}$$

(32)

where $C$ and $S$ are the arrays collecting the concentrations of the chromophores constituting the spectral mixture across all the specimens and their individual spectroscopic signatures (molar absorbivities), respectively [221]. Based on the relation described by Equation (31), the aim of curve resolution methods is to perform a data-driven unmixing of the matrix $X$, i.e., to estimate the number of constituents in the aforementioned mixture and obtain their pure concentration and spectral profiles from the information encoded in the recorded data. Among the possible approaches for accomplishing such a task, the most popular (due to its flexibility and to the possibility of being applied also to multi-set or multi-way data) is Multivariate Curve Resolution–Alternating Least Squares (MCR-ALS [222,223,224,225,226,227]). As the name suggests, MCR-ALS is an iterative algorithm based on the alternating least squares concept, which involves the following steps:

• The number of components (mixture constituents) is estimated (e.g., according to a priori knowledge of the systems under study or by SVD/PCA).
• A first guess of either $C$ or $S^{\mathrm{T}}$ is calculated by methodologies like Evolving Factor Analysis (EFA [228]) or SIMPLe-to-use Interactive Self-Modeling Analysis (SIMPLISMA [229]).
• $C$ and $S^{\mathrm{T}}$ are iteratively updated using alternating least squares under appropriate constraints (e.g., non-negativity of the values in $C$ and/or $S^{\mathrm{T}}$) as:

  $$\begin{array}{c}C=XS{\left(S^{\mathrm{T}}S\right)}^{-1}\\ S^{\mathrm{T}}={\left(C^{\mathrm{T}}C\right)}^{-1}{C}^{\mathrm{T}}X\end{array}$$

  (33)
• Step 3 is repeated until a certain convergence criterion is met.

It is evident for the algorithmic procedure outlined before that one of the main features of MCR-ALS is its “self-modeling” nature, i.e., in principle, it does not require any specific preliminary information about the data at hand: it is only needed that the bilinear model in Equation (31) is satisfied and that some very generic characteristics of the pure concentration or spectral profiles (e.g., their non-negativity) are known. If further information is also available, though, it can be actively exploited in the form of additional mathematical restrictions. Basically, it is the possibility of implementing in a rather straightforward way a large variety of these restrictions what has made and still makes MCR-ALS very popular compared to other curve resolution strategies. Indeed, they not only reduce the ambiguity associated with the final model (that is to say the range of possible $C$/$S^T$ sets returning the same fit when approximating $X$) but, at the same time, they render the resolved profiles more physico-chemically meaningful [230,231]. Examples of possible restrictions are unimodality (if the constrained profiles are expected to exhibit a unique global maximum), closure (if a mass or concentration balance among all or part of the mixture constituents holds), and selectivity (if some of the resolved species are absent in some experiments or do not absorb at specific wavelengths) [222,226,232,233,234].

Moreover, as already anticipated, MCR-ALS can easily be applied to multiple data matrices in a multi-set configuration. This provides a different way of analyzing multi-block data with respect to what discussed in Section 5, as here the information shared by the different investigated arrays is resorted to for improving the unmixing of the pure constituents. Concurrently, since some of these constituents can be present only in one or a few matrices and some others may be present in all, multi-set MCR-ALS is one of those multi-block techniques by which information about both common and distinctive components can be obtained [235].

Similarly to ICA, MCR-ALS has mainly been employed in order to retrieve the spectral contributions of particular compounds constituting complex food samples and to determine how their abundance/concentration varies during, e.g., a renneting process (see

Table 4. Examples of applications of curve resolution and analysis of variance-based techniques in the fields of spectroscopy-based food analysis. MCR-ALS and ASCA stand for Multivariate Curve Resolution-Alternating Least Squares and ANOVA-Simultaneous Component Analysis, respectively. The abbreviation ANOVA denotes the analysis of variance.

Aim	Method	Reference
Chemical characterization of milk lactic acid fermentation	MCR-ALS	[236]
Milk renneting characterization and monitoring	MCR-ALS	[237,238]
Chemical characterization of beer fermentation	MCR-ALS	[239]
Assessment of coconut oil purity/adulteration degree	MCR-ALS	[265]
Chemical characterization of chocolate samples	MCR-ALS	[266]
Egg-pasta characterization	ASCA	[76]
Coffee bean roasting monitoring	ASCA	[150]
Cheddar cheese ripening monitoring	ASCA	[264]
Eggplant chilling injury characterization	ASCA	[267]

). In [236,237,238], for instance, the dynamic evolution of various forms of milk (liquid, sol-gel, coagulated) over the progression of a lactic acid fermentation reaction was monitored and their individual FT-NIR fingerprints disentangled. In [239], a similar study was conducted on beer fermentation: here, maltose, maltotriose, fructose, sucrose, dextrins and ethanol spectral and concentration profiles were retrieved for enabling a better understanding of the various phenomena behind their biotransformation.

6.2. Analysis of Multivariate Designed Data

The rational design of the experiments to be conducted is a fundamental step in any scientific discipline and the only way in which it can be guaranteed that the information sought be present in the collected data. The design of experiments involves four steps: (i) identifying some critical variables (factors) which could affect one or more properties of interest (responses), (ii) defining their range of variability that could be worth inspecting, (iii) planning the trials to be performed (so as to be at the same time as informative and as low in number as possible), and (iv) interpreting the final results [240,241,242,243]. In this subsection, attention will be essentially focused on this last step, as it is the one where traditional statistical approaches commonly fail when instead of measuring a few (uncorrelated) responses for each experimental condition setting, a full spectroscopic signature is registered to characterize the concerned samples [244].

As summarized before, based on a rational design, experiments are conducted under different conditions, characterized by the fact that the factors whose effect is investigated are fixed (or controlled) at specific values (levels). In other words, each combination of factor levels is a so-called design point, at which the measurement of the responses is carried out. For instance, one could be interested in studying how the spectroscopic profile of a specific product varies as a function of temperature and pH, and then decide to run experimental trials only at three levels of temperature (25 °C, 40 °C and 55 °C) and at two levels of pH (5 and 9). In case all possible combinations of factor levels are explored, the design is called full-factorial. In the example sketched before, a full-factorial design would encompass six distinct experiments: 25 °C and pH = 5; 40 °C and pH = 5; 55 °C and pH = 5; 25 °C and pH = 9; 40 °C and pH = 9; 55 °C and pH = 9. At each condition setting, the desired responses should be measured; if the aim is a spectroscopic characterization, the response is highly multivariate and it is represented by the spectroscopic signature of the respective sample. Moreover, very often, the whole design is replicated at least twice to have an estimate of the variability not to be ascribed to the controlled factors.

Given a particular design and having collected the data corresponding to every one of its points, the successive step is to analyze such data to first evaluate whether any of the controlled factors exhibits a statistically significant effect on the responses and (if this is the case) to interpret and describe in greater detail how this effect occurs (for instance, which spectral regions change with the increase of the temperature and how). The traditional statistical tool used in this type of contexts is the analysis of variance—univariate (ANOVA [245,246,247,248]) or multivariate (MANOVA [58,241,244,249,250,251,252]), depending on the nature of the response(s) measured. However, when dealing with a high number of possibly very correlated descriptors (like spectral variables) and, therefore, with ill-conditioned data matrices, MANOVA suffers from the same drawbacks as MLR or LDA. To overcome these limitations, in the last 15 years, various approaches for the assessment of multivariate or megavariate data resulting from designed experiments have been proposed, all of them having in common the same initial partitioning of the experimental data array, which follows the classical ANOVA scheme based on the linear additivity of the effects. In particular, by assuming that the experimental design at hand involves two factors (generically named A and B), the data matrix $X$ made up of all the spectra collected at the different design points after mean centering (i.e., after subtraction of the mean spectrum calculated across the entire dataset) is partitioned as follows:

$$X=X_{\mathrm{A}}+X_{\mathrm{B}}+X_{\mathrm{AB}}+X_{\mathrm{res}}$$

(34)

where $X_{\mathrm{A}}$, $X_{\mathrm{B}}$, and $X_{\mathrm{AB}}$ account for the effect of the two factors and their interaction, respectively, while $X_{\mathrm{res}}$ contains the residual variance in $X$ not explained by the model. More specifically, $X_{\mathrm{A}}$, $X_{\mathrm{B}}$, and $X_{\mathrm{AB}}$ carry identical replicates of the mean spectra corresponding to the different factor or interaction levels. Starting from the common ground defined by Equation (33), the various methods developed for the analysis of multivariate designed data proceed along different directions. Among these techniques, the most popular is surely ANOVA–Simultaneous Component Analysis (ASCA [253,254]), due to its relative simplicity and remarkable interpretability properties. In ASCA, the effect of each design term (factor or interaction) is estimated through the sum of squares of the elements of the corresponding matrix in Equation (29). For example, the effect of factor A is quantified as:

$$SS{Q}_{\mathrm{A}}={{\displaystyle \sum }}_n{{\displaystyle \sum }}_j{x}_{n,j,\mathrm{A}}^2=X_{\mathrm{A}}{}^2$$

(35)

and its statistical significance is evaluated non-parametrically by comparing $SS{Q}_{\mathrm{A}}$ with its null-distribution retrieved by permutations [255,256]. The interpretation of the significant effects is finally carried out by applying PCA to the respective matrices.

Alternative techniques described in the literature are ANOVA-PCA (APCA [257,258]), ANOVA-Target Projection (ANOVA-TP [259,260,261]), ANOVA–Common Dimensions (AComDim [262]) and regularized MANOVA (rMANOVA [263]).

Considering their underlying model structure, the methodologies outlined in this subsection have been widely exploited to obtain insights into how certain technological factors affect food-related processes like ageing (see Table 4); ASCA, for instance, has often been coupled to MIR and/or NIR spectroscopy to determine and examine the influence of (i) drying temperature and time on dried egg-pasta manufacturing [76], (ii) varietal origin and roasting time on coffee bean spectral properties [150], and (iii) storage temperature and conditions on Cheddar cheese ripening [264]. In all these case studies, the application of ASCA has generally permitted to highlight in a rapid, non-destructive, relatively cheap and green fashion the occurrence of specific phenomena (e.g., proteo- and lypolisis) ongoing along with the evolution of the monitored process itself.

7. Additional Fundamental Aspects of Chemometric Modelling: Data Preprocessing and Validation

In the previous sections, the main chemometric approaches for tackling data exploration, calibration and classification were critically discussed and examples of their application were illustrated both for single- and multi-block cases. Nevertheless, data analysis is not limited to the model building phase only; rather, there exist at least two further steps (before and after the construction of a chemometric model itself) that play a fundamental role in determining the quality of the obtained outcomes and their reliability: preprocessing and validation.

7.1. Data Pre-Processing

Raw data are almost always not suited to be analyzed as such, i.e., directly in the form in which they are yielded by the instrument, and spectroscopic profiles in particular result from multiple contributions, and only a few of these contributions might be relevant for the problem under study. The combination of multiple sources of unwanted variability can have a significant impact both on the model performance and on its interpretation; for instance, bilinear approaches extract components according to some variance/covariance-based criteria which may lead to nuisance artifacts in the presence of a high amount of spurious/uninteresting information. Accordingly, data are preliminarily subjected to one or more transformations (preprocessing), so as to remove or reduce their so-called detrimental fraction [268]. When considering the possible ways of pretreating experimental data, it is necessary to point out that there are some strategies which are valid independently from the nature of the registered signal, while others are conceived for more specific types of instrumental responses. Within the first family of techniques, the two most commonly used operations are centering and scaling [269]. Centering aims at removing variable offsets; normally, mean-centering is adopted, i.e., the mean value computed across each column of a data matrix is subtracted from all the elements of that column. This makes it possible to discard what is shared by the entire set of samples and magnify the differences among individuals. On the other hand, scaling is frequently applied to guarantee that the scales of the different recorded variables are comparable as well as their a priori contribution to the chemometric model. Mathematically, it consists of dividing all the elements of each column of the data matrix by a constant term (usually its standard deviation), and it is highly recommended when the descriptors have been collected using distinct analytical platforms and, therefore, have been reported in different measurement units. It goes without saying that scaling is rarely performed on spectral profiles.

As mentioned before, there also exist several other preprocessing strategies which directly focus on the nature of the instrumental response itself and are designed to correct/remove unwanted sources of variability that are related to its specific characteristics, such as stochastic noise, non-linearities, baseline or wavelength shifts, or undesired variations due to particular chemical and physical phenomena (such as light scattering). The impact of stochastic noise can be reduced by the use of methods which filter high-frequency contributions out of the investigated signal. Fourier [270] or wavelet [271] transforms or Savitzy-Golay smoothing [272] are possible solutions to cope with such an issue. The Savitzky-Golay approach, in particular, is worth mentioning in greater detail, as it can also be exploited for signal or spectral differentiation: it fits successive subsets of adjacent variables with a low-order (usually second or third) polynomial and approximates the variable around which each subset is centered with the value resulting from the respective estimated function.

Differentiation methodologies make it possible to remove additive (first derivative) and multiplicative (second derivative) effects from the registered profiles and to deconvolve, at least partially, overlapping peaks. However, they exhibit the drawback of significantly decreasing the signal-to-noise ratio. This is why they often require a preliminary smoothing operation—Savitzky-Golay differentiation, for example, couples the smoothing procedure described before to the point-to-point calculation of the derivative of the interpolating polynomial.

The impact of light scattering on spectroscopic data, especially in the NIR domain, can also be regarded as a multiplicative effect and corrected by the use of second-order derivation. Alternatively, scatter correction can be achieved by techniques like Standard Normal Variate (SNV) [273], Multiplicative Scatter Correction (MSC) [274] and Extended Multiplicative Signal Correction (EMSC) [275], which enables the simultaneous removal of the effects of scattering, baseline and potential interferents (if their spectroscopic signatures are available).

Finally, for the correction of a non-constant baseline, detrending [273] and penalized asymmetric least squares regression [276] represent feasible options. The former corresponds to fitting a polynomial function to all the variables constituting the spectral signal, but in those cases in which the baseline does not account for a large portion of the overall data variance, severe artifacts may be produced. The latter is underlain by a more effective weighted fitting procedure that automatically identifies the variables most likely to be affected by the baseline contribution and weights them more for the calculation of the approximating curve.

7.2. Validation

Having built a chemometric model does not guarantee per se its quality and/or validity, as these two aspects are strictly influenced by user-dependent factors like the chosen analysis technique or the values of certain metaparameters [2]. To assess whether it can lead to solid conclusions, the outcomes/predictions resulting from it are reliable and generalizable and its interpretation is sound and meaningful, such a model needs, therefore, to be validated [277]. Operationally, validation consists in evaluating the model’s performance on a dataset (called a test set) different from the one used in the training or calibration phase (referred to as a training or calibration set, as also specified before), but reflecting as closely as possible the distribution of future individuals and accounting for all potential sources of expected variability (e.g., in terms of geographical origin, composition, manufacturing process, etc.)—it is fundamental not to use for validation the same data exploited for model building as this could lead to overoptimistic considerations (being data-driven, chemometric models are conceived so as to return the best performance on those data). Unfortunately, in many real-world scenarios, it is not possible to design additional sampling campaigns and assemble completely new test sets. In these circumstances, it is common practice to split all the recorded measurement observations into two different blocks, one for model building and the other for model validation. However, for the validation procedure to be effective, the sample splitting scheme should ensure that both training and test sets span a representative amount of the overall data variability: in this sense, one should avoid random partitioning and, instead, resort to targeted subset selection algorithms like Kennard-Stone [278], Duplex [279,280] or D-optimal design-based approaches [281].

An alternative is constituted by internal resampling strategies [282], like cross-validation [70], in which the original data are repeatedly divided into a training and a test set and the procedure is iterated either until each object has been left out at least once from the former or until a maximum number of computational cycles has been reached. Cross-validation is often exploited when small sample sizes are concerned but, due to its underlying principle, usually yields optimistically biased outputs, as training and test sets are never completely independent from one another. Nevertheless, cross-validation can be fruitfully adopted to estimate optimal values for model metaparameters like the number of latent variables in PLS regression or PLSDA classification.

8. Software

Nowadays, multiple software tools for chemometric modelling, both commercial and open-source, are available and easily accessible. A detailed list can be found in [1]. In particular, a collection of freely downloadable Matlab (The Mathworks Inc., Natick, MA, USA) functions for the implementation of most of the methods described in this review can be found at https://www.chem.uniroma1.it/romechemometrics/research/algorithms/.

9. Conclusions

This article offers a global overview of the chemometric approaches most commonly used in the field of spectroscopy-based food analysis and authentication. Three different scenarios were surveyed: data exploration, regression and classification. Basic and simple descriptions of the main multivariate techniques employed in such a domain along with a comprehensive outline of their most recent and (as far as the authors are concerned) interesting applications were provided. Data preprocessing- and model validation-related issues were also thoroughly covered.