Robust Improper Maximum Likelihood: Tuning, Computation, and a Comparison With Other Methods for Robust Gaussian Clustering

ABSTRACT

The two main topics of this article are the introduction of the “optimally tuned robust improper maximum likelihood estimator” (OTRIMLE) for robust clustering based on the multivariate Gaussian model for clusters, and a comprehensive simulation study comparing the OTRIMLE to maximum likelihood in Gaussian mixtures with and without noise component, mixtures of t-distributions, and the TCLUST approach for trimmed clustering. The OTRIMLE uses an improper constant density for modeling outliers and noise. This can be chosen optimally so that the nonnoise part of the data looks as close to a Gaussian mixture as possible. Some deviation from Gaussianity can be traded in for lowering the estimated noise proportion. Covariance matrix constraints and computation of the OTRIMLE are also treated. In the simulation study, all methods are confronted with setups in which their model assumptions are not exactly fulfilled, and to evaluate the experiments in a standardized way by misclassification rates, a new model-based definition of “true clusters” is introduced that deviates from the usual identification of mixture components with clusters. In the study, every method turns out to be superior for one or more setups, but the OTRIMLE achieves the most satisfactory overall performance. The methods are also applied to two real datasets, one without and one with known “true” clusters. Supplementary materials for this article are available online.

KEYWORDS:

• Cluster analysis
• EM-algorithm
• Improper density
• Mixture models
• Model-based clustering
• Robustness

1. Introduction

In this article, we introduce and investigate the “optimally tuned robust improper maximum likelihood estimator” (OTRIMLE), a method for robust clustering with clusters that can be approximated by multivariate Gaussian distributions. Its one-dimensional version was introduced in Coretto and Hennig (2010). We also present a simulation study comparing OTRIMLE and other approaches for (mostly robust) model-based clustering, which is, to our knowledge, the most comprehensive study in the field and involves a careful discussion of the issue of comparing methods based on different model assumptions.

The basic idea of OTRIMLE is to fit an improper density to the data that is made up by a Gaussian mixture density and a “pseudo mixture component” defined by a small constant density, which is meant to capture outliers in low density areas of the data. This is inspired by the addition of a uniform “noise component” to a Gaussian mixture (Banfield and Raftery 1993). Hennig (2004) showed that using an improper density improves the breakdown robustness of this approach. The OTRIMLE has been found to work well for one-dimensional data in Coretto and Hennig (2010).

As in many other statistical problems, violations of the model assumptions and particularly outliers may cause problems in cluster analysis. Our general attitude to the use of statistical models in cluster analysis is that the models should not be understood as reflecting some underlying but in practice unobservable “truth,” but rather as thought constructs implying a certain behavior of methods derived from them (e.g., maximizing the likelihood), which may or may not be appropriate in a given application (more details on the general philosophy of clustering can be found in Hennig and Liao 2013). Using a model such as a mixture of multivariate Gaussian distributions, interpreting every mixture component as a “cluster,” implies that we look for clusters that are approximately “Gaussian-shaped,” but we do not want to rely on whether the data really were generated iid by a Gaussian mixture. We are interested in the performance of such methods in situations where one may legitimately look for Gaussian-shaped clusters, even if some data points do not belong to such clusters (called “noise” in the following), and even if the clusters are not precisely Gaussian. This reflects the fact that in practice, for example, mixtures of t-distributions are used for clustering the same datasets to which Gaussian mixtures are fitted as well, interpreting the resulting clusters in the same way.

For illustration of the outlier problem in model-based clustering, we use a five-dimensional dataset in which the 170 districts of the German city of Dortmund are characterized by a number of variables, which is discussed in detail in Section 6.1. Fitting a plain Gaussian mixture with G = 4 to all five variables by R’s MCLUST package (Fraley et al. 2012), one cluster is a one-point cluster consisting only of an extreme outlier, and two further clusters fit two different varieties of moderate outliers. More than 120 districts are collected in a single cluster. The task of robust clustering is to avoid having many or even most clusters dominated by outliers, and to produce a meaningful clustering structure also among the main bulk of nonextreme observations.

A number of model-based clustering methods that can deal with outliers have been proposed in recent years. An overview of these methods is given in Section 2. The OTRIMLE is introduced and discussed in Section 3, starting from the “RIMLE,” in which the level of the improper constant density is a tuning constant. We then introduce a method for optimal tuning and discuss its computation. Section 5 presents a simulation study that uses a unified approach for defining elliptically shaped clusters with noise/outliers, presented in Section 4. The Dortmund dataset mentioned above is discussed in Section 6 along with a dataset of folk song melodies from two known regions. Some further issues including the estimation of the number of clusters are discussed in Section 7. Additional details about the example dataset (including full scatterplots), the simulation study, and computation of methods are provided in an online supplement (Coretto and Hennig in press b). Theoretical properties of the RIMLE with the tuning constant fixed are investigated in Coretto and Hennig (in press a) and cited here.

2. Methods from the Literature

In the following, assume an observed sample $\underset{\_}{x_n}=\{x_1,x_2,...,x_n\},$ where x_i is the realization of a random variable $X_i\in {\mathbb{R}}^p$ with p ⩾ 1; X₁, …, X_n iid. The goal is to cluster the sample points into G distinct groups.

2.1. Maximum Likelihood (ML) for Gaussian Mixtures (gmix)

Let φ(x; μ, Σ) be the density of a multivariate Gaussian distribution with mean vector $\mu \in {\mathbb{R}}^p$ and p × p covariance matrix Σ. Assume that the observed sample is iid drawn from the finite Gaussian mixture distribution having density (2.1) $$\begin{array}{c}\hfill m(x;\theta )=\sum _{j=1}^G{\pi }_j\varphi (x;{\mu }_j,{\Sigma }_j),\end{array}$$(2.1) where π_j ∈ [0, 1] for all j = 1, 2, …, G and ∑^G_{j = 1}π_j = 1, θ is the parameter vector containing the triplets π_j, μ_j, Σ_j for all j = 1, 2, …, G. Clustering coincides with assigning points to the mixture components based on ML parameter estimates. π_j can be interpreted as the expected proportion of points originated from the jth component. Let ${\theta }_n^{\text{ml}}$ be the ML estimator for θ, usually computed by the expectation–maximization (EM) algorithm (Dempster et al. 1977). The ML estimator under (2.1(2.1) $$\begin{array}{c}\hfill m(x;\theta )=\sum _{j=1}^G{\pi }_j\varphi (x;{\mu }_j,{\Sigma }_j),\end{array}$$(2.1) ) exists only under appropriate constraints on the covariances matrices. These constraints (which are also relevant for the methods introduced below) will be discussed in detail in Section 3. Let ${\tau }_{ij}^{\text{ml}}$ be the estimated posterior probability that the observed point x_i has been drawn from the jth mixture component, that is, (2.2) $$\begin{array}{c}\hfill {\tau }_{ij}^{\text{ml}}=\frac{{\pi }_{j,n}^{\text{ml}}\varphi (x_i;{\mu }_{j,n}^{\text{ml}},{\Sigma }_{j,n}^{\text{ml}})}{m(x_i;{\theta }_n^{\text{ml}})}\text{for}\text{all}j=1,2,...,G.\end{array}$$(2.2) The point x_i can then be assigned to the kth cluster if $k=arg{max}_{j=1,2,...,G}{\tau }_{ij}^{\text{ml}}$. This assignment method is common to all model-based clustering methods introduced here. gmix is implemented in R’s MCLUST package (Fraley et al. 2012). As illustrated in Section 1 and proven in Hennig (2004), the method can be strongly affected by outliers and deviations from the model assumptions, and we now turn to approaches that attempt to deal with this problem. For lack of space, we present these in detail in the online supplement and only give a short overview here.

2.2. ML-Type Estimator for Gaussian Mixtures with Uniform Noise (gmix.u)

