A robust approach to model-based classification based on trimming and constraints

TL;DR

This paper introduces a robust model-based classification method that employs trimming and eigenvalue constraints to improve classification accuracy in contaminated datasets, especially when training data is limited.

Contribution

It presents a novel robust modification to the Model-Based Classification framework using trimming and eigenvalue constraints, enhancing noise handling and reliability.

Findings

Effective in handling noise in response and explanatory variables

Provides reliable classification with contaminated datasets

Demonstrates benefits through experiments on real and simulated data

Abstract

In a standard classification framework a set of trustworthy learning data are employed to build a decision rule, with the final aim of classifying unlabelled units belonging to the test set. Therefore, unreliable labelled observations, namely outliers and data with incorrect labels, can strongly undermine the classifier performance, especially if the training size is small. The present work introduces a robust modification to the Model-Based Classification framework, employing impartial trimming and constraints on the ratio between the maximum and the minimum eigenvalue of the group scatter matrices. The proposed method effectively handles noise presence in both response and exploratory variables, providing reliable classification even when dealing with contaminated datasets. A robust information criterion is proposed for model selection. Experiments on real and simulated data, artificially adulterated, are provided to underline the benefits of the proposed method.

Tables

Table 1. Table 1: Nomenclature, covariance structure and number of free parameters in 𝚺 1 , … , 𝚺 G subscript 𝚺 1 … subscript 𝚺 𝐺 \bm{\Sigma}_{1},\ldots,\bm{\Sigma}_{G} : γ 𝛾 \gamma denotes the number of parameters related to the orthogonal rotation and δ 𝛿 \delta the number of parameters related to the eigenvalues. The last column indicates whether the eigenvalue-ratio (ER) constraint is required.

Model	${𝚺}_g$	$\gamma$	$\delta$	ER
EII	$\lambda 𝑰$	-	$1$	Not required
VII	${\lambda }_g𝑰$	-	$G$	Required
EEI	$\lambda 𝑨$	-	$p$	Not required
VEI	${\lambda }_g𝑨$	-	$G+p-1$	Required
EVI	$\lambda {𝑨}_g$	-	$Gp-\left(G-1\right)$	Required
VVI	${\lambda }_g{𝑨}_g$	-	$Gp$	Required
EEE	$\lambda 𝑫𝑨{𝑫}^{{}^{\prime }}$	$p\left(p-1\right)/2$	$p$	Not required
VEE	${\lambda }_g𝑫𝑨{𝑫}^{{}^{\prime }}$	$p\left(p-1\right)/2$	$G+p-1$	Required
EVE	$\lambda 𝑫{𝑨}_g{𝑫}^{{}^{\prime }}$	$p\left(p-1\right)/2$	$Gp-\left(G-1\right)$	Required
EEV	$\lambda {𝑫}_g𝑨{𝑫}_g^{{}^{\prime }}$	$Gp\left(p-1\right)/2$	$p$	Not required
VVE	${\lambda }_g𝑫{𝑨}_g{𝑫}^{{}^{\prime }}$	$p\left(p-1\right)/2$	$Gp$	Required
VEV	${\lambda }_g{𝑫}_g𝑨{𝑫}_g^{{}^{\prime }}$	$Gp\left(p-1\right)/2$	$G+p-1$	Required
EVV	$\lambda {𝑫}_g{𝑨}_g{𝑫}_g^{{}^{\prime }}$	$Gp\left(p-1\right)/2$	$Gp-\left(G-1\right)$	Required
VVV	${\lambda }_g{𝑫}_g{𝑨}_g{𝑫}_g^{{}^{\prime }}$	$Gp\left(p-1\right)/2$	$Gp$	Required

Table 2. Table 2: Average misclassification errors on B = 1000 𝐵 1000 B=1000 runs, varying method and contamination rate η 𝜂 \eta . Standard errors are reported in parenthesis.

$\eta$	$0$	$0.05$	$0.10$	$0.15$	$0.20$	$0.25$
EDDA	0.009	0.031	0.053	0.079	0.099	0.112
	(0.005)	(0.026)	(0.043)	(0.051)	(0.054)	(0.05)
UPCLASS	0.008	0.041	0.091	0.142	0.166	0.186
	(0.004)	(0.056)	(0.088)	(0.088)	(0.08)	(0.067)
RMDA	0.009	0.045	0.049	0.052	0.07	0.08
	(0.005)	(0.072)	(0.063)	(0.057)	(0.068)	(0.073)
RLDA	0.027	0.027	0.026	0.026	0.026	0.067
	(0.009)	(0.009)	(0.009)	(0.008)	(0.009)	(0.037)
REDDA	0.01	0.01	0.01	0.009	0.024	0.042
	(0.005)	(0.005)	(0.005)	(0.005)	(0.014)	(0.014)
RUPCLASS	0.01	0.009	0.009	0.008	0.019	0.044
	(0.005)	(0.005)	(0.005)	(0.005)	(0.013)	(0.014)

Table 3. Table 3: Misclassification rates in the unlabelled set for different classification methods. Average values for 50 random splits in training and validation (three proportions are considered), standard deviations reported in parentheses.

Error Rate	EDDA	UPCLASS	RMDA	RLDA	SVM	AdaBoost
50% Tr - 50% Te	0.033	0.065	0.291	0.1	0.025	0.036
	(0.012)	(0.049)	(0.091)	(0.02)	(0.008)	(0.011)
25% Tr - 75% Te	0.078	0.112	0.303	0.12	0.048	0.042
	(0.025)	(0.028)	(0.08)	(0.04)	(0.021)	(0.012)
10% Tr - 90% Te	0.24	0.126	0.375	0.157	0.109	0.058
	(0.031)	(0.023)	(0.065)	(0.08)	(0.036)	(0.021)

Table 4. Table 4: Misclassification rates in the unlabelled set, % of wrongly labelled samples correctly trimmed in the labelled set and % of those correctly trimmed observations properly a-posteriori assigned to the Beet Sucrose group. Average values for 50 random splits in training and validation (three proportions are considered), standard deviations reported in parentheses.

		RSIMCA	REDDA	RUPCLASS
50% Tr - 50% Te	Error Rate	0.069	0.05	0.029
		(0.029)	(0.013)	(0.01)
	% Correctly Trimmed	1	0.977	1
		(0)	(0.075)	(0)
	% Correctly Assigned	1	1	1
		(0)	(0)	(0)
25% Tr - 75% Te	Error Rate	0.075	0.053	0.032
		(0.038)	(0.034)	(0.009)
	% Correctly Trimmed	1	0.88	0.96
		(0)	(0.25)	(0.145)
	% Correctly Assigned	1	0.963	1
		(0)	(0.162)	(0)
10% Tr - 90% Te	Error Rate	0.111	0.121	0.053
		(0.051)	(0.039)	(0.038)
	% Correctly Trimmed	0.99	0.47	0.73
		(0.071)	(0.238)	(0.381)
	% Correctly Assigned	0.99	0.72	0.84
		(0.019)	(0.071)	(0.37)

Equations

