The Cellwise Minimum Covariance Determinant EstimatorThanks: To appear, Journal of the American Statistical Association.Thanks: Corresponding author, peter@rousseeuw.net .

Any practicing statistician or data scientist knows that real data sets often contain outliers. One definition of outliers says that they are cases that do not obey the fit suggested by the majority of the data, which raises suspicion that they may have been generated by a different mechanism. Since cases typically correspond to rows of the data matrix, they are sometimes called rowwise outliers. They may be the result of gross errors, but they can also be nuggets of valuable information. In either case, it is important to find them. In computer science this is called anomaly detection, and in some areas it is known as exception mining. In statistics several approaches were tried, such as testing for outliers and the computation of outlier diagnostics. In our experience the approach working best is that of robust statistics, which aims to fit the majority of the data first, and then flags outliers by their large deviation from that fit.

In this paper we focus on single-class multivariate numerical data without a response variable (although the results are relevant for classification and regression too). The goal is to robustly estimate the central location of the point cloud as well as its covariance matrix, and at the same time flag the outliers that may be present. The underlying model is that the data come from a multivariate Gaussian distribution, after which some data has been replaced by outliers that can be anywhere.

The Minimum Covariance Determinant (MCD) estimator introduced by Rousseeuw (1984); Rousseeuw (1985) is highly robust to casewise outliers. Its definition is quite intuitive. Take an integer $h$ that is at least half the sample size $n$. We then look for the subset containing $h$ cases such that the determinant of its usual covariance matrix is as small as possible. The resulting robust location estimate is then the mean of that subset, and the robust covariance matrix is its covariance matrix multiplied by a consistency factor. One can show that the estimates are not overly affected when there are fewer than $n-h$ outlying cases. The MCD became computationally feasible with the algorithm of Rousseeuw and Van Driessen (1999), followed by even faster algorithms by Hubert et al. (2012) and De Ketelaere et al. (2020). Copt and Victoria-Feser (2004) computed the MCD for incomplete data. The MCD has also been generalized to high dimensions (Boudt et al., 2020), and to non-elliptical distributions using kernels (Schreurs et al., 2021). For a survey on the MCD and its applications see Hubert et al. (2018). The MCD is available in the procedure ROBUSTREG in SAS, in SAS/IML, in Matlab’s PLS Toolbox, and in the R packages robustbase (Maechler et al., 2022) and rrcov (Todorov, 2012) on CRAN. In Python one can use MinCovDet in scikit-learn (Pedregosa et al., 2011).

In recent times a different outlier paradigm has gained prominence, that of cellwise outliers, first published by Alqallaf et al. (2009). It assumes that the data were generated from a certain distributional model, after which some individual cells (entries) were replaced by other values. The difference between the casewise and the cellwise paradigm is illustrated in Figure 1. In the left panel the outlying cases are shown as black rows. In the panel on the right the cellwise outliers correspond to fewer black squares in total, but together they contaminate over half of the cases, so the existing methods for casewise outliers may fail.

In reality we do not know in advance which cells in the right panel of Figure 1 are outlying (black), unlike the simpler problem of incomplete data where we do know which cells are missing. When the variables have substantial correlations, the cellwise outliers need not be marginally outlying, and then it can be quite hard to detect them. Van Aelst et al. (2011) proposed one of the first detection methods. Rousseeuw and Van den Bossche (2018) predict the values of all cells and flag the cells that differ much from their prediction.

There has been some work on estimating the underlying covariance matrix in the presence of cellwise outliers. One approach is to compute robust covariances between each pair of variables, and to assemble them in a matrix. To estimate these pairwise covariances, Öllerer and Croux (2015) and Croux and Öllerer (2016) use rank correlations. Tarr et al. (2016) instead use the robust pairwise correlation estimator of Gnanadesikan and Kettenring (1972) in combination with the robust scale estimator $Q_{n}$ of Rousseeuw and Croux (1993). As the resulting matrix is not necessarily positive semidefinite (PSD), they then compute the nearest PSD matrix by the method of Higham (2002). Raymaekers and Rousseeuw (2021a) obtain a PSD covariance matrix by transforming (‘wrapping’) the original data variables.

Many cellwise robust methods were developed for settings such as principal components (Hubert et al., 2019), discriminant analysis (Aerts and Wilms, 2017), clustering (García-Escudero et al., 2021), graphical models (Katayama et al., 2018), low-rank approximation (Maronna and Yohai, 2008), regression (Öllerer et al., 2016; Filzmoser et al., 2020), and variable selection (Su et al., 2021). Also, isolated outliers in functional data (Hubert et al., 2015) can be seen as cellwise outliers.

In the next section we introduce the cellwise MCD estimator. It is the first method with a single objective that combines detection and estimation, unlike some existing methods which do detection and estimation separately. Because of this cellMCD has provable cellwise breakdown properties, see section 3. There we also derive its consistency. Section 4 describes its algorithm, and proves that it converges. It is faster than the earlier methods. Some illustrations on real data are shown in section 5. The performance of the method is studied by simulation in section 6, indicating that it is very robust against adversarial contamination. Section 7 concludes with a discussion.