Banfield and Raftery (1993) added a uniform mixture component on the smallest hyperrectangle covering the data to (2.1(2.1) $$\begin{array}{c}\hfill m(x;\theta )=\sum _{j=1}^G{\pi }_j\varphi (x;{\mu }_j,{\Sigma }_j),\end{array}$$(2.1) ), calling it “noise component” to accommodate “noise.”

2.3. ML for Mixtures of Student t-Distributions (tmix)

McLachlan and Peel (2000) replaced the Gaussian densities in (2.1(2.1) $$\begin{array}{c}\hfill m(x;\theta )=\sum _{j=1}^G{\pi }_j\varphi (x;{\mu }_j,{\Sigma }_j),\end{array}$$(2.1) ) with multivariate Student-t densities, because they have heavier tails and can therefore accommodate outliers in a better way. Observations can be declared “noise” if they lie in a low density area of the t-distribution fitting their cluster. Hennig (2004) showed that neither tmix nor gmix.u is breakdown-robust.

2.4. TCLUST

TCLUST is based on maximizing a trimmed likelihood of a “fixed partition model” with cluster weights π_j. With R = ∪^G_{j = 1}R_j, $\#\left\{R\right\}=\left[n\right(1-\alpha \left)\right]$ the number of nontrimmed points: (2.3) $$\begin{array}{c}\hfill {\theta }^{\text{tclust}}:=\underset{\theta \in \Theta ,\#\left\{R\right\}=\left[n\right(1-\alpha \left)\right]}{\arg \max }\sum _{j=1}^G\sum _{i\in R_j}\left(log{\pi }_j+log\varphi (x;{\mu }_j,{\Sigma }_j)\right).\end{array}$$(2.3) For background, see Gallegos (2002), Gallegos and Ritter (2005), and García-Escudero et al. (2008). The TCLUST methodology is implemented in R’s TCLUST package by Fritz, García-Escudero, and Mayo-Iscar (2012). Partition methods with trimming started with the trimmed k-means proposal of Cuesta-Albertos, Gordaliza, and Matrán (1997).

2.5. Further Existing Work

More approaches to robust model-based clustering can be found in the literature. Neykov et al. (2007) proposed and implemented a trimmed likelihood method. Qin and Priebe (2013) introduced an EM-algorithm adapted to maximum L_q-likelihood estimation and studied its behavior under a gross error model. References to other approaches to robust clustering are given in García-Escudero et al. (2010).

3. Optimally Tuned Robust Improper Maximum Likelihood

3.1. Robust Improper Maximum Likelihood