$$\begin{array}{r}\mathcal{G}\sim Mul{t}_G(1;{\tau }_1,\dots ,{\tau }_G)\\ \mathcal{X}|\mathcal{G}=g\sim {\mathcal{N}}_p({\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\end{array}$$

$$f({\mathbf{x}}_n,{\mathbf{l}}_n;\mathbf{\Theta })=\prod _{g=1}^G{\left[{\tau }_g\varphi ({\mathbf{x}}_n;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right]}^{l_{ng}}$$

$$z_{mg}=\mathbb{P}(\mathcal{G}=g|\mathcal{X}={\text{y}}_m)=\frac{{\tau }_g\varphi ({\mathbf{y}}_m;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)}{{\sum }_{j=1}^G{\tau }_j\varphi ({\mathbf{y}}_m;{\mathbf{\mu }}_j,{\mathbf{\Sigma }}_j)}.$$

$${\mathbf{\Sigma }}_g={\lambda }_g{\mathbf{D}}_g{\mathbf{A}}_g{\mathbf{D}}_g^{{}^{\prime }}$$

$$\begin{array}{rl}\ell (\mathbf{\tau },\mathbf{\mu },\mathbf{\Sigma }|\mathbf{X},\mathbf{Y},\mathbf{l})& =\sum _{n=1}^N\sum _{g=1}^G{l}_{ng}\log \left[{\tau }_g\varphi ({\mathbf{x}}_n;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right]+\\ & +\sum _{m=1}^M\log \left[\sum _{g=1}^G{\tau }_g\varphi ({\mathbf{y}}_m;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right]\end{array}$$

$$\begin{array}{rl}{\ell }_C(\mathbf{\tau },\mathbf{\mu },\mathbf{\Sigma }|\mathbf{X},\mathbf{Y},\mathbf{l},\mathbf{z})& =\sum _{n=1}^N\sum _{g=1}^G{l}_{ng}\log \left[{\tau }_g\varphi ({\mathbf{x}}_n;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right]+\\ & +\sum _{m=1}^M\sum _{g=1}^G{z}_{mg}\log \left[{\tau }_g\varphi ({\mathbf{y}}_m;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right]\end{array}$$

$$\begin{array}{rl}{\ell }_{trim}(\mathbf{\tau },\mathbf{\mu },\mathbf{\Sigma }|\mathbf{X},\mathbf{Y},\text{l})& =\sum _{n=1}^N\zeta ({\mathbf{x}}_n)\sum _{g=1}^G{\text{l}}_{ng}\log \left[{\tau }_g\varphi ({\mathbf{x}}_n;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right]+\\ & +\sum _{m=1}^M\phi ({\mathbf{y}}_m)\log \left[\sum _{g=1}^G{\tau }_g\varphi ({\mathbf{y}}_m;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right]\end{array}$$

$${\text{M}}_n/{\text{m}}_n\le c$$

$$\begin{array}{rl}\hat{\mathbf{\theta }}& =\{{\hat{\tau }}_1,\dots ,{\hat{\tau }}_G,{\hat{\mathbf{\mu }}}_1,\dots ,{\hat{\mathbf{\mu }}}_G,{\hat{\mathbf{\Sigma }}}_1,\dots ,{\hat{\mathbf{\Sigma }}}_G\}=\\ & =\{{\hat{\tau }}_1^{(0)},\dots ,{\hat{\tau }}_G^{(0)},{\hat{\mathbf{\mu }}}_1^{(0)},\dots ,{\hat{\mathbf{\mu }}}_G^{(0)},{\hat{\mathbf{\Sigma }}}_1^{(0)},\dots ,{\hat{\mathbf{\Sigma }}}_G^{(0)}\}\end{array}$$

$$f({\mathbf{x}}_n|l_{ng}=1;\hat{\mathbf{\theta }})=\varphi \left({\mathbf{x}}_n;{\hat{\mathbf{\mu }}}_g,{\hat{\mathbf{\Sigma }}}_g\right)g=1,\dots ,G.$$

$${\hat{\tau }}_g=\frac{{\sum }_{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}}{\lceil N(1-{\alpha }_l)\rceil }g=1,\dots ,G$$

$${\hat{\mathbf{\mu }}}_g=\frac{{\sum }_{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}{\mathbf{x}}_n}{{\sum }_{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}}g=1,\dots ,G.$$

$$D\left({\mathbf{x}}_n;{\hat{\mathbf{\theta }}}^{(k)}\right)=\prod _{g=1}^G{\left[\varphi \left({\mathbf{x}}_n;{\hat{\mathbf{\mu }}}_g^{(k)},{\hat{\mathbf{\Sigma }}}_g^{(k)}\right)\right]}^{l_{ng}}n=1,\dots ,N$$

$$D\left({\mathbf{y}}_m;{\hat{\mathbf{\theta }}}^{(k)}\right)=\sum _{g=1}^G{\hat{\tau }}_g^{(k)}\varphi \left({\mathbf{y}}_m;{\hat{\mathbf{\mu }}}_g^{(k)},{\hat{\mathbf{\Sigma }}}_g^{(k)}\right)m=1,\dots ,M.$$

$${\hat{z}}_{mg}^{(k+1)}=\frac{{\hat{\tau }}_g^{(k)}\varphi \left({\mathbf{y}}_m;{\hat{\mathbf{\mu }}}_g^{(k)},{\hat{\mathbf{\Sigma }}}_g^{(k)}\right)}{D\left({\mathbf{y}}_m;{\hat{\mathbf{\theta }}}^{(k)}\right)}g=1,\dots ,G;m=1,\dots ,M.$$

$${\hat{\tau }}_g^{(k+1)}=\frac{{\sum }_{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}+{\sum }_{m=1}^M\phi ({\mathbf{y}}_m){\hat{z}}_{mg}^{(k+1)}}{\lceil N(1-{\alpha }_l)\rceil +\lceil M(1-{\alpha }_u)\rceil }g=1,\dots ,G$$

$${\hat{\mathbf{\mu }}}_g^{(k+1)}=\frac{{\sum }_{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}{\mathbf{x}}_n+{\sum }_{m=1}^M\phi ({\mathbf{y}}_m){\hat{z}}_{mg}^{(k+1)}{\mathbf{y}}_m}{{\sum }_{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}+{\sum }_{m=1}^M\phi ({\mathbf{y}}_m){\hat{z}}_{mg}^{(k+1)}}g=1,\dots ,G.$$

$$a^{(k)}=\frac{{\ell }_{trim}^{(k+1)}-{\ell }_{trim}^{(k)}}{{\ell }_{trim}^{(k)}-{\ell }_{trim}^{(k-1)}}$$

$${\ell }_{{\infty }_{trim}}^{(k)}={\ell }_{trim}^{(k)}+\frac{1}{1-a^{(k)}}\left({\ell }_{trim}^{(k+1)}-{\ell }_{trim}^{(k)}\right).$$

$$RBIC=2{\ell }_{trim}(\hat{\mathbf{\tau }},\hat{\mathbf{\mu }},\hat{\mathbf{\Sigma }})-v_{XXX}^c\log \left(\lceil N(1-{\alpha }_l)\rceil +\lceil M(1-{\alpha }_u)\rceil \right)$$

$$v_{XXX}^c=Gp+G-1+\gamma +(\delta -1)\left(1-\frac{1}{c}\right)+1.$$

$$\mathbf{\tau }=(0.3,0.2,0.5{)}^{\prime },{\mathbf{\mu }}_1=(0,0{)}^{\prime },{\mathbf{\mu }}_2=(4,-4{)}^{\prime },{\mathbf{\mu }}_3=(0,8{)}^{\prime }$$

$${\mathbf{\Sigma }}_1=\left[\begin{array}{cc}1& 0.3\\ 0.3& 1\end{array}\right]{\mathbf{\Sigma }}_2=\left[\begin{array}{cc}1& -0.3\\ -0.3& 1\end{array}\right]{\mathbf{\Sigma }}_3=\left[\begin{array}{cc}6.71& 2.09\\ 2.09& 6.71\end{array}\right].$$

$${\mathbf{x}}_i=\mathbf{\mu }+\mathbf{\Lambda }{\mathbf{u}}_i+{\mathbf{e}}_i$$

$${\hat{\mathbf{u}}}_i={\hat{\mathbf{\Lambda }}}^{{}^{\prime }}{\left(\hat{\mathbf{\Lambda }}{\hat{\mathbf{\Lambda }}}^{{}^{\prime }}+\hat{\mathbf{\Psi }}\right)}^{-1}({\mathbf{x}}_i-\bar{\mathbf{x}})$$

$$\begin{array}{rl}{E}_{{\hat{\mathbf{\theta }}}^{(k)}}({\mathcal{Z}}_{mg}|{\mathbf{y}}_m)& =\mathbb{P}\left({\mathcal{Z}}_{mg}=1|{\mathbf{y}}_m;{\hat{\theta }}^{(k)}\right)=\\ & =\frac{\mathbb{P}\left({\mathbf{y}}_m|{\mathcal{Z}}_{mg}=1;{\hat{\theta }}^{(k)}\right)\mathbb{P}\left({\mathcal{Z}}_{mg}=1;{\hat{\theta }}^{(k)}\right)}{{\sum }_{j=1}^G\mathbb{P}\left({\mathbf{y}}_m|{\mathcal{Z}}_{mj}=1;{\hat{\theta }}^{(k)}\right)\mathbb{P}\left({\mathcal{Z}}_{mj}=1;{\hat{\theta }}^{(k)}\right)}=\\ & =\frac{{\hat{\tau }}_g^{(k)}\varphi \left({\mathbf{y}}_m;{\hat{\mathbf{\mu }}}_g^{(k)},{\hat{\mathbf{\Sigma }}}_g^{(k)}\right)}{{\sum }_{j=1}^G{\hat{\tau }}_j^{(k)}\varphi \left({\mathbf{y}}_m;{\hat{\mathbf{\mu }}}_j^{(k)},{\hat{\mathbf{\Sigma }}}_j^{(k)}\right)}=\\ & ={\hat{z}}_{mg}^{(k+1)}g=1,\dots ,G;m=1,\dots ,M.\end{array}$$

$$\begin{array}{rl}Q(\mathbf{\theta };{\hat{\mathbf{\theta }}}^{(k)})& =\sum _{n=1}^N\zeta ({\mathbf{x}}_n)\sum _{g=1}^G{l}_{ng}\log \left[{\tau }_g\varphi ({\mathbf{x}}_n;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right]+\\ & +\sum _{m=1}^M\phi ({\mathbf{y}}_m)\sum _{g=1}^G{\hat{z}}_{mg}\log \left[{\tau }_g\varphi ({\mathbf{y}}_m;{\mathbf{\mu }}_g,{\mathbf{\Sigma }}_g)\right].\end{array}$$

$$\mathcal{L}(\mathbf{\theta },\kappa )=Q(\mathbf{\theta };{\hat{\mathbf{\theta }}}^{(k)})-\kappa \left(\sum _{j=1}^G{\tau }_j-1\right)$$

$$\frac{\partial }{\partial {\tau }_g}\mathcal{L}(\mathbf{\theta },\kappa )=\frac{{\sum }_{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}}{{\tau }_g}+\frac{{\sum }_{m=1}^M\phi ({\mathbf{y}}_m){\hat{z}}_{mg}}{{\tau }_g}-\kappa$$

$$\sum _{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}+\sum _{m=1}^M\phi ({\mathbf{y}}_m){\hat{z}}_{mg}-\kappa {\tau }_g=0.$$

$${\hat{\tau }}_g^{(k+1)}=\frac{{\sum }_{n=1}^N\zeta ({\mathbf{x}}_n)l_{ng}+{\sum }_{m=1}^M\phi ({\mathbf{y}}_m){\hat{z}}_{mg}^{(k+1)}}{\lceil N(1-{\alpha }_l)\rceil +\lceil M(1-{\alpha }_u)\rceil }g=1,\dots ,G.$$

Full text

1 Introduction

In statistical learning, we define classification as the task of assigning group memberships to a set of unlabelled observations. Whenever a labelled sample (i.e., the training set) is available, the information contained in such dataset is exploited to classify the remaining unlabelled observations (i.e., the test set), either in a supervised or in a semi-supervised manner, depending whether the information contained in the test set are included in building the classifier (e.g. McNicholas, 2016). Either way, the presence of unreliable data points can be detrimental for the classification process, especially if the training size is small (Zhu and Wu, 2004).

Broadly speaking, noise is anything that obscures the relationship between the attributes and the class membership (Hickey, 1996). In a classification context, Wu (1995) distinguishes between two types of noise: attribute noise and class noise. The former is related to contamination in the exploratory variables, that is when observations present unusual values on their predictors; whereas the latter refers to samples whose associated labels are wrong. Zhu and Wu (2004) and the recent work of Prati et al. (2018) offer an extensive review on the topic and the methods that have been proposed in the literature to deal with attribute noise and class noise, respectively. Generally, three main approaches can be employed when building a classifier from a noisy dataset: cleaning the data, modeling the noise and using robust estimators of model parameters (Bouveyron and Girard, 2009).

The approach presented in this paper is based on a robust estimation of a Gaussian mixture model with parsimonious structure, to account for both attribute and label noise. Our conjecture is that the contaminated observations would be the least plausible units under the robustly estimated model: the corrupted subsample will be revealed by detecting those observations with the lowest contributions to the associated likelihood. Impartial trimming (Gordaliza, 1991a, b; Cuesta-Albertos et al., 1997) is employed for robustifying the parameter estimates, being a well established technique to treat mild and gross outliers in the clustering literature (García-Escudero et al., 2010) and here used, for the first time, to additionally account for label noise in a classification framework. A semi-supervised approach is developed, where information contained in both labelled and unlabelled samples is combined for improving the classifier performance and for defining a data-driven method to identify outlying observations possibly present in the test set.

The rest of the manuscript is organized as follows. A brief review on model-based discriminant analysis and classification is given in Section 2, Section 3 introduces the robust updating classification rules, covering the model formulation, inference aspects and model selection. Simulation studies to compare the method introduced in Section 3 with other popular classification methods are reported in Section 4. Finally, in Section 5 our proposal is employed in performing classification and adulteration detection in a food authenticity context, dealing with contaminated samples of Irish honey. Concluding notes and further research directions are outlined in Section 6. The proof of Proposition 1, details on the parameters values for simulation study II and efficient algorithms for enforcing the eigen-ratio constraint for different patterned models are deferred respectively to appendices A, B and C.

2 Model-Based Discriminant Analysis and Classification

In this Section we review the main concepts of supervised classification based on mixture models, with particular focus on Eigenvalue Decomposition Discriminant Analysis and its semi-supervised formulation, as introduced in Dean et al. (2006). This approach is the basis of the novel robust semi-supervised classifier introduced in Section 3.

2.1 Eigenvalue Decomposition Discriminant Analysis

Model-based discriminant analysis (McLachlan, 1992; Fraley and Raftery, 2002) is a probabilistic approach for supervised classification, in which a classifier is built from a complete set of learning observations $\{(\mathbf{x}_{1},\mathbf{l}_{1}),\ldots,(\mathbf{x}_{N},\mathbf{l}_{N})\}$; where $\mathbf{x}_{n}$ and $\mathbf{l}_{n}$, $n=1,\ldots,N$, are independent realizations of random vectors $\mathcal{X}\in\mathbb{R}^{p}$ and $\mathcal{G}\in\{1,\ldots,G\}$, respectively. That is, $\mathbf{x}_{n}$ denotes a $p$-variate observation and $\mathbf{l}_{n}$ its associated class label, such that $l_{ng}=1$ if observation $n$ belongs to group $g$ and [math] otherwise, $g=1,\ldots,G$. Considering a Gaussian framework, the probabilistic mechanism that is assumed to have generated the data is as follows:

[TABLE]

where $\mathcal{G}$ is multinomially distributed with $\tau_{g}$ probability of observing class $g$ and the conditional density of $\mathcal{X}$ given $\mathcal{G}$ is multivariate normal with mean vector $\bm{\mu}_{g}$ and variance covariance matrix $\bm{\Sigma}_{g}$. Therefore, the joint density of $(\mathbf{x}_{n},\mathbf{l}_{n})$ is given by:

[TABLE]

where $\phi(\cdot;\bm{\mu}_{g},\bm{\Sigma}_{g})$ denotes the multivariate normal density and $\bm{\Theta}$ represents the collection of parameters to be estimated, $\bm{\Theta}=\{\tau_{1},\ldots,\tau_{G},\bm{\mu}_{1},\ldots,\bm{\mu}_{G},\bm{\Sigma}_{1},\ldots,\bm{\Sigma}_{G}\}$. Discriminant analysis makes use of data with known labels to estimate model parameters for creating a classification rule. The trained classifier is subsequently employed for assigning a set of unlabelled observations $\mathbf{y}_{m}$, $m=1,\ldots,M$ to the class $g$ with the associated highest posterior probability:

[TABLE]

using the maximum a posteriori (MAP) rule. The afore-described framework is widely employed in classification tasks, thanks to its probabilistic formulation and well-established efficacy.

The number of parameters in the component variance covariance matrices grows quadratically with the dimension $p$. Thus, Bensmail and Celeux (1996) introduced a parsimonious parametrization proposing to enforce additional assumptions on the matrices structure, based on the eigen-decomposition of Banfield and Raftery (1993) and Celeux and Govaert (1995):

[TABLE]

where $\bm{D}_{g}$ is an orthogonal matrix of eigenvectors, $\bm{A}_{g}$ is a diagonal matrix such that $|\bm{A}_{g}|=1$ and $\lambda_{g}=|\bm{\Sigma}_{g}|^{1/p}$. This elements correspond respectively to the orientation, shape and volume (alternatively called scale) of the different Gaussian components. Allowing each parameter in (4) to be equal or different across groups, Bensmail and Celeux (1996) define a family of 14 patterned models, listed in Table 1. Such class of models is particularly flexible, as it includes very popular classification methods like Linear Discriminant Analysis and Quadratic Discriminant Analysis as special cases for the EEE and VVV models, respectively (Hastie and Tibshirani, 1996). Eigenvalue Decomposition Discriminant Analysis (EDDA) is implemented in the mclust R package (Fop et al., 2016).

2.2 Updating Classification Rules

Exploiting the assumption that the data generating process outlined in (1) is the same for both labelled and unlabelled observations, Dean et al. (2006) propose to include also the data whose memberships are unknown in the parameter estimation. That is, information about group structure that may be contained in both labelled and unlabelled samples is combined in order to improve the classifier performance, in a semi-supervised manner.

Under the framework defined in Section 2.1, and given the set of available information $\{(\mathbf{x}_{n},\mathbf{l}_{n})|n=1,\ldots,N\}\cup\{\mathbf{y}_{m}|m=1,\ldots,M\}$, the observed log-likelihood is

[TABLE]

in which both labelled and unlabelled samples are accounted for in the likelihood definition. Treating the (unknown) labels $z_{mg}$, $m=1,\ldots,M$, $g=1,\ldots,G$ as missing data and including them in the likelihood specification defines the so called complete-data log-likelihood:

[TABLE]

Maximum likelihood estimates for (5) are obtained through the EM algorithm (Dempster et al., 1977), iteratively computing the expected value for the unknown labels given the current set of parameter estimates (E-Step), and employing (6) to find maximum likelihood estimates for the unknown parameters (M-Step). The unlabelled data are then classified according to $\hat{z}_{mg}$, using the MAP. The updating classification rules was demonstrated to give improved classification performance over the classical model-based discriminant analysis in some food authenticity applications, particularly when the training size is small. An implementation of this can be found in the upclass R package (Russell et al., 2014).

3 Robust Updating Classification Rules

We introduce here a Robust modification to the Updating Classification Rule described in Section 2.2, with the final aim of developing a classifier whose performance is not affected by contaminated data, either in the form of label noise and outlying observations.

3.1 Model Formulation

The main idea of the proposed approach is to employ techniques originated in the branch of robust statistics to obtain a model-based classifier in which parameters are robustly estimated and outlying observations identified. We are interested in providing a method that jointly accounts for noise on response and exploratory variables, where the former might be present in the labelled set and the latter in both the labelled and unlabelled sets. We propose to modify the log-likelihood in (5) with a trimmed mixture log-likelihood (Neykov et al., 2007) and to employ impartial trimming and constraints on the covariance matrices for achieving both robust parameter estimation and identification of the unreliable sub-sample. Impartial trimming is enforced by considering the distinct structure of the likelihoods associated to the labelled and unlabelled sets, accounting for the possible label noise that might be present in the labelled sample (see Section 3.2 for details). Following the same notation introduced in Section 2.1, we aim at maximizing the trimmed observed data log-likelihood:

[TABLE]

where $\zeta(\cdot)$, $\varphi(\cdot)$ are 0-1 trimming indicator functions, that express whether observation $\mathbf{x}_{n}$ and $\mathbf{y}_{m}$ are trimmed off or not. A fixed fraction $\alpha_{l}$ and $\alpha_{u}$ of observations, belonging to the labelled and unlabelled set respectively, is unassigned by setting $\sum_{n=1}^{N}\zeta(\mathbf{x}_{n})=\lceil N(1-\alpha_{l})\rceil$ and $\sum_{m=1}^{M}\varphi(\mathbf{y}_{m})=\lceil M(1-\alpha_{u})\rceil$. In this way, the less plausible samples under the currently estimated model are tentatively trimmed out at each step of the iterations that leads to the final estimate. The labelled trimming level $\alpha_{l}$ and the unlabelled trimming level $\alpha_{u}$ account for possible adulteration in both sets. At the end of the iterations, a value of $\zeta(\mathbf{x}_{n})=0$ or $\varphi(\mathbf{y}_{m})=0$ corresponds to identify $\mathbf{x}_{n}$ or $\mathbf{y}_{m}$, respectively, as unreliable observations. Notice that impartial trimming automatically deals with both class noise and attribute noise, as observations that suffer from either noise structure will give low contribution to the associated likelihood.

Maximization of (7) is carried out via the EM algorithm, in which an appropriate Concentration Step (Rousseeuw and Driessen, 1999) is performed in both labelled and unlabelled sets at each iteration to enforce the impartial trimming. In addition, we protect the parameter estimation from spurious solutions, that may arise whenever one component of the mixture fits a random pattern in the data. We consider the eigenvalue-ratio restriction:

[TABLE]

where $\textbf{M}_{n}=\max_{g=1\ldots G}\max_{l=1\ldots p}d_{lg}$ and $\textbf{m}_{n}=\min_{g=1\ldots G}\min_{l=1\ldots p}d_{lg}$, with $d_{lg}$, $l=1,\ldots,p$ being the eigenvalues of the matrix $\bm{\Sigma}_{g}$ and $c\geq 1$ being a fixed constant (Ingrassia, 2004). Constraint (8) simultaneously controls differences between groups and departures from sphericity, by forcing the relative length of the axes of the equidensity ellipsoids, based on the multivariate normal distribution, to be smaller than $\sqrt{c}$ (García-Escudero et al., 2014). Notice that the constraint in (8) is still needed whenever either the shape or the volume is free to vary across components (García-Escudero et al., 2017), that is for all models in Table 1 that present “Required” entry in the ER column. The considered approach is the (semi)-supervised version of the methodology proposed in Dotto and Farcomeni (2019), which is framed in a completely unsupervised scenario. Feasible and computationally efficient algorithms for enforcing the eigen-ratio constraint for different patterned models are reported in the Appendix C.

3.2 Estimation Procedure

The EM algorithm for obtaining Maximum Trimmed Likelihood Estimates of the robust updating classification rule involves the following steps:

• •

Robust Initialization: set $k=0$. Employing only the labelled data, we obtain robust starting values for the mean vector $\bm{\mu}_{g}$ and covariance matrix $\bm{\Sigma}_{g}$ of the multivariate normal density for each group $g$, $g=1,\ldots,G$, employing the following procedure:

For each class $g$, draw a random $(p+1)$-subset $J_{g}$ and compute its empirical mean $\hat{\bm{\mu}}^{(0)}_{g}$ and variance covariance matrix $\hat{\bm{\Sigma}}^{(0)}_{g}$ according to the considered parsimonious structure. This procedure yields better initial subsets than drawing random $\lceil N(1-\alpha_{l})\rceil$-subsets directly, because the probability of drawing an outlier-free $(p+1)$-subset is much higher than that of drawing an outlier-free $\lceil N(1-\alpha_{l})\rceil$-subset (Hubert et al., 2018). 2. 2.

Set

[TABLE]

where $\hat{\tau}_{1}^{(0)}=\ldots=\hat{\tau}_{G}^{(0)}=1/G$. 3. 3.

For each $\mathbf{x}_{n}$, $n=1,\ldots,N$, compute the conditional density

[TABLE]

$\lfloor N\alpha_{l}\rfloor\%$ of the samples with lowest value of (9) are temporarily discarded as possible outliers, namely label noise and/or attribute noise. That is, $\zeta(\mathbf{x}_{n})=0$ for such observations. 4. 4.

The parameter estimates are updated, based on the non-discarded observations:

[TABLE]

Estimation of $\bm{\Sigma}_{g}$ depends on the considered patterned model, details are given in Bensmail and Celeux (1996). 5. 5.

Iterate $3-4$ until the $\lfloor N\alpha_{l}\rfloor$ discarded observations are exactly the same on two consecutive iterations, then stop (usually, $\leq 3$ iterations are required).

The procedure described in steps $1-5$ is performed nsamp times, and the parameter estimates $\hat{\bm{\theta}}^{R}$ that lead to the highest value of the objective function $\ell_{trim}(\hat{\bm{\tau}},\hat{\bm{\mu}},\hat{\bm{\Sigma}}|\mathbf{X},\mathbf{\text{l}})=\sum_{n=1}^{N}\zeta(\mathbf{x}_{n})\sum_{g=1}^{G}\text{l}_{ng}\log{\left[\hat{\tau}_{g}\phi(\mathbf{x}_{n};\hat{\bm{\mu}}_{g},\hat{\bm{\Sigma}}_{g})\right]}$, out of nsamp repetitions, are retained. The afore-described procedure stems from the ideas of the FastMCD algorithm of Rousseeuw and Driessen (1999), here adapted for dealing with parsimonious structures in the covariance matrices. Retaining $\hat{\bm{\theta}}^{R}$ as final estimates leads to a fully supervised robust model-based method, called REDDA hereafter (see Section 4.1.1). Then, if the selected patterned model allows for heteroscedastic $\bm{\Sigma}_{g}$ and (8) is not satisfied, constrained maximization is enforced, see Appendix C for details.

• •

EM Iterations: denote by $\hat{\bm{\theta}}^{(k)}=\{\hat{\tau}_{1}^{(k)},\ldots,\hat{\tau}_{G}^{(k)},\hat{\bm{\mu}}_{1}^{(k)},\ldots,\hat{\bm{\mu}}_{G}^{(k)},\hat{\bm{\Sigma}}^{(k)}_{1},\ldots,\hat{\bm{\Sigma}}_{G}^{(k)}\}$ the parameter estimates at the $k$-th iteration of the algorithm.

• –

Step 1 - Concentration: the trimming procedure is implemented by discarding the $\lfloor N\alpha_{l}\rfloor$ observations $\mathbf{x}_{n}$ with smaller values of

[TABLE]

and discarding the $\lfloor M\alpha_{u}\rfloor$ observations $\mathbf{y}_{m}$ with smaller values of

[TABLE]

• –

Step 2 - Expectation: for each non-trimmed observation $\mathbf{y}_{m}$ compute the posterior probabilities

[TABLE]

• –

Step 3 - Constrained Maximization: the parameter estimates are updated, based on the non-discarded observations and the current estimates for the unknown labels:

[TABLE]

Estimation of $\bm{\Sigma}_{g}$ depends on the considered patterned model and on the eigenvalues-ratio constraint. Details are given in Bensmail and Celeux (1996) and, if (8) is not satisfied, in Appendix C.

• –

Step 4 - Convergence of the EM algorithm: check for algorithm convergence (see Section 3.3). If convergence has not been reached, set $k=k+1$ and repeat steps 1-4.

Notice how the trimming step differs between the labelled and unlabelled observations. We implicitly assume that a label in the training set conveys a sound meaning about the presence of a class of objects. Therefore, in the labelled set, we opted for trimming the samples with lowest conditional density $f(\mathbf{x}_{n}|l_{ng}=1;\hat{\bm{\theta}}^{(k)})=\phi\left(\mathbf{x}_{n};\hat{\bm{\mu}}^{(k)}_{g},\hat{\bm{\Sigma}}^{(k)}_{g}\right)$. The alternative choice of considering the joint density $f(\mathbf{x}_{n},l_{ng};\hat{\bm{\theta}}^{(k)})=\prod_{g=1}^{G}\left[\tau_{g}\phi(\mathbf{x}_{n};\bm{\mu}_{g},\bm{\Sigma}_{g})\right]^{l_{ng}}$ is instead prone to trim off completely groups with small prior probability $\tau_{g}$ for large enough value of $\alpha_{l}$, and should be discarded. Note that with (12) we are both discriminating label noise (i.e., observations that are likely to belong to the mixture model but whose associated label is wrong) and outliers. In the unlabelled set, on the other hand, trimming is based on the marginal density $f(\mathbf{y}_{m};\hat{\bm{\theta}}^{(k)})=\sum_{g=1}^{G}\hat{\tau}^{(k)}_{g}\phi\left(\mathbf{y}_{m};\hat{\bm{\mu}}^{(k)}_{g},\hat{\bm{\Sigma}}^{(k)}_{g}\right)$, having no prior information on the group membership of the samples.

Once convergence is reached, the estimated values $\hat{z}_{mg}$ provide a classification for the unlabelled observations $\mathbf{y}_{m}$, assigning observation $m$ into group $g$ if $\hat{z}_{mg}>\hat{z}_{mg^{{}^{\prime}}}$ for all $g^{{}^{\prime}}\neq g$. Final values of $\zeta(\mathbf{x}_{n})=0$, and $\varphi(\mathbf{y}_{m})=0$, classify $\mathbf{x}_{n}$ and $\mathbf{y}_{m}$ respectively, as outlying observations.

The routines for estimating the robust updating classification rules have been written in R language (R Core Team, 2018): the source code is available at https://github.com/AndreaCappozzo/rupclass. The estimation procedure detailed in this Section implies the monotonicity of the algorithm, according to:

Proposition 1: If the values $\zeta(\mathbf{x}_{n})$, $\varphi(\mathbf{y}_{m})$, $n=1,\ldots,N$, $m=1,\ldots,M$ are kept fixed, the EM algorithm described in Section 3.2 implies $\ell_{trim}(\hat{\bm{\theta}}^{(k+1)}|\mathbf{X},\mathbf{Y},\mathbf{\text{l}})\geq\ell_{trim}(\hat{\bm{\theta}}^{(k)}|\mathbf{X},\mathbf{Y},\mathbf{\text{l}})$ at any $k$.

The proof is reported in Appendix A. Furthermore, our estimation procedure reduces possible incorrect modes of the optimization function (spurious maximizers) and offers a constructive way to obtain a maximizer $\hat{\theta}_{n}$ for the sample problem, that converges to the global maximizer for the population, see García-Escudero et al. (2008) and García-Escudero et al. (2015).

3.3 Convergence Criterion

We assess whether the EM algorithm has reached convergence evaluating at each iteration how close the trimmed log-likelihood is to its estimated asymptotic value, using the Aitken acceleration (Aitken, 1926):

[TABLE]

where $\ell_{trim}^{(k)}$ is the trimmed observed data log-likelihood from iteration $k$. The asymptotic estimate of the trimmed log-likelihood at iteration $k$ is given by (Bohning et al., 1994):

[TABLE]

The EM algorithm is considered to have converged when $|\ell_{\infty_{trim}}^{(k)}-\ell_{trim}^{(k)}|<\varepsilon$; a value of $\varepsilon=10^{-5}$ has been chosen for the experiments reported in the next Sessions.

3.4 Model Selection

A robust likelihood-based criterion is employed for choosing the best model among the 14 patterned covariance structures listed in Table 1 and a reasonable value for the constraint $c$ in (8):

[TABLE]

where $\ell_{trim}(\hat{\bm{\tau}},\hat{\bm{\mu}},\hat{\bm{\Sigma}})$ denotes the maximized trimmed observed data log-likelihood and $v_{XXX}^{c}$ a penalty term whose definition is:

[TABLE]

That is, $v_{XXX}^{c}$ depends on the total number of parameters to be estimated: $\gamma$ and $\delta$ for every $XXX$ patterned model are given in Table 1. It also accounts for the trimming levels and for the eigen-ratio constraint $c$, according to Cerioli et al. (2018). Note that, when $c\rightarrow+\infty$ and $\alpha_{l}=\alpha_{u}=0$, (19) is the Bayesian Information Criterion (Schwarz, 1978).

4 Simulation studies

In this Section, we present two simulated data experiments: Simulation Study I compares the performances of several model-based classification methods in a low dimensional setting when dealing with noisy data at different contamination rates; Simulation Study II considers a higher dimensional scenario in which the accuracy performance of some popular classification methods is assessed, at a fixed contamination rate. In both scenarios we consider a joint noise structure on response and exploratory variables.

4.1 Simulation Study I

4.1.1 Experimental Setup

We consider a data generating process given by a mixture of $G=3$ components of bivariate normal distributions, according to the following parameters:

[TABLE]

[TABLE]

$600$ observations were generated from the model, randomly assigning $N=200$ to the labelled set and $M=400$ to the unlabelled set. The labelled set was subsequently adulterated with contamination rate $\eta$ (ranging from [math] to $0.25$), wrongly assigning $\lceil\eta/2N\rceil$ of the third group units to the first class and adding $\lceil\eta/2N\rceil$ randomly labelled points generated from a Uniform distribution on the square with vertices $\left[(-20,-20),(-20,20),(20,-20),(20,20)\right]$. The contamination is therefore twofold, involving jointly label switching and outliers for a total of $\eta N$ adulterated labelled units. Examples of labelled datasets with different contamination rates are reported in Figure 1. Performances of 6 model-based classification methods are considered:

• •

EDDA: Eigenvalue Decomposition Discriminant Analysis (Bensmail and Celeux, 1996)

• •

UPCLASS: Updating Classification Rules (Dean et al., 2006)

• •

RMDA: Robust Mixture Discriminant Analysis (Bouveyron and Girard, 2009)

• •

RLDA: Robust Linear Discriminant Analysis (Hawkins and McLachlan, 1997)

• •

REDDA: Robust Eigenvalue Decomposition Discriminant Analysis. This is the supervised version of the model described in Section 3, where only the labelled observations are used for parameter estimation obtained via the robust initialization detailed in Section 3.2.

• •

RUPCLASS: Robust Updating Classification Rules. The semi-supervised method described in Section 3.

To make a fair performance comparison, a level of $\alpha_{l}=0.15$ (REDDA and RUPCLASS) and $\alpha_{u}=0.05$ (RUPCLASS) have been kept fixed throughout the simulation study. Nevertheless, exploratory tools such as Density-Based Silhouette plot (Menardi, 2011) and trimmed likelihood curves (García-Escudero et al., 2011) could be employed to validate and assess the choice of $\alpha_{l}$ and $\alpha_{u}$. A more automatic approach, like the one introduced in Dotto et al. (2018), could also be adapted to our framework. This, however, goes beyond the scope of the present manuscript, it will nonetheless be addressed in the future. A value of $c=20$ was selected for the eigenvalue-ratio restriction in (8). Simulation study results are presented in the following subsections.

4.1.2 Classification Performance

Average misclassification errors for the different methods and for varying contamination rates are reported in Table 2 and in Figure 2. The error rate is computed on the unlabelled dataset and averaged over the $B=1000$ simulations. As expected, the misclassification error is fairly equal to all methods when there is no contamination rate, with the only exception being RLDA: this is due to the implicit model assumption that $\bm{\Sigma}_{1}=\bm{\Sigma}_{2}=\bm{\Sigma}_{3}$, which is not the case in our simulated scenario. As the contamination rate increases, so does the error rate for the non-robust methods (EDDA and UPCLASS), whereas for RLDA and RMDA it has a lower increment rate. Nevertheless, such methods fail to jointly cope with both sources of adulteration, namely class and attribute noise. Our proposals REDDA and RUPCLASS, thanks to the trimming step enforced in the estimation process, have always higher correct classification rates, on average, at any adulteration level. Notice that, to compare results of robust and non-robust methods, also the trimmed observations were classified a-posteriori according to the Bayes rule, assigning them to the component $g$ having greater value of $\hat{\tau}_{g}\phi(\mathbf{y}_{m};\hat{\bm{\mu}}_{g},\hat{\bm{\Sigma}}_{g})$.

On average, the robust semi-supervised approach performs better than the supervised counterpart, due to the information incorporated from genuine unlabelled data in the estimation process. Interestingly, the same behavior is not reflected in the non-robust counterparts, where the detrimental effect of contaminated labelled units magnifies the bias of the UPCLASS method. Therefore, robust solutions are even more paramount when a semi-supervised approach is considered.

4.1.3 Parameter Estimation

Figure 3 reports the box plots of the simulated estimation error over $B=1000$ Monte Carlo repetitions for the parameters of the mixture model, computing Euclidean norms for the proportion vector $\bm{\tau}$, the mean vectors $\bm{\mu}_{g}$ and covariance matrices $\bm{\Sigma}_{g}$, $g=1,\ldots,3$. The estimated values for the mixing proportion are mildly affected when increasing contamination is considered; conversely, the estimation of $\bm{\mu}_{2}$ is on average heavily influenced by the adulterating process, and also the robust methods fail to estimate it correctly as soon as the contamination rate $\eta$ is larger than the trimming level $\alpha_{l}=0.15$. Clearly, the estimation of the variance covariance matrices is as well badly affected in most extreme scenarios, where their entries are inflated in order to accommodate more and more bad points. Our robust proposals are less affected by the harmful effect of adding anomalous observations, also in the most adulterated scenario.

4.2 Simulation Study II

4.2.1 Experimental Setup

We consider here a simulating model with a larger number of features ($d=10$), where the data generating process is given by a mixture of $G=4$ components of a multivariate t-distribution with $\nu=6$ degrees of freedom. More details on the parameter values are contained in Appendix B. $1000$ observations were generated from the model, randomly assigning $N=250$ to the labelled set and $M=750$ to the unlabelled set. The training set was subsequently adulterated wrongly labelling $10$ units and adding $15$ randomly labelled outlying points, uniformly generated in the $d$-dimensional hypercube over $[10,15]^{10}$. We therefore consider a scenario in which $10\%$ of the learning units are contaminated, via both label and attribute noise.

Together with the model-based methods previously described in Section 4.1.1, we included in the performance evaluation widely used classification techniques that, even though not engineered to achieve robustness, are noise tolerant. Particularly, the ensemble learner AdaBoost (Freund and Schapire, 1997) and the kernel method Support Vector Machine (Cortes and Vapnik, 1995) were added to the comparison. Furthermore, the robust adaptation of the SIMCA method for high-dimensional classification (Vanden Branden and Hubert, 2005) was also considered. The classification performance of the afore-described techniques are tested against the proposed methodologies, under different combinations of $(c,\alpha_{l},\alpha_{u})$: accuracy results are reported in the next Section.

4.2.2 Classification Performance

Boxplots of the misclassification errors for the considered methods are reported in Figure 4. The error rate is computed on the $M=750$ units of the test set, under $B=1000$ repetitions of the generating process and subsequent adulteration scheme described in Section 4.2.1. As it was already apparent from the previous simulation study, accuracy for non-robust methods is badly affected by the contamination present in the learning set. Even though not specifically designed for dealing with adulterated datasets, SVM and AdaBoost perform better than the non-robust model-based approaches, thanks to their non-parametric nature and flexibility. As expected, the best classification accuracy are obtained by the robust methodologies, namely RLDA and our proposals REDDA and RUPCLASS. We also check the sensitivity of our techniques comparing different combinations of $(c,\alpha_{l},\alpha_{u})$. As it is easily visible in the boxplots, setting a smaller than needed labelled trimming level $\alpha_{l}$ leads to a loss in prediction accuracy, as a portion of adulterated units still affects the learning phase. Once the corrupted observations are correctly trimmed (i.e., $\alpha_{l}\text{ is set}\geq 0.1$), accuracy seems to remain stable with little influence induced by the choice of $c$ and $\alpha_{u}$, with only a slight preference for the semi-supervised RUPCLASS over its supervised version REDDA. This shows that setting a higher value of $\alpha_{l}$ is less detrimental than underestimating it, and that the impartial trimming almost exactly identifies the corrupted units when $\alpha_{l}=0.1$, that is the true adulteration proportion. The bad performance of RSIMCA is only due to the simulating process: given the fact that data truly lie on a 10-dimensional space, performing (robust) dimensional reduction prior to classification evidently leads to a concealment in the grouping structure.

The proposed methodologies were shown to be capable of dealing with data whose distribution is not exactly Gaussian, but where an effective robust decision rule can be built employing Gaussian mixture models.

5 Application to Midinfrared Spectroscopy of Irish Honey

The semi-supervised method introduced in Section 3 is employed in performing adulteration detection and classification in a food authenticity context: we consider the task of discriminating between pure and adulterated Irish Honey, where the training set itself contains unreliable samples.

5.1 Honey Samples

Honey is defined as “the natural sweet substance, produced by honeybees from the nectar of plants or from secretions of living parts of plants, or excretions of plant-sucking insects on the living parts of plants, which the bees collect, transform by combining with specific substances of their own, deposit, dehydrate, store and leave in honeycombs to ripen and mature” (Alimentarius, 2001). Being a relatively expensive commodity to produce and extremely variable in nature, honey is prone to adulteration for economic gain: in 2015 the European Commission organized an EU coordinated control plan to assess the prevalence on the market of honey adulterated with sugars and honeys mislabelled with regard to their botanical source or geographical origin. It is therefore of prime interest to employ robust analytical methods to protect food quality and uncover its illegal adulteration.

We consider here a dataset of midinfrared spectroscopic measurements of 530 Irish honey samples. Midinfrared spectroscopy is a fast, non-invasive method for examining substances that does not require any sample preparation, it is therefore an effective procedure for collecting data to be subsequently used in food authenticity studies (Downey, 1996). The spectra measurements lie in the wavelength range of $3700~{}\mbox{nm}$ and $13600~{}\mbox{nm}$, recorded at intervals of $35~{}\mbox{nm}$, with a total of 285 absorbance values. The dataset contains 290 Pure Honey observations, while the rest of the samples are honey diluted with adulterant solutions: 120 with Dextrose Syrup and 120 with Beet Sucrose, respectively. Kelly et al. (2006) gives a thorough explanation of the adulteration process. The aim of the study is to discriminate pure honey from the adulterated samples, when varying sample size of the labelled set whilst including a percentage of wrongly labelled units. Such a scenario is plausible to be encountered in real situations, since in a context in which the final purpose is to detect potential adulterated samples it may happen that the learning data is itself not fully reliable. An example of the data structure is reported in Figure 5.

5.2 Robust Dimensional Reduction

Prior to perform classification and adulteration detection, a preprocessing step is needed due to the high-dimensional nature of the considered dataset ($p=285$ variables). To do so, we robustly estimate a factor analysis model, retaining a set $d$ of factors, $d\ll p$, to be subsequently employed with the Robust Updating Classification Rules. Formally, for each Honey sample $\mathbf{x}_{i}$, we postulate a factor model of the form:

[TABLE]

where $\bm{\mu}$ is a $p\times 1$ mean vector, $\bm{\Lambda}$ is a $p\times d$ matrix of factor loadings, $\bm{u}_{i}$ are the unobserved factors, assumed to be realizations of a $d$-variate standard normal and the errors $\bm{e}_{i}$ are independent realizations of $\mathcal{N}(\mathbf{0},\bm{\Psi})$, with $\bm{\Psi}$ a $p\times p$ diagonal matrix. In such a way, the observed variables are assumed independent given the factors. For a general review on factor analysis, see for example Chapter 9 in (Mardia et al., 1979). Parameters in (21) are estimated employing a robust procedure based on trimming and constraints (García-Escudero et al., 2016), yielding dimensionality reduction at the same time. Given the robustly estimated parameters, the latent traits are computed using the regression method (Thomson, 1939):

[TABLE]

The estimated factors scores $\hat{\bm{u}}_{i}$ will be used for the classification task reported in the upcoming Section. For the considered dataset, after a graphical exploration of Cattell’s scree plot for the correlation matrix robustly estimated via MCD (Rousseeuw and Driessen, 1999), reported in Figure 6, we deem sufficient to set the number $d$ of latent factors equal to $10$. Parameters were estimated setting a trimming level $\alpha=0.1$ and $c_{noise}=1000$.

5.3 Classification Performance

SVM, AdaBoost and RSIMCA are designed to optimally perform in a high dimensional setting. Therefore, to respect the specificity of each family of methodologies, we directly applied SVM, AdaBoost and RSIMCA on the whole spectra. For EDDA, UPCLASS, RMDA, RLDA, REDDA and RUPCLASS we preprocessed the data with the dimension reduction method described in Section 5.2. To discriminate between pure and adulterated honey samples, we divided the available data into a training (labelled) sample and a validation (unlabelled) sample. We investigated the effect of having different sample sizes in the labelled set, both in terms of classification accuracy and adulteration detection. Particularly, 3 proportions have been considered: $50\%$ - $50\%$ , $25\%$ - $75\%$ and $10\%$ - $90\%$ for splitting data into training and validation set, respectively, within each group. For each split, $10\%$ of the Beet Sucrose adulterated samples were incorrectly labelled as Pure Honey in the training set, adding class noise in the discrimination task. The trimming levels $\alpha_{l}$ and $\alpha_{u}$ were set equal to $0.12$ and $0.05$, respectively. Table 3 and 4 summarize the accuracy results employing different classification approaches under the described scenarios. Careful investigation has been dedicated to measuring the ability of the robust methodologies in correctly determining (i.e., trimming) the $10\%$ of incorrectly labelled samples, that is, units adulterated with Beet Sucrose and erroneously labelled as Pure Honey: such information, only relevant for RSIMCA, REDDA and RUPCLASS models, is reported in Table 4. % Correctly Trimmed indicates the class noise percentage correctly detected by the impartial trimming. For the recognized class noise, % Correctly Assigned indicates the percentage of units properly a-posteriori assigned to the Beet Sucrose group. RSIMCA performs remarkably well in identifying the adulterated units, even though the classification accuracy is lower than the one obtained employing RUPCLASS model. As expected, the semi-supervised approach performs much better in terms of classification rate when the labelled sample size is small. Comparing the error rate of the robust techniques with the other methods in Table 3 we notice how powerful classifiers like SVM and AdaBoost work well also in dealing with adulterated datasets: SVM error rate in the 50% Tr - 50% Te is on average lower than the one obtained with RUPCLASS. However, when the labelled sample size decreases a semi-supervised approach is preferable: RUPCLASS reports the lowest error rate for both 25% Tr - 75% Te and 10% Tr - 90% Te scenarios. VEV and VVV models have been almost always chosen: model selection was performed through the Robust criteria defined in Section 3.4.

Results in Table 4 show that the proposed methodology is effective not only for accurately robustifying the parameter estimates, but also for efficiently detecting observations affected by class noise, firstly by trimming and subsequently by correctly assigning them: a critical information that cannot be obtained with standard classification methods like SVM and AdaBoost.

6 Concluding Remarks

In this paper we have proposed a robust modification to a family of semi-supervised patterned models, for performing classification in presence of both class and attribute noise.

We have shown that our methodology effectively addresses the issues generated by these two noise types, by identifying wrongly labelled units (noise in the response variable) and corrupted attributes in units (noise in the explanatory variables). Robust parameter estimates can therefore be obtained by excluding the noisy observations from the estimation procedure, both in the training set, and in the test set. Our proposal has been based on incorporating impartial trimming and eigenvalue-ratio constraints in previous semi-supervised methods. We have adapted the trimming procedure to the two different frameworks, i.e., for the labelled units and the unlabelled ones. After completing the robust estimation process, trimmed observations can be classified as well, by the usual Bayes rule. This final step allows the researcher to detect whether one observation is indeed extreme in terms of its attributes or it has been wrongly assigned to a different class. Such feature seems particularly desirable in food authenticity applications, where, due to imprecise readings and fraudulent units, it is likely to have label noise also within the labelled set. Some simulations, and a study on real data from pure and adulterated Honey samples, have shown the effectiveness of our proposal.

As an open point for further research, an automatic procedure for selecting reasonable values for the labelled and unlabelled trimming levels, along the lines of Dotto et al. (2018), is under study. Additionally, a robust wrapper variable selection for dealing with high-dimensional problems could be useful for further enhancing the discriminating power of the proposed methodology.

Acknowledgements

The authors are very grateful to Agustin Mayo-Iscar and Luis Angel García Escudero for both stimulating discussion and advices on how to enforce the eigenvalue-ratio constraints under the different patterned models. Andrea Cappozzo deeply thanks Michael Fop for his endless patience and guidance in helping him with methodological and computational issues encountered during the draft of the present manuscript. Brendan Murphy’s work is supported by the Science Foundation Ireland Insight Research Centre (SFI/12/RC/2289_P2)

Appendix A

Proof of Proposition 1: Considering the random variable $\mathcal{Z}_{mg}$ corresponding to $z_{mg}$, the E-step on the $(k+1)$th iteration requires the calculation of the conditional expectation of $\mathcal{Z}_{mg}$ given $\mathbf{y}_{m}$:

[TABLE]

Therefore, the Q function, to be maximized with respect to $\bm{\theta}$ in the M-step, is given by

[TABLE]

The maximization of (24) according to the mixture proportion $\tau_{g}$, $\sum_{j=1}^{G}\tau_{j}=1$ is solved considering the Lagrangian $\mathcal{L}(\bm{\theta},\kappa)$:

[TABLE]

with $\kappa$ the Lagrangian coefficient. The partial derivative of (25) with respect to $\tau_{g}$ has the form:

[TABLE]

and setting (26) equal to [math] for all $g=1,\ldots,G$ we obtain:

[TABLE]

Summing (27) over $g$, $g=1,\ldots,G$, provides the value of $\kappa=\lceil N(1-\alpha_{l})\rceil+M(1-\alpha_{u})\rceil$ and substituting it in the previous expression yields the ML estimate for $\tau_{g}$:

[TABLE]

The partial derivative of (24) with respect to the mean vector $\bm{\mu}_{g}$ reads:

[TABLE]

Equating (29) to [math] and rearranging terms provides the ML estimate of $\bm{\mu}_{g}$:

[TABLE]

Discarding quantities that do not depend on $\bm{\Sigma}_{g}$, we can rewrite (24) as follows:

[TABLE]

where $\bm{W}_{g}^{X}=\sum_{n=1}^{N}\zeta(\mathbf{x}_{n})l_{ng}\left[\left(\mathbf{x}_{n}-\bm{\mu}_{g}\right)\left(\mathbf{x}_{n}-\bm{\mu}_{g}\right)^{\prime}\right]$ and $\bm{W}_{g}^{Y}=\sum_{m=1}^{M}\varphi(\mathbf{y}_{m})\hat{z}_{mg}\left[\left(\mathbf{y}_{m}-\bm{\mu}_{g}\right)\left(\mathbf{y}_{m}-\bm{\mu}_{g}\right)^{\prime}\right]$. Finally, considering the eigenvalue decomposition $\bm{\Sigma}_{g}=\lambda_{g}\bm{D}_{g}\bm{A}_{g}\bm{D}^{{}^{\prime}}_{g}$, (31) simplifies to:

[TABLE]

The partial derivative of (32) with respect to $\left(\lambda_{g},\bm{A}_{g},\bm{D}_{g}\right)$ depends on the considered patterned structure: for a thorough derivation the reader is referred to Bensmail and Celeux (1996). If (8) is not satisfied, the constraints are enforced as detailed in Appendix C. Lastly, notice that in performing the concentration step the optimal observations of both training and test sets are retained, i.e. the ones with the highest contribution to the objective function.

The afore-described procedure falls within the structure of a general EM algorithm, for which the likelihood function does not decrease after an EM iteration, as shown in Dempster et al. (1977) and reported in page 78 of McLachlan and Krishnan (2008).

$\square$

Appendix B

This appendix details the structure of the Simulation Study in Section 4.2.1. We consider a data generating process given by a mixture of $G=4$ components of multivariate t-distributions (McLachlan and Peel, 1998; Peel and McLachlan, 2000), according to the following parameters:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A generalized pairs plot of contaminated labelled units under the afore-described Simulation Setup is reported in Figure 7.

Appendix C

This final Section presents feasible and computationally efficient algorithms for enforcing the eigenvalue-ratio constraint according to the different patterned models in Table 1. At the $k-$th iteration of the M step, the goal is to update the estimates for the variance-covariance matrices $\hat{\bm{\Sigma}}_{g}^{(k+1)}=\hat{\lambda}_{g}^{(k+1)}\hat{\bm{D}}_{g}^{(k+1)}\hat{\bm{A}}_{g}^{(k+1)}\hat{\bm{D}}^{{}^{\prime}(k+1)}_{g}$, $g=1,\ldots,G$ such that,

[TABLE]

where $\hat{a}_{lg}^{(k+1)}$ indicates the diagonal entries of matrix $\hat{\bm{A}}_{g}^{(k+1)}$. Denote with $\hat{\Sigma}_{g}^{U}=\hat{\lambda}_{g}^{U}\hat{\bm{D}}_{g}^{U}\hat{\bm{A}}_{g}^{U}\hat{\bm{D}}_{g}^{{}^{\prime}U}$ the estimates for the variance covariance matrices obtained following Bensmail and Celeux (1996) without enforcing the eigenvalues-ratio restriction in (33). Lastly, denote with $\hat{\bm{\Delta}}^{U}_{g}=\hat{\lambda}_{g}^{U}\hat{\bm{A}}_{g}^{U}$ the matrix of eigenvalues for $\hat{\bm{\Sigma}}_{g}^{U}$, with diagonal entries $\hat{d}_{lg}^{U}=\hat{\lambda}_{g}^{U}\hat{a}_{lg}^{U}$, $l=1,\ldots,p$.

Constrained maximization for VII, VVI and VVV models

Compute $\bm{\Delta}_{g}$ applying the optimal truncation operator defined in Fritz et al. (2013) to $\left\{\hat{\bm{\Delta}}^{U}_{1},\ldots,\hat{\bm{\Delta}}^{U}_{G}\right\}$, under condition (33) 2. 2.

Set $\hat{\lambda}_{g}^{(k+1)}=|\bm{\Delta}_{g}|^{1/p}$, $\hat{\bm{A}}_{g}^{(k+1)}=\frac{1}{\hat{\lambda}_{g}^{(k+1)}}\bm{\Delta}_{g}$, $\hat{\bm{D}}_{g}^{(k+1)}=\hat{\bm{D}}_{g}^{U}$

Constrained maximization for VVE model

Compute $\bm{\Delta}_{g}$ applying the optimal truncation operator defined in Fritz et al. (2013) to $\left\{\hat{\bm{\Delta}}^{U}_{1},\ldots,\hat{\bm{\Delta}}^{U}_{G}\right\}$, under condition (33) 2. 2.

Given $\bm{\Delta}_{g}$, compute the common principal components $\bm{D}$ via, for example, a majorization-minimization (MM) algorithm (Browne and McNicholas, 2014) 3. 3.

Set $\hat{\lambda}_{g}^{(k+1)}=|\bm{\Delta}_{g}|^{1/p}$, $\hat{\bm{A}}_{g}^{(k+1)}=\frac{1}{\hat{\lambda}_{g}^{(k+1)}}\bm{\Delta}_{g}$, $\hat{\bm{D}}_{g}^{(k+1)}=\bm{D}$

Constrained maximization for EVI, EVV models

Compute $\bm{\Delta}_{g}$ applying the optimal truncation operator defined in Fritz et al. (2013) to $\left\{\hat{\bm{\Delta}}^{U}_{1},\ldots,\hat{\bm{\Delta}}^{U}_{G}\right\}$, under condition (33) 2. 2.

Compute $\bm{\Delta}^{\star}_{g}$ constraining $\bm{\Delta}_{g}$ such that $\bm{\Delta}^{\star}_{g}=\lambda^{\star}\bm{A}_{g}^{\star}$. That is, constraining $|\bm{\Delta}^{\star}_{g}|$ to be equal across groups (Maronna and Jacovkis, 1974; Gallegos, 2002). Details are given in Section 3.2 of Fritz et al. (2012) 3. 3.

Iterate $1-2$ until (33) is satisfied 4. 4.

Set $\hat{\lambda}_{g}^{(k+1)}=\lambda^{\star}$, $\hat{\bm{A}}_{g}^{(k+1)}=\bm{A}_{g}^{\star}$, $\hat{\bm{D}}_{g}^{(k+1)}=\hat{\bm{D}}_{g}^{U}$

Constrained maximization for EVE model

Compute $\bm{\Delta}_{g}$ applying the optimal truncation operator defined in Fritz et al. (2013) to $\left\{\hat{\bm{\Delta}}^{U}_{1},\ldots,\hat{\bm{\Delta}}^{U}_{G}\right\}$, under condition (33) 2. 2.

Compute $\bm{\Delta}^{\star}_{g}$ constraining $\bm{\Delta}_{g}$ such that $\bm{\Delta}^{\star}_{g}=\lambda^{\star}\bm{A}_{g}^{\star}$. Details are given in Section 3.2 of Fritz et al. (2012) 3. 3.

Iterate $1-2$ until (33) is satisfied 4. 4.

Given $\bm{A}^{\star}_{g}$, compute the common principal components $\bm{D}$ via, for example, a majorization-minimization (MM) algorithm (Browne and McNicholas, 2014) 5. 5.

Set $\hat{\lambda}_{g}^{(k+1)}=\lambda^{\star}$, $\hat{\bm{A}}_{g}^{(k+1)}=\bm{A}_{g}^{\star}$, $\hat{\bm{D}}_{g}^{(k+1)}=\bm{D}$

Constrained maximization for VEI, VEV models

Set $\bm{\Delta}_{g}=\hat{\bm{\Delta}}^{U}_{g}$ 2. 2.

Set $\lambda_{g}^{\star}=\hat{\lambda}_{g}^{U}$, $g=1,\ldots,G$ 3. 3.

Compute $\bm{\Delta}^{\star}_{g}$ applying the optimal truncation operator defined in Fritz et al. (2013) to $\left\{\bm{\Delta}_{1},\ldots,\bm{\Delta}_{G}\right\}$, under condition (33) 4. 4.

Compute $\bm{A}^{\star}=\left.\sum_{g=1}^{G}\frac{1}{\lambda_{g}^{\star}}\bm{\Delta}^{\star}_{g}\middle/\left|\sum_{g=1}^{G}\frac{1}{\lambda_{g}^{\star}}\bm{\Delta}^{\star}_{g}\right|^{1/p}\right.$ 5. 5.

Compute $\lambda_{g}^{\star}=\frac{1}{p}tr\left(\bm{\Delta}^{\star}_{g}{\bm{A}^{\star}}^{-1}\right)$ 6. 6.

Set $\bm{\Delta}_{g}=\lambda_{g}^{\star}\bm{A}^{\star}$ 7. 7.

Iterate $3-6$ until (33) is satisfied 8. 8.

Set $\hat{\lambda}_{g}^{(k+1)}=\lambda^{\star}_{g}$, $\hat{\bm{A}}_{g}^{(k+1)}=\bm{A}^{\star}$, $\hat{\bm{D}}_{g}^{(k+1)}=\hat{\bm{D}}_{g}^{U}$

Constrained maximization for VEE model

Set $\bm{K}_{g}=\hat{\bm{\Sigma}}^{U}_{g}$ 2. 2.

Set $\lambda_{g}^{\star}=\hat{\lambda}_{g}^{U}$, $g=1,\ldots,G$ 3. 3.

Compute $\bm{K}_{g}^{\star}$ applying the optimal truncation operator defined in Fritz et al. (2013) to $\left\{\bm{K}_{1},\ldots,\bm{K}_{G}\right\}$, under condition (33) 4. 4.

Compute $\bm{C}^{\star}=\left.\sum_{g=1}^{G}\frac{1}{\lambda_{g}^{\star}}\bm{K}^{\star}_{g}\middle/\left|\sum_{g=1}^{G}\frac{1}{\lambda_{g}^{\star}}\bm{K}^{\star}_{g}\right|^{1/p}\right.$ 5. 5.

Compute $\lambda_{g}^{\star}=\frac{1}{p}tr\left(\bm{K}^{\star}_{g}{\bm{C}^{\star}}^{-1}\right)$ 6. 6.

Set $\bm{K}_{g}=\lambda_{g}^{\star}\bm{C}^{\star}$ 7. 7.

Iterate $3-6$ until (33) is satisfied 8. 8.

Considering the spectral decomposition for $\bm{C}^{\star}=\bm{D}^{\star}\bm{A}^{\star}{\bm{D}^{\star}}^{{}^{\prime}}$, set $\hat{\lambda}_{g}^{(k+1)}=\lambda^{\star}_{g}$, $\hat{\bm{A}}_{g}^{(k+1)}=\bm{A}^{\star}$, $\hat{\bm{D}}_{g}^{(k+1)}=\bm{D}^{\star}$