We first note that the casewise MCD can be reformulated in terms of likelihood. The likelihood of a $d$-variate Gaussian distribution is

where $\bm{\mu}$ is a column vector, $\bm{\Sigma}$ is a positive definite matrix, and the Mahalanobis distance is $\MD(\bm{x},\bm{\mu},\bm{\Sigma})=\sqrt{(\bm{x}-\bm{\mu})^{\top}\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu})}$. For a sample $\bm{x}_{1},\ldots,\bm{x}_{n}$ we put $L(\bm{x}_{i},\bm{\mu},\bm{\Sigma}):=-2\ln(f(\bm{x}_{i},\bm{\mu},\bm{\Sigma}))$ so the maximum likelihood estimator (MLE) of $(\bm{\mu},\bm{\Sigma})$ minimizes

Let us now look for a subset $H\subset\{1,...,n\}$ with $h$ elements which minimizes (2) where the sum is only over $i$ in $H$. We can also write this with weights $w_{i}$ that are 0 or 1 in the objective $\sum_{i=1}^{n}w_{i}L(\bm{x}_{i},\bm{\mu},\bm{\Sigma})$, so we minimize

For the minimizing set of weights $w_{i}$ we know from maximum likelihood that $\bm{\widehat{\mu}}$ is the mean of the $\bm{x}_{i}$ in $H$, so it is the weighted mean of all $\bm{x}_{i}$ , and similarly

But then the third term of (3) becomes

which is constant, and so is the second term. Therefore minimizing (3) is equivalent to minimizing the determinant of (4), which is the definition of the casewise MCD.

In the context of incomplete data, Dempster et al. (1977) and others defined the observed likelihood. Let us denote the missingness pattern of the $n\times d$ data matrix $\bm{X}$ by the $n\times d$ matrix $\bm{W}$ with entries $w_{ij}$ that are 0 for missing $x_{ij}$ and 1 otherwise. Its rows $\bm{w}_{i}$ take the place of the scalar weights $w_{i}$ in (3). For the Gaussian model the observed likelihood of the $i$th observation (Little and Rubin, 2020) is given by:

in which

is called the partial Mahalanobis distance by Danilov et al. (2012). Here $\bm{x}_{i}^{(\bm{w}_{i})}$ is the vector with only the entries for which $w_{ij}=1$, and similarly for $\bm{\mu}^{(\bm{w}_{i})}$. The matrix $\bm{\Sigma}^{(\bm{w}_{i})}$ is the submatrix of $\bm{\Sigma}$ containing only the rows and columns of the variables $j$ with $w_{ij}=1$. Finally, $d^{(\bm{w}_{i})}$ is the dimension of $\bm{x}_{i}^{(\bm{w}_{i})}$, i.e. the number of non-missing entries of $\bm{x}_{i}$ . By convention, a case $\bm{x}_{i}$ consisting exclusively of NA’s has $d^{(\bm{w}_{i})}=0$, $\MD(\bm{x}_{i},\bm{w}_{i},\bm{\mu},\bm{\Sigma})=0$ and $|\bm{\Sigma}^{(\bm{w}_{i})}|=1$. Putting $L(\bm{x}_{i},\bm{w}_{i},\bm{\mu},\bm{\Sigma}):=-2\ln(f(\bm{x}_{i},\bm{w}_{i},\bm{\mu},\bm{\Sigma}))$ we see that maximizing the observed likelihood of the entire data set comes down to minimizing

This maximum likelihood estimate of $(\bm{\mu},\bm{\Sigma})$ is typically computed by the EM algorithm (Dempster et al., 1977).

When constructing a cellwise MCD, the matrix $\bm{W}$ now describes which cells are flagged: a flagged cell $x_{ij}$ gets $w_{ij}=0$. The notations $\bm{x}_{i}^{(\bm{w}_{i})}$, $d^{(\bm{w}_{i})}$, $\bm{\mu}^{(\bm{w}_{i})}$, and $\bm{\Sigma}^{(\bm{w}_{i})}$ are interpreted analogously. The matrix $\bm{W}$ is not given in advance, but will be obtained through the estimation procedure. Now $h$ can no longer apply to the number of unflagged cases. Instead, we apply it to the number of unflagged cells per column. We could minimize