The robust improper maximum likelihood estimator (RIMLE) is based on the “noise component”-idea for robustification (gmix.u). The main idea is to use a pseudo-model where the noise is represented by an improper constant density over the whole Euclidean space: (3.1) $$\begin{array}{c}\hfill {\psi }_{\delta }(x,\theta )={\pi }_0\delta +\sum _{j=1}^G{\pi }_j\varphi (x;{\mu }_j,{\Sigma }_j),\end{array}$$(3.1) with π₀, π_j ∈ [0, 1] for j = 1, 2, …, G, ${\pi }_0+{\sum }_{i=1}^G{\pi }_j=1$. δ > 0 is the improper constant density (icd). The parameter vector θ contains all Gaussian parameters plus all proportion parameters including π₀, i.e., θ = (π₀, π₁, …, π_G, μ₁, …, μ_G, Σ₁, …, Σ_G). Given the sample improper pseudo-log-likelihood function (3.2) $$\begin{array}{c}\hfill l_n\left(\theta \right)=\frac{1}{n}\sum _{i=1}^{n}log{\psi }_{\delta }(x_i,\theta ),\end{array}$$(3.2) the RIMLE is defined as (3.3) $$\begin{array}{c}\hfill {\theta }_n\left(\delta \right)=\underset{\theta \in \Theta }{\arg \max }{l}_n\left(\theta \right),\end{array}$$(3.3) where Θ is an appropriate constrained parameter space discussed below. θ_n(δ) is then used to cluster points as for model-based clustering methods. Define pseudo posterior probabilities in analogy with (2.2(2.2) $$\begin{array}{c}\hfill {\tau }_{ij}^{\text{ml}}=\frac{{\pi }_{j,n}^{\text{ml}}\varphi (x_i;{\mu }_{j,n}^{\text{ml}},{\Sigma }_{j,n}^{\text{ml}})}{m(x_i;{\theta }_n^{\text{ml}})}\text{for}\text{all}j=1,2,...,G.\end{array}$$(2.2) ): $$\begin{array}{c}\hfill {\tau }_j(x_i,\theta ):=\left\{\begin{array}{cc}\hfill {\displaystyle \frac{{\pi }_0\delta }{{\psi }_{\delta }(x_i,\theta )}}& \hfill \text{if}j=0\\ \hfill {\displaystyle \frac{{\pi }_j\varphi (x_i,{\mu }_j,{\Sigma }_j)}{{\psi }_{\delta }(x_i,\theta )}}& \hfill \text{if}j=1,2,...,G;\end{array}\right.\text{for}i=1,2,...,n,\end{array}$$and assign the points based on the following rule (3.4) $$\begin{array}{c}\hfill J(x_i,\theta ):=\underset{j\in \left\{0,1,2,...,G\right\}}{\arg \max }{\tau }_j(x_i,\theta ).\end{array}$$(3.4) Fixing δ, (3.1(3.1) $$\begin{array}{c}\hfill {\psi }_{\delta }(x,\theta )={\pi }_0\delta +\sum _{j=1}^G{\pi }_j\varphi (x;{\mu }_j,{\Sigma }_j),\end{array}$$(3.1) ) does not define a proper probability model, but (3.3(3.3) $$\begin{array}{c}\hfill {\theta }_n\left(\delta \right)=\underset{\theta \in \Theta }{\arg \max }{l}_n\left(\theta \right),\end{array}$$(3.3) ) yields a useful procedure for data modeled as a proportion of (1 − π₀) of a mixture of Gaussian distributions plus a proportion of π₀ points not assigned to any meaningful cluster. Regions of high density are rather associated with clusters than with noise, so the noise regions should be those with the lowest density. This could be achieved by using the uniform density as in gmix.u, but for this the presence of gross outliers the dependence of the uniform distribution on the convex hull of the data still causes a robustness problem (Hennig 2004).

The optimization problem in (3.3(3.3) $$\begin{array}{c}\hfill {\theta }_n\left(\delta \right)=\underset{\theta \in \Theta }{\arg \max }{l}_n\left(\theta \right),\end{array}$$(3.3) ) requires that Θ is suitably defined, otherwise θ_n(δ) may not exist. As discovered by Day (1969), the Gaussian mixture likelihood can degenerate. This problem extends to (3.1(3.1) $$\begin{array}{c}\hfill {\psi }_{\delta }(x,\theta )={\pi }_0\delta +\sum _{j=1}^G{\pi }_j\varphi (x;{\mu }_j,{\Sigma }_j),\end{array}$$(3.1) ) as well. Let λ_{k, j} be an eigenvalue of Σ_j for some k = 1, 2, …, p and j = 1, 2, …, G. Take a sequence ${\left({\theta }_m\right)}_{m\in \mathbb{N}}$ such that λ_{k, j, m}→0 and μ_{j, m} = x₁, then l_n(θ_m) → +∞. There are various ways to avoid this issue. Let λ_max(θ) and λ_min(θ) be, respectively, the maximum and the minimum eigenvalues of the covariance matrices in θ, Coretto and Hennig (in press a) adopted the “eigenratio constraint” (3.5) $$\begin{array}{c}\hfill {\lambda }_{max}\left(\theta \right)/{\lambda }_{min}\left(\theta \right)\le \gamma <+\infty \end{array}$$(3.5) with fixed γ ⩾ 1. γ = 1 constrains all component covariance matrices to be spherical and equal, as in k-means clustering, while γ > 1 restricts the relative scatter discrepancy among clusters. This type of constraint has been proposed by Dennis (1981) and studied by Hathaway (1985) for one-dimensional Gaussian mixtures. EM-algorithms for computing the ML of multivariate Gaussian mixtures under (3.5(3.5) $$\begin{array}{c}\hfill {\lambda }_{max}\left(\theta \right)/{\lambda }_{min}\left(\theta \right)\le \gamma <+\infty \end{array}$$(3.5) ) have been studied by Ingrassia (2004) and Ingrassia and Rocci (2007), although asymptotic properties of the corresponding MLE have not been proved. The same constraints are used for TCLUST by García-Escudero et al. (2008). There are a number of alternative constraints, see Ingrassia and Rocci (2011) and Gallegos and Ritter (2009).

Although (3.5(3.5) $$\begin{array}{c}\hfill {\lambda }_{max}\left(\theta \right)/{\lambda }_{min}\left(\theta \right)\le \gamma <+\infty \end{array}$$(3.5) ) prevents the unboundness of the likelihood in standard mixture models and TCLUST, for RIMLE this is not enough. Points not fitted by any of the Gaussian components can still be fitted by the improper uniform component. Therefore, Coretto and Hennig (in press a) proposed an additional “noise proportion constraint,” (3.6) $$\begin{array}{c}\hfill \frac{1}{n}\sum _{i=1}^n{\tau }_0(x_i,\theta )\le {\pi }_{\text{max}},\end{array}$$(3.6) for fixed 0 < π_max < 1. The quantity n^{− 1}∑ⁿ_{i = 1}τ₀(x_i, θ) can be interpreted as the estimated proportion of noise points. Setting π_max = 0.5 just implements a familiar condition in robust statistics that at most half of the data should be classified as “outliers/noise.” The resulting restricted parameter space for RIMLE is then (3.7) $$\begin{array}{ccc}\hfill \Theta :& =& \left\{\theta :{\pi }_j\ge 0\forall j\ge 1;{\pi }_0+\sum _{j=1}^G{\pi }_j=1;\frac{1}{n}\sum _{i=1}^n{\tau }_0(x_i,\theta )\right.\hfill \\ & \le & \left.{\pi }_{\text{max}};\frac{{\lambda }_{max}\left(\theta \right)}{{\lambda }_{min}\left(\theta \right)}\le \gamma \right\}.\hfill \end{array}$$(3.7) Coretto and Hennig (in press a) showed that θ_n(δ) exists for any δ ⩾ 0 if $\#\left(\underset{\_}{x_n}\right)>G+\lceil n{\pi }_{max}\rceil$ and that θ_n(0) exists under the milder condition that $\#\left(\underset{\_}{x_n}\right)>G$. For δ = 0, the RIMLE reduces to ML for plain Gaussian mixtures. Let E_Pf(x) be the expectation of f(x) under x ∼ P. The RIMLE functional is defined as (3.8) $$\begin{array}{c}\hfill {\theta }^{★}\left(\delta \right)=\underset{\theta \in {\Theta }_G}{\arg \max }{E}_{P}log{\psi }_{\delta }(x;\theta ).\end{array}$$(3.8) Existence of (3.8(3.8) $$\begin{array}{c}\hfill {\theta }^{★}\left(\delta \right)=\underset{\theta \in {\Theta }_G}{\arg \max }{E}_{P}log{\psi }_{\delta }(x;\theta ).\end{array}$$(3.8) ), consistency of θ_n(δ) on the quotient space topology identifying all log-likelihood maxima and its breakdown point are shown in Coretto and Hennig (in press a).

3.2. Optimal Improper Density Level

Occasionally, subject matter knowledge may be available aiding the choice of δ, but such situations are rather exceptional. Here we suggest a data-dependent choice of δ. Note that δ is not treated as a model quantity to be estimated here, but rather as a tuning device to enable a good robust clustering. The aim of the RIMLE is to approximate the density of clustered regions of points when these regions look like those produced by a Gaussian distribution. We define the “optimal” δ value as minimizer of a criterion function measuring the discrepancy of the found clusters from the Gaussian prototype. Given θ_n(δ), define the clusterwise squared Mahalanobis distances to clusters’ centers as $$\begin{array}{c}\hfill d_{i,j,n}={(x_i-{\mu }_{j,n})}^{\text{'}}{\Sigma }_{j,n}^{-1}(x_i-{\mu }_{j,n}),\end{array}$$and the clusterwise weighted empirical distribution of d_{i, j, n}, (3.9) $$\begin{array}{ccc}\hfill {\mathbb{M}}_j(t;\delta )& =& \frac{1}{{\sum }_{i=1}^n{\tau }_j(x_i,{\theta }_n\left(\delta \right))}\sum _{i=1}^n{\tau }_j(x_i,{\theta }_n\left(\delta \right))\mathbf{\text{1}}\{d_{i,j,n}\le t\},\hfill \\ \hfill j& =& 1,2,...,G.\hfill \end{array}$$(3.9) In ${\mathbb{M}}_j$, the ith point’s distance is weighted according to the pseudo posterior probability that the ith observation has been generated by the jth mixture component. If the jth cluster is approximately Gaussian and μ_{j, n} and Σ_{j, n} are good approximations of its location and scatter, we expect that squared Mahalanobis distances to μ_{j, n} of the points indeed belonging to mixture component no. j (for which τ_j( · ) indicates the estimated probability) will approximate a χ²_p distribution. With χ²_p(a) being the value of the cdf of the χ²_p distribution at a, define the Kolmogorov-type distance for the jth cluster (3.10) $$\begin{array}{c}\hfill {\text{K}}_j\left(\delta \right)=\underset{i=1,...,n}{max}\left|{\mathbb{M}}_j(d_{i,j,n};\delta )-{\chi }_p^2\left(d_{i,j,n}\right)\right|.\end{array}$$(3.10) The quality of the overall Gaussian approximation is then evaluated by weighting ${\text{K}}_j(·)$ with the estimated component proportion π_{j, n}: (3.11) $$\begin{array}{c}\hfill \text{D}\left(\delta \right)=\frac{1}{{\sum }_{j=1}^G{\pi }_{j,n}}\sum _{j=1}^G{\pi }_{j,n}{\text{K}}_j\left(\delta \right).\end{array}$$(3.11) For a given constant β ⩾ 0, define the optimal icd level as (3.12) $$\begin{array}{c}\hfill {\delta }_n=\underset{\delta \in [0,{\delta }_{max}]}{\arg \min }\left[\text{D}\left(\delta \right)+\beta {\pi }_{0,n}\right].\end{array}$$(3.12) The corresponding optimally tuned RIMLE (OTRIMLE) will be denoted as θ_n(δ_n). Existence and uniqueness of δ_n are not trivial, see Section 3.3. Section 5.4 is about how D(δ) behaves as a function of δ.

The straightforward choice for β, formalizing the “Gaussian cluster” concept, is β = 0. However, often in practice it is not so important that the clusters are of as precise Gaussian shape as possible. β > 0 (but normally smaller than 1) formalizes that less Gaussian shapes of clusters are tolerated if this brings the estimated noise proportion down. As can be seen in Section 5, choosing $\beta =\frac{1}{3}$ leads to improvements if the true clusters are t-distributed. Section 5.5 gives more details on the effect of choosing β > 0.

3.3. Computation

For a fixed δ, the RIMLE can be appropriately computed using an expectation–maximization (EM) algorithm. See Coretto and Hennig (in press a) for details, and the online supplement for details and background of the following. The outcome of the EM-algorithm depends on the initialization. We used an initialization method inspired by the MCLUST software. To avoid spurious clusters, we consider as valid initial partitions only those containing at least min.pr× n observations in each cluster (min.pr=0.005, say). As a first attempt to find such a valid partition, nearest neighbor-based clutter/noise detection proposed by Byers and Raftery (1998) is applied to identify an initial noise guess. Agglomerative hierarchical clustering based on ML criteria for Gaussian mixture models proposed by Banfield and Raftery (1993) is then used for finding initial Gaussian clusters among the nonnoise. See the online supplement for the case that the found partition is not “valid.”

The OTRIMLE can be found by computing RIMLEs on a grid of δ values ranging from zero to some large enough δ_max. In practice, we solve the program (3.12(3.12) $$\begin{array}{c}\hfill {\delta }_n=\underset{\delta \in [0,{\delta }_{max}]}{\arg \min }\left[\text{D}\left(\delta \right)+\beta {\pi }_{0,n}\right].\end{array}$$(3.12) ) by the “golden section search” of Kiefer (1953) over the candidate set δ ∈ [0, δ_max]. In most numerical experiments, we found that no more than 30 RIMLE evaluations are required. δ_max can be chosen as highest density value occurring within an initialized cluster, discarding δ-values for which the RIMLE-solution ends up at the border of the parameter space.

4. Definition of “True” Clusters and Misclassification

In most simulation studies in cluster analysis, in which data are generated from mixture (or fixed partition) models, it is assumed that the “true” clusters are identified with the mixture components, and methods can then be evaluated by misclassification rates. But this can be problematic. Consider the comparison of ML estimators for Gaussian mixtures and for mixtures of t-distributions. In most applications, both approaches would be considered as potentially appropriate for doing the same thing, namely, finding clusters that are unimodal and elliptical. In applications in which the clustering is of main interest (as opposed to parameter estimation), researchers would not mind much whether the density around their cluster cores rather looks like a Gaussian or a t-distribution. But implications for which points are considered as “true outliers” versus “truly belonging to a cluster” would be different, because some points generated by a t-distribution with low degrees of freedom are indeed outlying with respect to the core of the t-distribution from which they are generated.

More generally, the identification of clusters and mixture components cannot be taken for granted. Hennig (2010) illustrated that the interpretation of Gaussian mixture components as clusters depends on whether the components are separated enough. And in robust clustering, one would often interpret a group of a few points with low density as “noise” even if they were generated by a Gaussian distribution.

We now define a “reference truth” for the mixture models that are used in our simulation study, from which misclassification rates then can be computed. For motivation consider Figure 1, which shows an artificial dataset drawn from a mixture of two Gaussians and a uniform distribution. Figure 1(a) shows the unlabeled dataset on which cluster analysis operates. Figure 1(b) shows the points labeled by the mixture components that generated them. Observe that there are three red stars (generated from the uniform noise) in the middle of the region where the two Gaussians have most of their mass. Furthermore, there are green points from the left Gaussian component with lower density that fall into the dense blue region. No method can be expected to reconstruct all the cluster memberships in such overlapping regions.

Figure 1. An artificial dataset consisting of 950 points drawn from two two-dimensional Gaussian distributions and 50 points from a uniform distribution on the square [ − 10, 10] × [ − 5, 15]. (a) unlabeled points, (b) colored according to the mixture components, (c) colors represent the two GR_α (with α = 10^{− 4}); red stars belong to NR_α.

Figure 1(c) shows what we define as the “reference truth,” defined next.

The idea is that we choose probability measures P₁, P₂, …, P_G to correspond to the G “true clusters” (implying that they “cluster,” that is, generate clearly distinguishable, although not necessarily nonoverlapping, data patterns). For each of these, we define a region of points that can be considered as nonoutliers based on a mean and covariance matrix functional, which are Fisher-consistent at the Gaussian distribution, but robust and existing for other distributions, too. This defines a region of nonoutliers of Gaussian shape. We consider all points “noise” that are outliers according to this definition for all P₁, P₂, …, P_G. Quadratic discriminant analysis assigns points to clusters that are nonoutliers with respect to more than one of the P_j. This means that points are assigned to clusters by optimal classification boundaries under the Gaussian assumption, even if the components are in fact not Gaussian. This formalizes using “Gaussian cluster prototypes” without assuming that clusters really have to be Gaussian.

To this end, let m_j and S_j be the minimum covariance determinant (MCD) center and scatter functional at P_j (Rousseeuw 1985). Cator and Lopuhaä (2012) proved existence of the MCD functional for a wide class of probability measures. The S_j can be corrected for achieving consistency at P_j equal to the Gaussian distribution (Croux and Haesbroeck 1999; Pison et al. 2002), so that when P_j is Gaussian, m_j and S_j are the corresponding mean vector and covariance matrix. Let π_j be the expected proportion of points generated from P_j. We allow ∑^G_{j = 1}π_j ⩽ 1, so that points could be generated by (noise-)distributions other than P₁, P₂, …, P_G. Define the quadratic discriminant score for assigning the point $y\in {\mathbb{R}}^p$ to the jth cluster by maximizing $$\begin{array}{ccc}\hfill \text{qs}(y;{\pi }_j,m_j,S_j):& =& log\left({\pi }_j\right)-\frac{1}{2}log(det\left(S_j\right))\hfill \\ & & -\frac{1}{2}{(y-m_j)}^{\text{'}}{S}_j^{-1}(y-m_j),\hfill \\ & & \text{for}j=1,2,...,G.\hfill \end{array}$$If clusters are indeed Gaussian, this is equivalent to (3.4(3.4) $$\begin{array}{c}\hfill J(x_i,\theta ):=\underset{j\in \left\{0,1,2,...,G\right\}}{\arg \max }{\tau }_j(x_i,\theta ).\end{array}$$(3.4) ). Consider $$\begin{array}{c}\hfill {\mathcal{E}}_{\alpha }(m_j,S_j):=\{y:{(y-m_j)}^{\text{'}}{S}_j^{-1}(y-m_j)\le {\chi }_p^2(1-\alpha )\},\end{array}$$where χ²_p(1 − α) is the 1 − α quantile of the χ²_p distribution. For a fixed α, the ellipsoid ${\mathcal{E}}_{\alpha }(m_j,S_j)$ defines the subset of ${\mathbb{R}}^p$ that hosts the jth cluster. The size of this ellipsoid is defined in terms of χ²_p(1 − α), because for Gaussian P_j, $P_j({\mathbb{R}}^p\setminus {\mathcal{E}}_{\alpha }(m_j,S_j))=\alpha$. For a fixed level α, the α-Gaussian region is defined as the union of these ellipsoids: $$\begin{array}{c}\hfill {\text{GR}}_{\alpha }:=\bigcup _{j=1}^G{\mathcal{E}}_{\alpha }(m_j,S_j),\end{array}$$and the noise region is given by ${\text{NR}}_{\alpha }:={\mathbb{R}}^p\setminus {\text{GR}}_{\alpha }$.

Definition 1

(α-Gaussian cluster memberships). Given α ∈ [0, 1), a data-generating process (DGP) with cluster parameters θ_C ≔ {(π_j, m_j, S_j), j = 1, 2, …, G}, and a dataset $\underset{\_}{x_n}:=\{x_1,x_2,...,x_n\}$; the α-Gaussian cluster memberships are given by (4.1) $$\begin{array}{ccc}\hfill {\text{AGR}}_{\alpha }(x_i;{\theta }_C):& =& \mathbf{\text{1}}\{x_i\in {\text{GR}}_{\alpha }\}\times \underset{j=1,2,...,G}{\arg \max }\hfill \\ & & j\mathbf{\text{1}}\{j={\arg \max }_{g=1,...,G}\text{qs}(y;{\pi }_g,m_g,S_g)\}.\hfill \end{array}$$(4.1)

${\text{AGR}}_{\alpha }(x_i,{\theta }_C)=0$ means that $x_i\in {\text{NR}}_{\alpha }$.

This definition is inspired by the definition of outliers with respect to a reference model as in Davies and Gather (1993) and Becker and Gather (1999). A difference is that here the parameter α does not directly control the probability of the noise region. Once α and the triples (π_j, m_j, S_j) are fixed for all j = 1, 2, …, G, the size of the noise region will depend on the degree of overlap and Gaussianity of the ellipsoids in ${\text{GR}}_{\alpha }$. α needs to be small because the idea of an outlier implies that under the Gaussian distributions outliers are very rare. We choose α = 10^{− 4}, which implies that the probability that there is at least one outlier in n = 500 iid Gaussian observations is 0.0488.

The different robust clustering methods have different implicit ways of classifying points as “noise” (noise component, trimming, outlier identification in t-distributions). To make them comparable, we use (4.1(4.1) $$\begin{array}{ccc}\hfill {\text{AGR}}_{\alpha }(x_i;{\theta }_C):& =& \mathbf{\text{1}}\{x_i\in {\text{GR}}_{\alpha }\}\times \underset{j=1,2,...,G}{\arg \max }\hfill \\ & & j\mathbf{\text{1}}\{j={\arg \max }_{g=1,...,G}\text{qs}(y;{\pi }_g,m_g,S_g)\}.\hfill \end{array}$$(4.1) ) to unify the point assignment of the methods by computing ${\text{AGR}}_{\alpha }(·)$ based on the parameters estimated by the methods, from which we assume that estimators of the triples (π_j, m_j, S_j) (cluster proportion, center, and covariance matrix) can be computed (see Section 5.2). Let the estimated cluster parameters be ${\widehat{\theta }}_{C,n}$. A misclassification rate can then be computed by applying an optimal permutation σ{ · } of cluster labels: (4.2) $$\begin{array}{ccc}\hfill \text{mcr}({\theta }_C,{\widehat{\theta }}_{C,n})& =& \underset{\sigma }{\arg \min }\frac{1}{n}\sum _{i=1}^n\mathbf{\text{1}}\{{\text{AGR}}_{\alpha }(x_i;{\theta }_C)\hfill \\ & =& \sigma \left\{{\text{AGR}}_{\alpha }(x_i;{\widehat{\theta }}_{C,n})\right\}\}.\hfill \end{array}$$(4.2) See the online supplement for computation of the MCD functional for nonnormal distributions.

5. A Comparative Simulation Study

Here, we present a comprehensive simulation study comparing the OTRIMLE with the methods introduced in Section 2.

5.1. Data-Generating Processes

The methods are compared on a total of 24 DGPs with 1000 Monte Carlo replicates each. Half of the DGPs produce two-dimensional datasets. The remaining 12 DGPs are 20-dimensional versions that are constructed adding independent 18-dimensional uncorrelated zero-means unit-variance Gaussian and/or Student-t marginals. Therefore, clusters are always only defined on the first two marginals. Note that the aim of the simulation study is not variable selection; we designed the DGPs so that clustering information is only in the first two dimensions to be able to visualize and control the clustering patterns, but we compare clustering methods that use all variables (trying variable selection methods is beyond the scope of this article). We do not think that variable selection or dimension reduction is mandatory in clustering, because the meaning of the clusters is determined by the involved variables, see the discussion rejoinder in Hennig and Liao (2013). We choose n = 1000 for two-dimensional designs and n = 2000 for the 20-dimensional versions. DGPs have been designed to test a variety of “noise patterns,” numbers of clusters G, and patterns of separation/overlap between different clusters.

We consider two main classes of DGP, namely, DGPs with a uniform noise component on the first two marginals, and DGPs that do not have a noise component. The first group includes the following setups, all of which have clusters generated from Gaussian distributions, and for p = 20 the 18 uninformative variables are Gaussian: (i) for “WideNoise” DGPs, the uniform noise component produces points that are widespread but overlap with the clustered regions entirely; (ii) for “SideNoise” DGPs the uniform noise component spreads points on a wide region that overlaps slightly with some of the clusters; (iii) in “SunSpot” DGPs there is a uniform component that produces few extremely outlying points. On the other hand, we consider DGPs that do not include a noise component (i.e., π₀ = 0). This second group can be divided into three further subgroups: (i)/(ii) in “GaussT” and “TGauss” DGPs, multivariate Student-t distributions with three degrees of freedom are used. In “GaussT” these are used as uninformative distributions for p = 20, whereas the first two clustered dimensions use Gaussians; in “TGauss” the clusters are generated by noncentral multivariate t₃-distributions and for p = 20 the 18 uninformative variables are Gaussian; (iii) in “Noiseless” DGPs, all points are drawn from Gaussian distributions.

For each of the six setups, there are variants with a lower and a higher number of clusters G, p = 2 (denoted by “l”) and p = 20 (denoted by “h”), adding up to 24 DGPs. The nomenclature used in the following puts these at the end of the setup name, that is, “TGauss.5h” refers to the “TGauss”-setup with higher G = 5 and higher P = 20. For “WideNoise,” “SideNoise,” and “GaussT,” the lower G was 2 and the higher G was 3. For “SunSpot,” “TGauss,” and “Noiseless,” the lower G was 3 and the higher G was 5. The overlap between clusters as well as the combinations of cluster shapes varied between DGPs. Full details of the definition of the DGPs are given in the online supplement together with exemplary scatterplots of the first two dimensions of a dataset from every setup.

5.2. Implementation of Methods

Table 1 summarizes settings for the compared methods. TCLUST and RIMLE/OTRIMLE are based on eigenratio constraints, but this is not the case for the MCLUST software (Fraley et al. 2012) and the available implementation of mixtures of t-distributions (McLachlan and Peel 2000). To have full comparability of the solutions, the eigenratio constraints (see Section 3.3) have been implemented by us for OTRIMLE/RIMLE, gmix, tmix, and gmix.u; the latter is computed by use of the same routine that is used for RIMLE /OTRIMLE. For TCLUST, the constraints are in the original R-package.

Table 1. Summary of the main settings for the methods under comparison.

Method	Setup
gmix.u	RIMLE with δ = 1/Vn, where Vn is volume of the smallest hyperrectangle that contains the data.
ot.rimle	OTRIMLE without penalty (β = 0).
ot.rimle.p	OTRIMLE with penalty term β = 1/3.
tclust	TCLUST with fixed trimming level set at 0.1.
ot.tclust	TCLUST with trimming level selected by the OTRIMLE criterion without penalty (β = 0).
ot.tclust.p	TCLUST with trimming level selected by the OTRIMLE criterion with penalty term β = 1/3.
tmix	ML for the Student-t mixture model (2.1(2.1) $$\begin{array}{c}\hfill m(x;\theta )=\sum _{j=1}^G{\pi }_j\varphi (x;{\mu }_j,{\Sigma }_j),\end{array}$$(2.1) ) with v = 3 for all components plus eigenratio constraints.
gmix	ML for the Gaussian mixture model plus eigenratio constraints.

For^{Figure 2}all methods, the eigenratio constraint has been set equal to 20 for each of the 24 DGPs. The latter choice is motivated by the fact that 20 is larger than the maximum true eigenratio across the designs involved in the comparison, and it is in general a value that often enables rather smooth optimization (obviously in reality the true eigenvalue ratios are not known, but for variables with comparable scales and value ranges, 20 gives the covariance matrices enough flexibility for most applications). OTRIMLE has been tested with and without the penalty term $\beta =\frac{1}{3}$ in (3.12(3.12) $$\begin{array}{c}\hfill {\delta }_n=\underset{\delta \in [0,{\delta }_{max}]}{\arg \min }\left[\text{D}\left(\delta \right)+\beta {\pi }_{0,n}\right].\end{array}$$(3.12) ), denoted by ot.rimle (β = 0) and ot.rimle.p, respectively. Other values for β, ranging from 0.1 to 0.5, have been tried, and results did not change much. A difficulty with TCLUST is that an automatic data-driven choice of the trimming level is currently not available. In tclust.f we set the trimming level to 10%. This choice is motivated by the fact the DGPs produced an average proportion of points belonging to the NR_α set in the range [0%, 23%] (see Table 1 in Coretto and Hennig in press b). Furthermore, since the trimming level plays a role similar to δ in RIMLE/OTRIMLE, there are two versions of TCLUST for which the same idea for automatic decisions the trimming level is used as proposed here for OTRIMLE, see Section 3. In ot.tclust and ot.tclust.p, the trimming level has been selected using (3.12(3.12) $$\begin{array}{c}\hfill {\delta }_n=\underset{\delta \in [0,{\delta }_{max}]}{\arg \min }\left[\text{D}\left(\delta \right)+\beta {\pi }_{0,n}\right].\end{array}$$(3.12) ) with the trimming level playing the role of δ, and the weights τ_j( · ) in ${\mathbb{M}}_j(·)$ replaced by the 0–1 crisp weights of TCLUST, again using β = 0 or $\beta =\frac{1}{3}$. For t-mixtures, in the original proposal by McLachlan and Peel (2000), the degrees of freedom are assumed to be equal across the mixture components and are estimated from the data and covariance matrix constraints are not considered. In tmix, we fix the degrees of freedom to 3 and we incorporate eigenratio constraints, as for the other methods. This is motivated as follows: (i) for some of the DGPs (particularly the SunSpot designs) constraints were needed to avoid spurious solutions; (ii) for some of the sampling designs not based on Student-t distributions, estimation of the degrees of freedom of the t-distribution produced an extremely large variability in the misclassification rates; (iii) since a number of designs are based on Student-t with 3 degrees of freedom, the decision to fix this parameter to 3 gives the t-mixture a slight advantage, which seems fair given that the majority of setups rather seem to favor noise component/trimming.

For some of the DGPs, the experience suggested that solutions may depend strongly on the initialization. To reduce the bias introduced by different initializations, all methods are initialized from the same partition, see Section 3.3 (an additional set of R functions has been provided by TCLUST’s authors to allow for this).

5.3. Results

The methods are compared using misclassification rates as defined in (4.2(4.2) $$\begin{array}{ccc}\hfill \text{mcr}({\theta }_C,{\widehat{\theta }}_{C,n})& =& \underset{\sigma }{\arg \min }\frac{1}{n}\sum _{i=1}^n\mathbf{\text{1}}\{{\text{AGR}}_{\alpha }(x_i;{\theta }_C)\hfill \\ & =& \sigma \left\{{\text{AGR}}_{\alpha }(x_i;{\widehat{\theta }}_{C,n})\right\}\}.\hfill \end{array}$$(4.2) ). These are more relevant in clustering tasks than parameter estimates. The results are graphically summarized in Figure 2, while average misclassification rates with standard errors are given in Table 2. Each square in the plot is a color-coded representation of the misclassification rate averaged over the 1000 Monte Carlo replicates for a given method-DGP pair. Further details about average misclassification rates are given in the online supplement. It also contains boxplots of the misclassification rates for all method-DGP pairs. Figure 2 shows clear evidence that using robust methods is important. The gmix method only performs well for Noiseless and some DGPs with Gaussians and t-distributions, but most other methods work well (although slightly worse at times) for these DGPs, too. tmix works well for most DGPs involving t-distributions, but for the other DGPs with noise/outliers, it is often seriously worse than gmix.u, OTRIMLE, and TCLUST.

Figure 2. Level plot representing the sample mean of the Monte Carlo distribution of misclassification rates (percentage scale) for each DGP-method pair. Each square of the plot represents the average misclassification according to the bottom gray color scale.

Table 2. Monte Carlo average misclassification rates (%) with their standard errors in brackets. Misclassification rates are computed as in (4.1(4.1) $$\begin{array}{ccc}\hfill {\text{AGR}}_{\alpha }(x_i;{\theta }_C):& =& \mathbf{\text{1}}\{x_i\in {\text{GR}}_{\alpha }\}\times \underset{j=1,2,...,G}{\arg \max }\hfill \\ & & j\mathbf{\text{1}}\{j={\arg \max }_{g=1,...,G}\text{qs}(y;{\pi }_g,m_g,S_g)\}.\hfill \end{array}$$(4.1) ). Notice that both averages (and standard errors) are reported in percentage scale.

	Method
DGP	gmix.u	ot.rimle	ot.rimle.p	tclust	ot.tclust	ot.tclust.p	tmix	gmix
WideNoise.2l	1.33(0.02)	5.00(0.29)	1.35(0.02)	12.70(0.33)	13.11(0.32)	10.20(0.32)	17.48(0.10)	18.40(0.04)
WideNoise.2h	17.36(0.04)	5.28(0.05)	5.40(0.06)	22.94(0.32)	16.20(0.33)	12.85(0.25)	18.90(0.03)	17.43(0.04)
WideNoise.3l	0.34(0.01)	0.42(0.01)	0.40(0.01)	0.44(0.01)	0.38(0.01)	0.39(0.01)	3.10(0.02)	7.79(0.13)
WideNoise.3h	7.59(0.18)	3.18(0.09)	4.27(0.09)	0.72(0.01)	1.31(0.02)	3.04(0.06)	5.91(0.04)	8.56(0.15)
SideNoise.2l	0.01(0.00)	0.03(0.01)	0.04(0.02)	0.18(0.01)	0.02(0.00)	0.01(0.00)	10.76(0.11)	16.68(0.16)
SideNoise.2h	0.03(0.01)	0.09(0.02)	0.15(0.03)	0.20(0.01)	0.20(0.01)	0.66(0.02)	20.03(0.03)	12.61(0.15)
SideNoise.3l	0.06(0.00)	0.08(0.01)	0.11(0.03)	0.19(0.01)	0.07(0.00)	0.06(0.00)	18.89(0.47)	36.60(0.40)
SideNoise.3h	0.13(0.00)	0.19(0.03)	0.24(0.04)	0.25(0.01)	0.21(0.01)	0.53(0.02)	34.23(0.16)	31.59(0.31)
Sunspot.3l	0.11(0.00)	0.12(0.00)	0.11(0.00)	0.21(0.00)	0.13(0.01)	2.17(0.36)	34.65(0.11)	35.38(0.19)
Sunspot.3h	0.24(0.00)	0.24(0.00)	0.24(0.00)	21.38(0.95)	27.19(1.00)	36.83(1.01)	34.75(0.06)	40.00(0.37)
Sunspot.5l	3.32(0.12)	3.39(0.12)	3.37(0.12)	11.25(0.13)	10.01(0.14)	9.53(0.14)	10.01(0.10)	7.92(0.16)
Sunspot.5h	11.72(0.10)	11.63(0.10)	11.65(0.10)	11.88(0.10)	11.63(0.10)	11.61(0.10)	13.84(0.12)	13.30(0.13)
TGauss.3l	1.84(0.02)	4.56(0.09)	0.90(0.01)	0.81(0.01)	1.63(0.03)	0.95(0.01)	3.28(0.02)	3.33(0.06)
TGauss.3h	0.49(0.01)	0.53(0.01)	0.51(0.01)	0.77(0.01)	0.51(0.01)	0.49(0.01)	0.83(0.01)	0.85(0.02)
TGauss.5l	1.77(0.02)	4.65(0.07)	1.20(0.02)	1.21(0.01)	2.55(0.04)	1.32(0.02)	3.37(0.02)	3.59(0.05)
TGauss.5h	0.92(0.02)	0.97(0.01)	0.94(0.01)	1.39(0.01)	0.96(0.01)	0.89(0.01)	1.43(0.02)	1.65(0.04)
GaussT.2l	0.57(0.01)	0.57(0.01)	0.56(0.01)	28.77(0.33)	21.50(0.48)	16.52(0.48)	0.83(0.02)	0.57(0.01)
GaussT.2h	32.45(0.05)	33.60(0.14)	44.91(0.10)	42.42(0.13)	41.72(0.15)	43.13(0.13)	31.32(0.16)	33.17(0.03)
GaussT.3l	0.11(0.00)	0.11(0.00)	0.11(0.00)	0.27(0.01)	0.13(0.00)	0.12(0.00)	0.17(0.01)	0.11(0.00)
GaussT.3h	12.32(0.12)	3.69(0.03)	3.44(0.03)	9.87(0.03)	12.22(0.04)	6.85(0.08)	18.34(0.03)	36.23(0.49)
Noiseless.3l	0.12(0.00)	0.12(0.00)	0.12(0.00)	0.28(0.01)	0.15(0.01)	0.13(0.00)	0.18(0.01)	0.12(0.00)
Noiseless.3h	0.25(0.00)	0.25(0.00)	0.25(0.00)	0.38(0.00)	0.28(0.00)	0.28(0.00)	0.29(0.00)	0.25(0.00)
Noiseless.5l	4.53(0.08)	4.47(0.08)	4.43(0.08)	14.43(0.13)	13.33(0.13)	13.14(0.14)	8.23(0.14)	4.58(0.09)
Noiseless.5h	12.39(0.10)	12.63(0.10)	12.62(0.10)	16.69(0.12)	16.13(0.12)	16.10(0.12)	11.92(0.08)	12.37(0.10)

gmix.u performs relatively well, although for a number of DGPs it suffers strongly from high dimensionality. For WideNose in two dimensions, gmix.u will equal in many cases a proper ML estimator, so the method should be advantageous here, and gmix.u is indeed best for these DGPs. However, for the 20-dimensional WideNoise, taking a uniform distribution over the smallest hyperrectangle containing the data can no longer be the ML estimate for the noise-generating mixture component, and its performance deteriorates strongly.

Regarding the TCLUST methods, even though the automatic trimming level of ot.tclust and ot.tclust.p did not always improve the results, it demonstrated to provide a reasonable choice for the trimming level. In fact, for all situations where the true average noise proportion is about equal to the trimming level of tclust, performances of tclust, ot.tclust, and ot.tclust.p are very similar, meaning that the OTRIMLE criterion is a good starting point for fixing the trimming rate. Compared to the RIMLE-type methods, the performance of TCLUST suffers in situations where there is a considerable degree of overlap between clusters. For DGPs with overlap (such as WideNoise.2, SunSpot.5, GaussT.2, and Noiseless.5), the misclassification rate of TCLUST is completely dominated by misclassifications between clusters (to see this consider Tables 2 and 4 in the online supplement). The reason is that the TCLUST parameters are based on a classification-type likelihood, which relies on the separation between clusters. The TCLUST also seems to not tolerate the large number of Student-t marginals of GaussT.2h and GaussT.3h. TCLUST performs well for a number of DGPs and is clearly best in WideNoise.2h.

The OTRIMLE methods show a very good overall performance. They produce high misclassification rates only for some 20-dimensional DGPs for which all methods are in trouble (performances for GaussT.2h are generally bad with even tmix, the best method there, producing an average misclassification rate of more than 30%), and they are best for a number of DGPs, particularly 20-dimensional WideNoise and some TGauss-DGPs. The comparison between ot.rimle and ot.rimle.p is mixed (as between ot.tclust and ot.tclust.p), with $\beta =\frac{1}{3}$ improving matters clearly for TGauss.2l and TGauss.3l (the shape of the t-distribution encourages ot.rimle to assign too many points to the noise; see Tables 2 and 3 in the online supplement, but being significantly worse for WideNoise.3h.

Comparing OTRIMLE with gmix.u, there are a number of DGPs for which gmix.u has a slightly lower misclassification rate than one or both of ot.rimle and ot.rimle.p. In all of these DGPs, all of gmix.u, ot.rimle, and ot.rimle.p basically produce the same clustering structure, with disagreement only about the classification of some borderline points. Differences are more substantial in the setups in which gmix.u is worse. For WideNoise.2h, gmix.u in most cases does not detect any noise, so that one of the clusters consists mainly of noise. For WideNoise.3h, sometimes all or much noise is merged into a cluster, with some impact on the clustering structure. For GaussT.3h, substantial amounts of outliers are integrated into the clusters.

5.4. Behavior of D(δ)

Here, we investigate the behavior of $\text{D}\left(\delta \right)$ as a function of δ via Monte Carlo experiments under various DGPs from Section 5. For each of these DGPs, we produced 100 independent samples. We computed $\text{D}\left(\delta \right)$ for a grid of δ values taken from an interval [2.22 × 10^{− 308}, 1], adding δ = 0. In Figure 3^{Figure 4}, we report Monte Carlo averages±standard errors for $\text{D}\left(\delta \right)$ for two selected DGPs: WideNoise.3h and GaussT.3h. These are defined in Section 5.1, both with p = 20. The main difference is that noise is produced by a two-dimensional uniform distribution in WideNoise.3h, whereas in GaussT.3h dimensions 3–20 are from centered t₃ distributions, generating some outliers. Figure 3 reports the behavior of $\text{D}\left(\delta \right)$ for log (δ) > −200. For smaller values of δ (including δ = 0) the behavior of the curves was basically constant. In both graphs there is a clear minimum, although for GaussT.3h this minimum lies on the border. For WideNoise.3h the OTRIMLE criterion has a nice convex behavior around its minimum. In GaussT.3h, in dimensions 3–20 the distributional shape of the clusters deviated from Gaussianity by heavier tails, although the core of the distribution looks similar to a Gaussian. D( · ) then enforces a Gaussian shape by assigning many points to the noise component. The result is that π_{0, n} becomes large, and the optimal δ happens at point where larger values of δ do not produce parameter estimates within the constrained set anymore (unless the constraint is enforced). The latter is a reason for the use of the penalty term in (3.12(3.12) $$\begin{array}{c}\hfill {\delta }_n=\underset{\delta \in [0,{\delta }_{max}]}{\arg \min }\left[\text{D}\left(\delta \right)+\beta {\pi }_{0,n}\right].\end{array}$$(3.12) ) (see also Section 5.5). For the remaining 22 DGPs of Section 5.1, we found similar patterns (mostly similar to Widenoise.3h here). Another observation in Figure 3 is that around the minimum D(δ) seems to be quite stable for different datasets from the same design.

Figure 3. Monte Carlo average for OTRIMLE criterion D(δ) (blue solid line) ± standard errors (dotted lines) computed over a grid of values for δ ∈ {0}∪[2.22 × 10^{− 308}, 1]. The two plots refer to DGPs “WideNoise.3h” and “GaussT.3h” in Section 5.

5.5. Effect of β

For all 24 DGPs considered in Section 5.1, we also investigated the behavior of the noise proportion π_{0, n} as a function of β, see (3.12(3.12) $$\begin{array}{c}\hfill {\delta }_n=\underset{\delta \in [0,{\delta }_{max}]}{\arg \min }\left[\text{D}\left(\delta \right)+\beta {\pi }_{0,n}\right].\end{array}$$(3.12) ). For each design, we produced 100 independent samples, and for each of these we computed the OTRIMLE solution θ_n(δ_n) for various values of β ∈ [0, 1]. Figure 4 reports the Monte Carlo average of the estimated noise proportion π_{0, n} ± standard error. For both WideNoise.3h and GaussT.3h, an increase in β reduces smoothly the estimated noise proportion. There is some difference, however, in scales. The impact of β is much stronger for GaussT.3h, and the same happens for all those sampling designs in which within-cluster distributions deviate from Gaussianity. Results for other DGPs are quite similar.

Figure 4. Monte Carlo average for the estimated noise proportion π_{0, n} (blue solid line) ± standard errors (dotted lines) computed over a grid of values for β ∈ [0, 1]. The two plots refer to DGPs called “WideNoise.3h” and “GaussT.3h” in Section 5.

The discovery of the overall clustering structure was never affected by changing β between 0 and 0.5 for the DGPs from Section 5. In the majority of cases, there was no change of the clustering at all. The only difference was that larger β sometimes produced a lower percentage of data classified as noise. In the case that t₃-distributions were involved, this was good, because with β = 0 only fairly small central cores (60% –70%) of the points generated by t-distributions were assigned to the clusters, whereas for larger β only points were classified as noise that were really quite “outlying.” On the other hand, with a larger β, in data from DGPs with Gaussian clusters and some noise that was not clearly separated from the clusters, more “true” noise points were assigned to the clusters. There is no objective rule for what percentage of points from a t-distribution should be considered as “truly” outlying, and therefore it needs to be decided by the user how “tolerant” OTRIMLE is desired to be toward heavier distributional tails than the Gaussian ones.

Figure 5^{Figure 6}shows a situation in which the choice of β affects the clustering a lot. The dataset consists of 100 observations each of $\mathcal{N}(0,1)$ and $\mathcal{N}(3,1)$ along the x-axis and 12 observations from $\mathcal{N}(12,25)$. OTRIMLE was fitted with G = 2. The first two mixture components are not very well separated. β = 0 does not penalize noise, and therefore the observations from the third component are declared “noise” and the first two components are separated. However, $\beta =\frac{1}{3}$ merges the first two components and declares the third one the second cluster. The “switching point” between these two ways of “interpreting” the clustering structure is at about β = 0.3; larger values of β do not change the clustering anymore. Despite the “true” G being 3 here, regarding interpretation, depending on the application it may well make sense to either treat the smallest mixture component as noise/outliers, or to merge the first two components to a single cluster. β tunes the method to rather produce noise, or to rather tolerate not-so-Gaussian components in cases like this.

Figure 5. Clusterings with β = 0 and $\beta =\frac{1}{3}$ on artificial data example.

6. Applications

In this section, we apply OTRIMLE and the alternative methods mentioned before to two real datasets. The first example does not come with any “ground truth,” whereas the second one has “true” classes.

6.1. Dortmund Data

We here analyze a dataset giving information about 170 districts of the German city of Dortmund, which is described in Sommerer and Weihs (2005). We used five sociological key variables and transformed them in such a way that fitting Gaussian distributions within clusters makes sense. The resulting variables are the logarithm of the unemployment rate (“unemployment”), the birth/death balance divided by number of inhabitants (“birth.death”), the migration balance divided by number of inhabitants (“moves.in.out”), the logarithm of the rate of employees paying social insurance (“soc.ins.emp”), and the square root of the number of inhabitants (“inhabitants”).

All variables were centered and standardized by the median absolute deviation. Figure 6 shows a scatterplot of birth.death and moves.in.out. To deal with overplotting, an existing extreme outlier with values ≈ ( − 200, 50) is not shown. Figure 7 shows a scatterplot of unemployment and soc.ins.emp. Figure 7 shows some more moderate outliers. The left side of Figure 6 shows a clustering from fitting a plain Gaussian mixture with G = 4 to all five variables by R’s MCLUST package. Cluster 4 is a one-point cluster made of the extreme outlier. Clusters 1 and 3 basically fit two different varieties of moderate outliers, whereas all the more than 120 districts that are not extreme regarding these two variables are put together in a single cluster. Clearly, it would be more desirable to have a clustering that is not dominated so much by a few odd districts, given that there is some meaningful structure among the other districts. Such a clustering is produced by the OTRIMLE method, shown on the right side of Figure 6 and on the left side of Figure 7. The clustering is nicely interpretable with cluster no. 3 collecting a group of districts with higher migration balance and very scattered birth/death rate, cluster no. 1 being a high variation cluster characterized by high unemployment or high number of employees paying social insurance, cluster no. 2 being a homogenous group with medium number of employees paying social insurance and rather high but not very high unemployment, and cluster no. 4 collecting most districts with low values on both of these variables. For this dataset, values between β = 0 and β = 0.5 yield the same clustering (with eigenratio constraint γ = 20), despite the fact that β = 0 leads to log (δ) = −11.5 and π_{0, n} = 0.055, whereas for β = 0.5, (3.12(3.12) $$\begin{array}{c}\hfill {\delta }_n=\underset{\delta \in [0,{\delta }_{max}]}{\arg \min }\left[\text{D}\left(\delta \right)+\beta {\pi }_{0,n}\right].\end{array}$$(3.12) ) produces log (δ) = −11.9 and π_{0, n} = 0.054.

Figure 6. Scatterplot of birth.death and moves.in.out from Dortmund dataset with MCLUST clustering (left) and OTRIMLE clustering (right) with G = 4.

Figure 7. Scatterplot of soc.ins.emp and unemployment from Dortmund dataset with OTRIMLE clustering (left) and TCLUST clustering with trimming rate 7% (right) with G = 4.

Looking at the methods introduced in Section 2, it turns out that substantial disagreement may exist between different robust clustering methods. Applying the methods implemented in tclust and discussed in García-Escudero et al. (2011), α = 0.07 was found to be a good trimming rate for TCLUST with G = 4 here. The resulting clustering is compared with OTRIMLE’s in Figure 7. Although what was trimmed is almost identical to OTRIMLE’s “noise,” the clustering is somewhat different, with TCLUST’s clusters no. 1 and 4 ranging into the area of “high variation characterized by high unemployment or high number of employees paying social insurance” as mainly represented by cluster no. 3 of TCLUST and cluster no. 1 of OTRIMLE. It is hard to interpret OTRIMLE’s cluster no. 4 using any pair of variables. This can be seen in the online supplement (Coretto and Hennig in press b), as well as solutions of the other methods.

6.2. Folk Song Data

The second dataset was provided by Daniel Müllensiefen. The observations are 776 folk song melodies, 586 of which are from Luxembourg and the remaining 190 are from Warmia in Poland. These are the “true” classes. The melodies are originally from the ESAC melody database (Schaffrath 1992). The 18 features (see the online supplement (Coretto and Hennig in press b) for a list) were computed by the software “FANTASTIC” (Müllensiefen 2009).

Visual inspection reveals that there are many unusual melodies, that is, outliers in the dataset. The main bulks of melodies from Luxembourg and Warmia differ systematically from each other, although there is much overlap and no strong separation. For measuring to what extent clusterings computed with G = 2 coincided with the two regions, we used the adjusted Rand index (ARI; Hubert and Arabie 1985) with an expected value of 0 for two random clusterings and a maximum of 1 for perfect agreement. OTRIMLE with settings as above (β = 0, γ = 20) classified 36.9% of the observations as “noise.” The ARI between the OTRIMLE solution and the original regions is 0.155. For this (as for the other clustering methods), the OTRIMLE solution was interpreted as a three-cluster solution with “noise” as third cluster. Default MCLUST yields an ARI = −0.045, MCLUST with noise yields ARI = −0.017, ot.tclust yields ARI = 0.016 (the original TCLUST function with trimming rate 0.369 as suggested by OTRIMLE above achieves ARI = 0.089), and tmix yields ARI = 0.083. OTRIMLE’s ARI-value, though clearly better than that of the other methods, is not particularly high, but computing the ARI only on the observations that were not classified as noise achieves ARI = 0.392 (the solution with $\beta =\frac{1}{3}$ is slightly worse here but slightly better above regarding the ARI including the noise points), which suggests that there is a clear correspondence between OTRIMLE’s clustering and melodies that are typical for the regions.

7. Concluding Remarks

Despite our effort to make the simulation study fair, ultimately it would be good to have comparisons of methods run by researchers who did not have their hand in the design of any of the methods. Every method was best for certain DGPs in the simulation study, and simulation studies could be designed that make any method “win.” Readers need to make up their own mind about to what extent our study covered situations that are important to them. One of our major aims was to confront all methods with DGPs that do not exactly match their model assumptions, but for which the methods nevertheless could be legitimately used. In fact, we incorporated crucial ideas from both MCLUST (initialization) and TCLUST (eigenvalue ratio constraints), and the combination of these ideas used here could actually be beneficial for all methods. The methods presented in this article will soon be available in the new R-package OTRIMLE.

The problem of determining the number of clusters G is very important in practice. Here are two possible approaches for RIMLE. First, the most popular approach for fitting plain mixture models, namely, the Bayesian Information Criterion, can be used (treating RIMLE/OTRIMLE improper constant density as a proper one), Fraley and Raftery (1998). Second, G could be decided in an exploratory way by monitoring the changes of the pseudo-likelihood over different values of G in a similar way to what is done for TCLUST in García-Escudero et al. (2011). Investigating these approaches in depth is beyond the scope of this article.

Supplementary Materials

Supplementary materials contain: (i) additional details about methods and algorithms; (ii) detailed definitions of the sampling designs for the simulation study along with boxplots and summary tables for misclassification rates (iii) additional plots for real data applications; (iv) R software with instructions to reproduce the comparative simulation study and applications to real data.

Funding

The authors gratefully acknowledge financial support from MIUR research grant PRIN 2010J3LZEN, and the use of the high-performance computing infrastructure funded by University of Salerno (research program ASSA098434). The authors gratefully acknowledge the support from the EPSRC grant EP/K033972/1.