over $(\bm{\mu},\bm{\Sigma},\bm{W})$. The first constraint says that the smallest eigenvalue of $\bm{\Sigma}$ is at least as large as a number $a>0$, where the eigenvalues of $\bm{\Sigma}$ are denoted as $\lambda_{1}(\bm{\Sigma})\geqslant\ldots\geqslant\lambda_{d}(\bm{\Sigma})$. This ensures that $\bm{\Sigma}$ is nonsingular, which is required to compute Mahalanobis distances. In the second constraint, $||\bm{W}_{.j}||_{0}$ is the number of nonzero entries in the $j$-th column of $\bm{W}$. Note that we should not choose $h$ too low. Whereas for the casewise MCD we can take $h$ as low as $0.5n$, that would be ill-advised here because it could happen that two variables $j$ and $k$ do not overlap in the sense that $w_{ij}w_{ik}=0$ for all $i$, making it impossible to estimate their covariance. We will impose that $h\geqslant 0.75n$ throughout.

However, minimizing (8) typically treats too many cells as outlying. This is because a value of $h$ that is suitable for one variable may be too low for another, and we do not know ahead of time which variables have many outlying cells and which have few or none. To avoid flagging too many cells, we add a penalty counting the number of flagged cells in each column. The objective function of the cellwise MCD (cellMCD) then becomes

The notation $||\bm{1}_{d}-\bm{W}_{.j}||_{0}$ stands for the number of nonzero elements in this vector, so the number of zero weights in column $j$ of $\bm{W}$, i.e. the number of flagged cells in column $j$ of $\bm{X}$. The constants $q_{j}$ for $j=1,\ldots,d$ are computed (in Section 4) from the desired percentage of flagged cells in the absence of contamination. At the same time we keep the robustness constraint that $||\bm{W}_{.j}||_{0}\geqslant h$. Combining a penalty term with a $||.||_{0}$ constraint is not new, see the work of She et al. (2022) on casewise robust regression. In our context, the constraint $||\bm{W}_{.j}||_{0}\geqslant h$ will ensure the robustness of the estimator (through Proposition 2 below), whereas the penalty term $\sum_{j}q_{j}||\bm{1}_{d}-\bm{W}_{.j}||_{0}$ discourages flagging too many cells, which improves the estimation accuracy at clean data as seen in simulations.

The cellMCD method is the first cellwise robust technique that combines the fitting of the parameters and the flagging of outlying cells $(\bm{W})$ in one objective function. The constraint $||\bm{W}_{.j}||_{0}\geqslant h$ for $j=1,\ldots,d$ says that we require at least $h$ unflagged cells in each column. In order to avoid a singular covariance matrix, we obviously need $h>d$. Combining these inequalities we obtain $n>4d/3$ . But the curse of dimensionality implies that many spurious structures can be found in increasing dimensions, so we want a more comfortable ratio of cases per dimension. For the casewise MCD the rule of thumb is $n/d\geqslant 5$ (Rousseeuw and van Zomeren, 1990), and we will require that here too.

The cellMCD method defined by (9) is equivariant for permuting the cases, for shifting the data, and for multiplying the variables by nonzero constants. But unlike the casewise MCD it is not equivariant under general nonsingular linear transformations, or even orthogonal transformations. This is because cells are intimately tied to the coordinate system, and an orthogonal transformation changes the cells. This is an important difference between the casewise and cellwise approaches. For instance, consider the standard multivariate Gaussian model in dimension $d=4$ with the suspicious point $(10,0,0,0)$. By an orthogonal transformation of the data, this point can be moved to $(\sqrt{50},\sqrt{50},0,0)$ or to $(5,5,5,5)$. The casewise MCD is equivariant to such transformations and will still flag the same case. But in the cellwise paradigm $(10,0,0,0)$ has one outlying cell, $(\sqrt{50},\sqrt{50},0,0)$ has two, and $(5,5,5,5)$ has four, so cellMCD will react differently, as it should.

Alqallaf et al. (2009) define the cellwise breakdown value of a location estimator. Here we will focus on finite-sample breakdown values in the sense of Donoho and Huber (1983) and Lopuhaä and Rousseeuw (1991). The finite-sample cellwise breakdown value of an estimator $\bm{\widehat{\mu}}$ at a dataset $\bm{X}$ is given by the smallest fraction of cells per column that need to be replaced to carry the estimate outside all bounds. Formally, let $\bm{X}$ be a dataset of size $n$, and denote by $\bm{X}^{m}$ any corrupted sample obtained by replacing at most $m$ cells in each column of $\bm{X}$ by arbitrary values. Then the finite-sample cellwise breakdown value of a location estimator $\bm{\widehat{\mu}}$ at $\bm{X}$ is given by

Analogously to the casewise setting, we can also define the cellwise explosion breakdown value of a covariance estimator $\bm{\widehat{\Sigma}}$ as

Moreover, we define the cellwise implosion breakdown value of $\bm{\widehat{\Sigma}}$ as

The definitions of the corresponding casewise breakdown values are very similar, the only difference being that the corrupted samples, let us call them $\bm{\widetilde{X}}^{m}$, are obtained by replacing at most $m$ rows of $\bm{X}$ by arbitrary rows. If we denote the casewise breakdown values by $\delta^{*}_{n}$, $\delta^{+}_{n}$ and $\delta^{-}_{n}$ we can formulate the following simple but useful result:

The proof consists of realizing that the casewise contaminated samples $\bm{\widetilde{X}}^{m}$ can be seen as cellwise contaminated samples $\bm{X}^{m}$. It is thus generally true that the cellwise breakdown value is less than or equal to the casewise breakdown value. Therefore, all upper bounds on casewise breakdown values in the literature also hold for cellwise breakdown values.

When proving breakdown values one often assumes that the original data set $\bm{X}$ is in general position, meaning that no more than $d$ points lie in any $d-1$ dimensional affine subspace. In particular, no three points lie on a line, no 4 points lie on a plane, and so on. When the data are drawn from a continuous distribution, it is in general position with probability 1. Real data have a limited precision, so they are not always in general position.

The inequalities in Proposition 1 can be strict. For instance, the casewise implosion breakdown value of the classical covariance matrix Cov at a dataset in general position is very high, in fact it is $(n-d)/n$ which goes to 1 for increasing sample size $n$. This is because whenever $d+1$ of the original data points are kept, Cov remains nonsingular. In stark contrast, its cellwise implosion breakdown value is quite low:

To see why, let us pick $d$ points of $\bm{X}$ which lie on a hyperplane that is not parallel to any coordinate axis. In the remaining $n-d$ rows we can then replace a single cell such that all of the resulting points lie on the same hyperplane, so Cov becomes singular. We can do this by replacing no more than $\lceil(n-d)/d\rceil$ cells in each variable, which is a fraction $\lceil(n-d)/d\rceil/n$ of its $n$ cells.

Raymaekers and Rousseeuw (2023) recently derived a similar upper bound for all affine equivariant estimators $\bm{\widehat{\Sigma}}$. In order to obtain a higher cellwise breakdown value we are thus forced to leave the realm of affine equivariance. In fact, the constraint $\lambda_{d}(\bm{\widehat{\Sigma}})\geqslant a>0$ in the definition (9) of cellMCD is not affine invariant, but it keeps $\bm{\widehat{\Sigma}}$ from imploding. Therefore the cellwise implosion breakdown value of cellMCD is 1.

We also want to know the breakdown value of its location estimate $\bm{\widehat{\mu}}$ and the explosion breakdown value of $\bm{\widehat{\Sigma}}$. These naturally depend on the choice of $h$.

Proposition 2 shows that cellMCD is highly robust. Its proof is in Section A.1 of the Supplementary Material. By Proposition 1, it follows that these lower bounds also hold for the casewise breakdown values. This also implies that the method works on a mix of cellwise and casewise outliers as well. We do not actually recommend to choose $h$ as low as the proposition allows: as explained before this could lead to some poorly defined covariances and numerical instability. We stick with our earlier recommendation of $h\geqslant 0.75n$, and in fact $h=0.75n$ is the default in our implementation.

Let us now turn to the asymptotic behavior of cellMCD. At the uncontaminated model distribution and for large $n$ only a small fraction of cells is actually discarded, due to our choice of the constants $q_{j}$ in the penalty term. In that situation the large-sample behavior of cellMCD is therefore the same as without the columnwise constraint on $\bm{W}$. The cellMCD objective can then be written as

where

in which $\bm{q}=(q_{1},\ldots,q_{d})$ and $w=(w_{1},\ldots,w_{d})$. The cellMCD estimate is then

with $F_{n}$ the empirical distribution and $\Theta$ the parameter space of $(\mu,\Sigma)$, which incorporates the condition $\lambda_{d}(\Sigma)\geqslant a$. Denote the set of minimizers as $\Theta^{*}$. In section A.2 of the Supplementary Material the following Wald-type consistency result is shown, using work of Van der Vaart (2000):

The population minimizer for $\Sigma$ is not quite the underlying parameter, since a small fraction of cells is always given weight zero due to the penalty term in the objective. But for the location $\mu$ we can prove that the unique minimizer is indeed the underlying parameter vector, so the cellMCD functional for location is Fisher consistent:

In the algorithm we will need the following result about decomposing the Mahalanobis distance and the likelihood.

The proof can be found in section A.3 in the Supplementary Material. The proposition can be interpreted as follows. Take a case $\bm{x}_{i}$ with some but not all cells missing, and for simplicity assume that its missing components come first. Then put $\bm{x}_{1}=\bm{x}_{i}^{(\bm{1}-\bm{w}_{i})}$ and $\bm{x}_{2}$ the remainder. If $(\bm{\mu},\bm{\Sigma})$ are the true underlying parameters, $\bm{\widehat{x}}_{1}$ is the conditional expectation $E[\bm{X}_{1}|\bm{X}_{2}=\bm{x}_{2}]$ and $\bm{C}_{1}$ is the conditional covariance matrix $\mbox{Cov}[\bm{X}_{1}|\bm{X}_{2}=\bm{x}_{2}]$. The additivity in (16) and (17) justifies the use of the partial Mahalanobis distances and the observed likelihood in our setting. Moreover, the fact that the difference of two ‘nested’ $\MD^{2}$ is again an $\MD^{2}$ and hence non-negative implies that the $\MD^{2}$ is monotone for nested sets of variables. In particular, if $\bm{x}$ is observed fully we can write

where each time $s^{2}$ is the matrix $\bm{C}_{1}$ (which is a scalar here) and the residuals are$r(x_{1}|x_{2},\ldots,x_{d})=x_{1}-\widehat{x}_{1}(x_{2},\ldots,x_{d})$ and so on. Note that (18) holds for any order of the $d$ variables. However, in each order the relative contribution of variable $j$ to the total $\MD^{2}(\bm{x},\bm{\mu},\bm{\Sigma})$ may be different. For the likelihood we obtain similarly

in which the terms do not need to be positive.

If we set $q_{j}=0$ in the objective function (9) of cellMCD and use casewise weights, i.e. casewise constant $w_{ij}$ , we recover the objective function (3) of the original casewise MCD. The latter is not convex in $\bm{\mu}$ and $\bm{\Sigma}$, so neither is (9). The crucial ingredient in the algorithm for the casewise MCD is the concentration step (C-step) of Rousseeuw and Van Driessen (1999). After each C-step the new objective value is less than or equal to the old objective value, so iterating C-steps always converges to a stationary point. We will now construct a C-step for cellMCD with the same properties. Let us denote the current solution of cellMCD by $\bm{\widehat{\mu}}^{(k)}$, $\bm{\widehat{\Sigma}}^{(k)}$, and $\bm{W}^{(k)}$. Then the new C-step proceeds as follows.

For the proof see section A.3 in the Supplementary Material, which also contains the pseudocode of the algorithm. Many variations of the C-step are possible, such as cycling through the columns of $\bm{\widetilde{W}}$ in a different order. We could also cycle through the columns of $\bm{\widetilde{W}}$ more than once in part (a), and/or run more than one EM-step in part (b). But experiments in section A.6 of the Supplementary Material show that these changes have a negligible and non-systematic effect on estimation accuracy, so we stay with the current version which is the fastest.

Note that cellMCD can still be used when the data contains missing cells, indicated by $u_{ij}$ which are 0 for missing cells and 1 elsewhere. In that situation we first have to remove variables with more than $n-h$ missing values. In the C-step it then suffices to force $w_{ij}=0$ whenever $u_{ij}=0$.

In order to start our C-steps we need an initial estimator. In our experiments we found that the DDCW estimator of Raymaekers and Rousseeuw (2021b) gives good results and is very fast. It is a combination of the DetectDeviatingCells (DDC) method of Rousseeuw and Van den Bossche (2018) and the fast correlation method in (Raymaekers and Rousseeuw, 2021a). DDCW is described in section A.4 of the Supplementary Material. Instead of starting from a single initial estimate, one could also start from several initial estimates. Iterating C-steps from each (with the same $q_{j}$ and $a>0$) until convergence, one can then keep the solution with the lowest objective (9).

The only remaining question is how to select the constants $q_{j}$ but this is quite simple, we do not need cross-validation or an information criterion. In (20) the term $(x_{ij}-\widehat{x}_{ij})^{2}/C_{ij}$ is the square of the residual $x_{ij}-\widehat{x}_{ij}$ standardized robustly. For inlying cells this should be below a cutoff, for which we take the chi-squared quantile $\chi^{2}_{1,p}$ with one degree of freedom and probability $p$. The term $\ln(C_{ij})$ is approximated by using the conditional variance of variable $j$ in the initial estimate $\bm{\widehat{\Sigma}}_{0}$ , given by $C_{j}:=1/(\bm{\widehat{\Sigma}}_{0}^{-1})_{jj}$ . So we set each $q_{j}$ equal to

The effect of this choice is that a cell $x_{ij}$ is flagged iff it lies outside a robust tolerance interval around its predicted value $\widehat{x}_{ij}$ with coverage probability $p$. Therefore we only have to choose a single cutoff probability $p$ to generate all $q_{j}$ automatically. From simulations and examples we found that $p=0.99$ was a good choice overall, so it is set as the default. Section A.5 provides more information on the $q_{j}$ and the choice of $p$.

The algorithm has been implemented as the R function cellMCD(). It starts by checking the data for non-numerical variables, cases with too many NA’s and so on. Next, it robustly standardizes the variables, and then computes the initial estimator followed by C-steps until convergence. The constraint $\lambda_{d}(\bm{\widehat{\Sigma}})\geqslant a$ is applied to the standardized data, with default $a=10^{-4}$. The function also reports the number of flagged cells in each variable. All the plots in the next section were made by the companion function plot_cellMCD(). Both functions have been included in the R package cellWise on CRAN.

We will illustrate cellMCD on the cars data obtained from the Top Gear website by Alfons (2016), focusing on the 11 numerical variables price, displacement, horsepower, torque, acceleration time, top speed, miles per gallon, weight, length, width, and height. This dataset is popular because both the variables and the cases (the cars) can easily be interpreted. After removing two cars with mostly NA’s we have $n=295$. We also replaced the highly right-skewed variables price, displacement, horsepower, torque, and top speed by their logarithms. On these data we ran cellMCD in its default version.

To visualize the results, we first look by variable. Consider variable $j$, say horsepower. Its $i$-th cell has observed value $x_{ij}$ as well as its prediction $\widehat{x}_{ij}$ obtained from the unflagged cells in the same row $i$, as in (20). In (20) we also see the conditional variance $C_{ij}$ of this cell. It is then natural to plot the standardized cellwise residual

which is NA when $x_{ij}$ is missing. The left panel of Figure 2 shows the standardized residuals of the variable horsepower versus the index (case number) $i$. This plot was made by the function plot.cellMCD(), which also draws a horizontal tolerance band given by $\pm\,c$ where $c=\sqrt{\chi^{2}_{1,0.99}}\approx 2.57$ . Here, some residuals stick out below the tolerance band. The Renault Twizy and Citroen DS3 are energy savers, whereas the Caterham is a super lightweight fun car. The most extreme outlier is the Chevrolet Volt with a standardized residual below $-8$. Top Gear lists this car’s power as 86 hp, which cellMCD says is very low compared to what would be expected from the other 10 characteristics of this car. Looking it up revealed that the Volt actually has 149 hp. As far as we know this data error was not detected before.

The right panel of Figure 2 plots the standardized residuals of the variable length versus the observed length itself. The vertical lines are at $T\pm\,cS$ where $T$ and $S$ are robust univariate location and scale estimates of length, obtained from the function estLocScale() in the R package cellWise. The points to the left and right of such a vertical tolerance band are marginally outlying, i.e. their length stands out by itself without regard to the other variables. In the bottom left region of the plot we see five cars that are marginal outliers to the left and at the same time have outlying negative residuals, so they are short in absolute terms, as well as relative to what would be expected from their other characteristics. The Smart fortwo, Renault Twizy, Toyota IQ and Aston Martin Cygnet are indeed tiny.

However, not all cellwise outliers are marginal outliers. In the middle bottom part of the plot we marked three cars whose length is not unusual by itself, but that are short relative to what would be expected based on their other 10 variables. They are sports cars, often built small to achieve high speeds. Note that there could also be points that lie inside the horizontal band but (slightly) outside the vertical band. They would correspond to cells that look a bit unusual in the variable $j$, but whose observed value $x_{ij}$ is not that far from the predicted $\widehat{x}_{ij}$ based on its other variables.

The left panel of Figure 3 plots the standardized residual of each car’s weight versus its prediction. Since all the points lie within the vertical tolerance band, no predictions are outlying. But we do see some outlying residuals, most of which can easily be explained. The Bentley is a heavy luxury car, and the Mercedes-Benz G an all-terrain vehicle. Below the horizontal tolerance band we see four lightweight sports cars. What remains is the Peugeot 107 which is small but not sporty at all. Top Gear reports its weight as 210 kg, which seems much too light for a car. Based on its other characteristics, cellMCD predicts its weight as 757 kg with a standard error of 89.5 kg. Looking up this car, its actual weight turns out to be 800 kg, so the value in the Top Gear dataset was mistaken.

The right panel of Figure 3 shows the observed value of top speed versus its prediction. Below the superimposed $y=x$ line we find some electric cars (BMW i3, Vauxhall Ampera) and some small cars (Smart fortwo and Renault Zoe). The one standing out most is the Renault Twizy, a tiny electric one-seater vehicle. Above the line we see some extremely fast sports cars. Also note that some points appear to lie on a horizontal line. Top Gear reports their top speed as 155 mph, corresponding to 250 km/hour. Many of these cars were produced by Audi, BMW and Mercedes with a built-in 250 km/hour speed limiter.

The four plot types in Figures 2 and 3 all focus on a single variable. It can also be instructive to look at a pair of variables, say $j$ and $k$. Figure 4 shows the variables width versus acceleration. The points for which $w_{ij}=0$ or $w_{ik}=0$ or both are automatically plotted in red. The figure also contains an ellipse, given by

where $q$ is the 0.99 quantile of the $\chi_{2}^{2}$ distribution with two degrees of freedom. Note that outlyingness in this type of plot differs from cellwise outlyingness, since the former refers to two variables only, whereas the latter uses all 11 variables. So it is not unusual to see some red points inside the ellipse, and some black points outside it.

The width of the Land Rover is flagged as this is a wide all terrain vehicle. The red vertical line connects the observed point $(x_{ij},x_{ik})$ to its predicted point $(\widehat{x}_{ij},\widehat{x}_{ik})$ plotted in blue. That the line is vertical means that the width cell was flagged whereas the acceleration cell was not, that is, $w_{ik}=0$ and $w_{ij}=1$ . The acceleration of the Ssangyong Rodius and Lotus Elise is outlying on the left. In fact, Top Gear lists their acceleration time as 0 which is physically impossible: presumably the true value was missing and encoded as 0 instead of NA. The same happens for the Renault Twizy. Note that also the width cell of the Twizy is flagged, so the red line to its predicted point is slanted instead of horizontal. The Caterham also has both cells flagged, as seen from its slanted line.

In this section we evaluate the performance of cellMCD by a simulation study. The clean data is generated as $n$ points from a $d$-variate Gaussian distribution with mean $\bm{\mu}=\bm{0}$. Since there is no affine equivariance, letting $\bm{\Sigma}$ be the identity matrix is not sufficient. Instead we use the types “A09” and “ALYZ”. The entries of the A09 correlation matrix are given by $\bm{\Sigma}_{ij}=0.9^{|i-j|}$, yielding both small and large correlations. The ALYZ type are randomly generated correlation matrices following the procedure of Agostinelli et al. (2015) and typically have mostly small absolute correlations. We consider three combinations of sample size and dimension $(n,d)$: $(100,10)$, $(400,20)$, and $(800,40)$.

In these clean data, we then replace a fraction $\varepsilon$ in $\{0.1,0.2\}$ of cells by contaminated cells. These are generated as follows. First, for each column in the data matrix we randomly sample $n\varepsilon$ indices of cells to be contaminated. In each row, say $(z_{1},\ldots,z_{d})$, we then collect the indices of the cells to be contaminated. Denote this set of size $k$ by $K=\{j_{1},\ldots,j_{k}\}$. We next replace the cells $(z_{j_{1}},\ldots,z_{j_{k}})$ by the $k$-dimensional vector $\gamma\sqrt{k}\,\bm{v}_{K}/\mbox{MD}(\bm{v}_{K},\bm{\mu}_{K},\bm{\Sigma}_{K})$ where $\bm{\mu}_{K}$ and $\bm{\Sigma}_{K}$ are $\bm{\mu}$ and $\bm{\Sigma}$ restricted to the indices in $K$. The scalar $\gamma>0$ quantifies the distance of the outlying cells to the center of the distribution, and we vary $\gamma$ over $1,\ldots,10$. The vector $\bm{v}_{K}$ is the normed eigenvector of $\bm{\Sigma}_{K}$ with the smallest eigenvalue. In each row, the outlying cells are thus structurally outlying in the subspace generated by the variables in $K$. Therefore, these cells will often not be marginally outlying, especially when $|K|$ is large and $\gamma$ is relatively small, which makes them hard to detect. The R-package cellWise (Raymaekers and Rousseeuw, 2022) contains the function generateData which generates the contaminated data according to this procedure.

We compare the proposed method cellMCD to the following alternative estimators:

• •

  Grank, Spearman: the Gaussian and Spearman rank-based estimators used in Öllerer and Croux (2015) and Croux and Öllerer (2016);

• •

  GKnpd: the Gnanadesikan-Kettenring estimator used in Tarr et al. (2016);

• •

  2SGS: the two-step generalized S-estimator of Agostinelli et al. (2015);

• •

  DI: the detection-imputation algorithm of Raymaekers and Rousseeuw (2021b).

In order to evaluate the performance of the different estimators, we compute the Kullback-Leibler discrepancy between the estimated $\widehat{\bm{\Sigma}}$ and the true $\bm{\Sigma}$ given by

For each setting of the simulation parameters we generate 100 random datasets, and average the Kullback-Leibler discrepancy over these 100 replications. (For the variability around these averages see subsection A.6.1.)

Figure 5 presents the results for $d=10$, $n=100$ and $\varepsilon=0.1$. (The results for $\varepsilon=0.2$ were similar.) Both cellMCD and DI perform well, as does 2SGS provided $\gamma\geqslant 4$. As expected, the classical covariance matrix (Cov) and the casewise MCD (labeled caseMCD) were not robust to these adversarial cellwise outliers. Note that the performances of Grank, Spearman and GKnpd do not improve as $\gamma$ increases. While these estimators bound the influence that a single cell can have on the estimation, the effect remains substantial as the cell becomes more outlying. This is in contrast to 2SGS, DI and cellMCD in which far outliers get a zero weight.

The top panels of Figure 6 show the results for $n=400$ and $d=20$. The relative performances are similar to Figure 5. The 2SGS method still does well when $\gamma>4$, but now suffers more for low $\gamma$. The performances of DI and cellMCD are again very close, with cellMCD often doing slightly better.

The lower panels with $n=800$ and $d=40$ are similar, with cellMCD performing best for all values of $\gamma$ while DI is quite close, and 2SGS only doing well for higher $\gamma$.

Table 1 lists the computation times of the methods in the simulation, in seconds. The first five methods are fast but they performed poorly. The bottom three methods did better. In dimensions 20 and 40 the cellMCD method was the fastest among them.

We are also interested in the performance of these methods on data without outliers. For this we repeated the simulation with $\varepsilon=0$, again with 100 replications. The variability of each entry of the covariance matrix was measured taking the Fisher information of that entry into account. These results were then averaged over the upper triangular matrix entries including the diagonal. Next we divided the MSE of the classical MLE estimator by that of each robust method, yielding the finite-sample efficiencies in Table 2.

We see that the efficiency of cellMCD averages over 90%, which is excellent for a highly robust covariance estimator. This is similar to 2SGS, and outperforms DI. As expected Grank has a high efficiency, but we just saw that it performed poorly under contamination, as did GKnpd and Spearman. The finite-sample efficiency of cellMCD is much higher than that of the casewise MCD with the same coverage parameter $h=0.75n$, which is under 0.70 for this range of dimensions $d$. This is due to the penalty term in (9), which made the number of actually discarded cells much smaller than $0.25\,n$.

We conclude that cellMCD is about equally robust as DI but with better efficiency, and is about as efficient as 2SGS but with better robustness at contaminated data. Moreover, it does substantially better at contaminated data than the remaining methods.

The cellMCD method proposed here has an elegant formulation based on a single objective function, making it easier to understand than the earlier 2SGS and DI methods. We proved its good breakdown properties and consistency, and like the casewise MCD it can be computed by an algorithm based on C-steps that always lower the objective function and is guaranteed to converge. We have illustrated cellMCD on a real data set where the accompanying graphical displays revealed interesting aspects of the data that aided interpretation. Simulations indicate that cellMCD outperforms earlier cellwise methods, while being conceptually simple and rather fast to compute.

CellMCD is cellwise robust and incorporates a kind of sparsity penalty (on $\bm{1}-\bm{W}$). This naturally brings to mind the work of Candès et al. (2011). The goals are clearly related, but there are also some differences. The first is that their work assumes that the cellwise outlier pattern $\bm{W}$ is drawn uniformly at random, whereas we adopt the robustness paradigm that the outliers may be placed adversarially. Secondly, the method of Candès et al. (2011) is equivariant for transposing the data matrix, so it treats cases and variables in the same way, whereas in our setting they have to be treated differently. We do allow for some rows being flagged entirely, whereas we cannot allow flagging an entire column as this would make $\bm{\mu}$ and $\bm{\Sigma}$ not identifiable, which motivates our constraint $||\bm{W}_{.j}||_{0}\geqslant h$ for $j=1,\ldots,d$ .

The fact that implosion breakdown can happen easily in the cellwise setting, see (13), was not mentioned in the literature before. We feel that, apart from cellMCD, also other cellwise robust covariance estimators could benefit from a constraint such as $\lambda_{d}(\bm{\widehat{\Sigma}})\geqslant a$, or similarly from a formulation in which $\bm{\widehat{\Sigma}}$ is a convex combination of two matrices, one of which is a small multiple of the identity matrix.

The casewise MCD is typically followed by a reweighting step. This works as follows. First, the estimated covariance matrix $\bm{\widehat{\Sigma}}$ is multiplied by a correction factor $c_{n,d,h}$ such that $c_{n,d,h}\bm{\widehat{\Sigma}}$  is roughly unbiased when the original data are generated from a Gaussian distribution. Next, one computes the squared robust distances of the data points, given by $\RD^{2}_{i}=(\bm{x}_{i}-\bm{\widehat{\mu}})^{\top}(c_{n,d,h}\bm{\widehat{\Sigma}})^{-1}(\bm{x}_{i}-\bm{\widehat{\mu}})$. Each case $\bm{x}_{i}$ then gets a weight $w_{i}$ depending on its $\RD^{2}_{i}$ . Typically, the weight is set to 1 when $\RD_{i}^{2}$ is below some quantile of the $\chi_{d}^{2}$ distribution with $d$ degrees of freedom, and to 0 otherwise. The final estimates are then the weighted mean and the weighted covariance matrix (4). This reweighting step increases the finite-sample efficiency of the estimator.

For cellMCD, the analogous reweighting step would compute the standardized residual (22) of every cell $x_{ij}$ and compare its square to a quantile of the $\chi_{1}^{2}$ distribution with 1 degree of freedom, yielding zero-one weights $w_{ij}$. With these $w_{ij}$ one would then run the EM algorithm on the original data. But in fact, the result is not very different from the cellMCD result. This is because all the ingredients are already used in cellMCD, which contains the squared standardized residual in (20), the $\chi_{1}^{2}$ quantile in (21), and the partial likelihood on which EM is based in (9). So in some sense the components of a reweighting step are already built into cellMCD itself. This explains its rather high finite-sample efficiency in Table 2.