1. Introduction

Most datasets come in the form of a rectangular matrix $\mathbf{X}$ with n rows and d columns, where n is called the sample size and d the dimension. The rows of $\mathbf{X}$ correspond to the cases, whereas the columns are the variables. Both n and d can be large, and it may happen that d > n. There exist many data models as well as techniques to fit them, such as regression and principal component analysis.

It is well-known that many datasets contain outliers. Depending on the circumstances, outliers may be (a) undesirable errors, which can adversely affect the data analysis, or (b) valuable nuggets of unexpected information. Either way, it is important to be able to detect the outliers, which can be hard for high d.

In statistics and data analysis, the word outlier typically refers to a row of the data matrix, as in the left panel of Figure 1. There has been much research since the 1960s to develop fitting methods that are less sensitive to such outlying rows, and that can detect the outliers by their large residual (or distance) from that fit. This topic goes by the name of robust statistics, see, for example, the books by Maronna, Martin, and Yohai (2006) and Rousseeuw and Leroy (1987).

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 1. The rowwise and cellwise outlier paradigms (black means outlying).

Display full size

Recently, researchers have come to realize that the outlying rows paradigm is no longer sufficient for modern high-dimensional datasets. It often happens that most data cells (entries) in a row are regular and just a few of them are anomalous. The first article to formulate the cellwise paradigm was Alqallaf et al. (2009). They noted how outliers propagate: given a fraction ϵ of contaminated cells at random positions, the expected fraction of contaminated rows is (1) $$1-{(1-\epsilon )}^d$$(1) which quickly exceeds 50% for increasing ϵ and/or increasing dimension d, as illustrated in the right panel of Figure 1. This is fatal because rowwise methods cannot handle more than 50% of contaminated rows if one assumes some basic invariance properties, see, for example, Lopuhaä and Rousseeuw (1991).

The two paradigms are quite different. The outlying row paradigm is about cases that do not belong in the dataset, for instance, because they are members of a different population. Its basic units are the cases, and if you interpret them as points in d-dimensional space they are indeed indivisible. But if you consider a case as a row in a matrix then it can be divided, into cells. The cellwise paradigm assumes that some of these cells deviate from the values they should have had, perhaps due to gross measurement errors, whereas the remaining cells in the same row still contain useful information.

The anomalous data cells problem has proved to be quite hard, as the existing tools do not suffice. Recent progress was made by Danilov (2010),  Van Aelst, Vandervieren, and Willems (2012), Agostinelli et al. (2015),  Öllerer, Alfons, and Croux (2016), and Leung, Zhang, and Zamar (2016).

Note that in the right panel of Figure 1 we do not know at the outset which (if any) cells are black (outlying), in stark contrast with the missing data framework. To illustrate the difficulty of detecting deviating data cells, let us look at the artificial bivariate (d = 2) example in Figure 2. Case 1 lies within the near-linear pattern of the majority, whereas cases 2, 3, and 4 are outliers. The first coordinate of case 2 lies among the majority of the x_i1 (depicted as tickmarks under the horizontal axis) while its second coordinate x₂₂ is an outlying cell. For case 3, the first cell x₃₁ is outlying. But what about case 4? Both of its coordinates fall in the appropriate ranges. Either coordinate of case 4 could be responsible for the point standing out, so the problem is not identifiable. We would have to flag the entire row 4 as outlying.

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 2. Illustration of bivariate outliers.

Display full size

But now suppose that we obtain five more variables (columns) for these cases, which are correlated with the columns we already have. Then it may be possible to see that x₄₁ perfectly agrees with the values in the new cells of row 4, whereas x₄₂ does not. Then we would conclude that the cell x₄₂ is the culprit. Therefore, this is one of the rare situations where a higher dimensionality is not a curse but can be turned into an advantage. In fact, the method proposed in the next section typically performs better when the number of dimensions increases.

This example shows that deviating data cells do not need to stand out in their column: it suffices that they disobey the pattern of relations between this variable and those correlated with it. Here, we propose the first method for detecting deviating data cells that takes the correlations between the variables into account. Unlike existing methods, our DetectDeviatingCells (DDC) method is not restricted to data with over 50% of clean rows. Other advantages are that it can deal with high dimensions and that it provides predicted values of the outlying cells. As a byproduct it imputes missing values in the data. For software availability, see Section 8.

2. Brief Sketch of the Method

The (possibly idealized) model states that the rows ${\mathbf{x}}_i$ were generated from a multivariate Gaussian distribution with unknown d-variate mean $\mu$ and positive semidefinite covariance matrix $\Sigma$, after which some cells were corrupted or became missing. However, the algorithm will still run if the underlying distribution is not multivariate Gaussian.

The method uses various devices from robust statistics and is described more fully in Section 4. We only sketch the main ideas here. The method starts by standardizing the data and flagging the cells that stand out in their column. Next, each data cell is predicted based on the unflagged cells in the same row whose column is correlated with the column in question. Finally, a cell for which the observed value differs much from its predicted value is considered anomalous. The method then produces an imputed data matrix, which is like the original one except that cellwise outliers and missing values are replaced by their predicted values. The method can also flag an entire row if it contains too much anomalous behavior, so that we have the option of downweighting or dropping that row in subsequent fits to the data.

This^{Figure 3}method looks unusual but it successfully navigates around the many pitfalls inherent in the problem, and worked the best out of the many approaches we tried. Its main feature is the prediction of individual cells. This can be characterized as a locally linear fit, where “locally” is not meant in the usual sense (of Euclidean distance) but instead refers to the space of variables endowed with a kind of correlation-based metric.

Along the way DDC imputes missing data values by their predicted values. This is far less efficient than the EM algorithm of Dempster, Laird, and Rubin (1977) when the data are Gaussian and outlier-free, but it is more robust against cellwise outliers.

The DDC method has the natural equivariance properties. If we add a constant to all the values in a column of X, or multiply any column by a nonzero factor, or reorder the rows or the columns, the result will change as expected.

3. Examples

The first dataset was scraped from the website of the British television show Top Gear by Alfons (2016). It contains 11 objectively measured numerical variables about 297 cars. Five of these variables (such as Price) were rather skewed so we logarithmically transformed them.

The left panel of Figure 3 shows the result of applying a simple univariate outlier identifier to each of the columns of this dataset separately. For column j, this method first computes a robust location m_j (such as the median) of its values, as well as a robust scale s_j and then computes the standardized values z_ij = (x_ij − m_j)/s_j. A cell is then flagged as outlying if |z_ij| > c where the cutoff c = 2.576 corresponds to a tolerance of 99% so that under ideal circumstances about 1% of the cells would be flagged. Most of the cells were yellow, indicating they were not flagged. Here, we only show some of the more interesting rows out of the 297. Deviating cells whose z_ij is positive are colored red, and those with negative z_ij are in blue. Missing data are shown in white and labeled NA (from “Not Available”).

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 3. Selected rows from cell maps of Top Gear data: (left) when detecting anomalous data cells per column; (right) when using the proposed method.

Display full size

The right-hand panel of Figure 3 shows the corresponding cell map obtained by running DetectDeviatingCells on these data. In the algorithm, we used the same 99% tolerance setting yielding the same cutoff c = 2.576. The color of the deviating cells now reflects whether the observed cell value is much higher or much lower than the predicted cell value.

The new method shows much more structure than looking by column only. The high gas mileage (MPG) of the BMW i3 is no longer the only thing that stands out, and in fact it is an electric vehicle with a small additional gasoline engine. The columnwise method does not flag the Corvette’s displacement, which is not that unusual by itself, but it is high in relation to the car’s other characteristics. None of the properties of the Land Rover Defender (an all-terrain vehicle) stand out on their own, but their combination does. The weight of 210 kg listed for the Peugeot 107 is clearly an error, which gets picked up by both methods. The univariate method only flags the height of the Ssangyong Rodius, whereas DetectDeviatingCells also notices that its acceleration time of zero seconds from 0 to 62 mph is too low compared to its other properties, and in fact it is physically impossible.

As a second example we take the Philips dataset, which measures d = 9 characteristics of n = 677 TV parts created by a new production line. The left panel of Figure 4 is a rotated version of Figure 9(a) in Rousseeuw and Van Driessen (1999) and shows, from top to bottom, the robust distance ${\text{RD}}_i$ for i = 1, …, n. These were computed as $${\text{RD}}_i=\sqrt{{({\mathbf{x}}_i-\mathbf{T})}^{\text{'}}{\mathbf{S}}^{-1}({\mathbf{x}}_i-\mathbf{T})},$$

where $\mathbf{T}$ and $\mathbf{S}$ are the estimates of location and scatter obtained by the minimum covariance determinant (MCD) from (Rousseeuw 1985), a rowwise robust method. It shows that the process was out of control in the beginning (roughly rows 1–90) and near the end (around rows 480–550). However, by itself this does not yet tell us what happened in those products.

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 4. Left: robust distances of the rows of the Philips data. Right: map of DetectDeviatingCells on these data, combined per blocks of 15 rows.

Display full size

In the right panel, we see the cell map of DDC. The analysis was carried out on all 677 rows, but to display the results we created blocks of 15 rows and 1 column. The color of each cell is the “average color” over the block. The resulting cell map indicates that in the beginning variable 6 had lower values than would be predicted from the remaining variables, whereas the later outliers were linked with lower values of variable 8. (In between, the cell map gives hints about some more isolated deviations.) We do not claim that this is the complete answer, but at least it gives an initial indication of what to look for. (DDC also flagged some rows, but this is not shown in the cell map to avoid confusion.) Note that these data are not literally a multivariate sample because the rows were provided consecutively, but both the MCD and DDC methods ignore the time order of the rows.

The next example is the mortality by age for males in France, from 1816 to 2010 as obtained from the Human Mortality Database (2015). An earlier version was analyzed by Hyndman and Shang (2010). Each case corresponds to the mortalities in a given year. The left panel in Figure 5 was obtained by ROBPCA (Hubert, Rousseeuw, and Vanden Branden 2005), a rowwise robust method for principal component analysis. Such a method can only detect entire rows, and the rows it has flagged are represented as black horizontal lines. A black row does not reveal any information about its cells. The analysis was carried out on the dataset with the individual years and the individual ages, but as this resolution would be too high to fit on the page, we have combined the cells into 5 × 5 blocks afterward. The combination of some black rows with some yellow ones has led to gray blocks. We can see that there were outlying rows in the early years, the most recent years, and during two periods in between.

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 5. Male mortality in France in 1816–2010: (left) detecting outlying rows by a robust PCA method; (right) detecting outlying cells by DetectDeviatingCells. After the analysis, the cells were grouped in blocks of 5 × 5 for visibility.

Display full size

By contrast, the right panel in Figure 5 identifies a lot more information. The outlying early years saw a high infant mortality. During the Prussian war and both world wars there was a higher mortality among young adult men. And in recent years mortality among middle-aged and older men has decreased substantially, perhaps due to medical advances.

Our final example consists of spectra with d = 750 wavelengths collected on n = 180 archeological glass samples (Lemberge et al. 2000). Here, the number of dimensions (columns) exceeds the number of cases (rows). DetectDeviatingCells has no problems with this, and the R implementation took about 1 min on a laptop with Intel i7-4800MQ 2.7GHz processor for this dataset. The top panel in Figure 6 shows the rows detected by the robust principal component method. The lower panel is the cell map obtained by DetectDeviatingCells on this 180 × 750 dataset. After the analysis, the cells were again grouped in 5 × 5 blocks. We now see clearly which parts of each spectrum are higher/lower than predicted. The wavelengths of these deviating cells reveal the chemical elements responsible.

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 6. Cell maps of n = 180 archeological glass spectra with d = 750 wavelengths. The positions of the deviating data cells reveal the chemical contaminants.

Display full size

In general, the result of DetectDeviatingCells can be used in several ways:

In 2 and 3, one has the option to drop the flagged rows before taking the step. If that step is carried out by a sparse method such as the Lasso (Tibshirani 1996; Städler, Stekhoven, and Bühlmann 2014) or another form of variable selection, one would of course look more closely at the deviating cells in the variables that were selected.

Remark. In the above examples, the cellwise outliers happened to be mainly concentrated in fewer than half of the rows, allowing the rowwise robust methods MCD and ROBPCA to detect most of those rows. This may give the false impression that it would suffice to apply a rowwise robust method and then to analyze the flagged rows further. But this is not true in general: by (1) there will typically be too many contaminated rows, so rowwise robust methods will fail.

4. Detailed Description of the Algorithm

The DetectDeviatingCells (DDC) algorithm has been implemented in R and Matlab (see Section 8). As described before, the underlying model is that the data are generated from a multivariate Gaussian distribution $N(\mu ,\Sigma )$ but afterward some cells were corrupted or omitted. The variables (columns) of the data should thus be numerical and take on more than a few values. Therefore, the code does some preprocessing: it identifies the columns that do not satisfy this condition, such as binary dummy variables, and performs its computations on the remaining columns.

These remaining variables should then be approximately Gaussian in their center, that is, apart from any outlying values. This could be verified by QQ plots, and one could flag highly skewed or long-tailed variables. It is recommended to transform very non-Gaussian variables to approximate Gaussianity, like we took the logarithm of the right-skewed variable Price in the cars data in Section 3. More general tools are the Box–Cox and Yeo–Johnson (Yeo and Johnson 2000) transformations that have parameters, which need to be estimated as done by Bini and Bertacci (2006) and Riani (2008). We chose not to automate this step because we feel it is best carried out by a person with subject matter knowledge about the data. From here onward we will assume that the variables are approximately Gaussian in their center. The algorithm does not require joint normality of variables to run.

The actual algorithm then proceeds in eight steps.

Step 1: standardization. For each column j of $\mathbf{X}$, we estimate (2) $$m_j={\text{robLoc}}_i\left(x_{ij}\right)\text{and}{s}_j={\text{robScale}}_i(x_{ij}-m_j),$$(2) where robLoc is a robust estimator of location, and robScale is a robust estimator of scale, which assumes its argument has already been centered. The actual definitions of these estimators are given in the Appendix. Next, we standardize $\mathbf{X}$ to $\mathbf{Z}$ by (3) $$z_{ij}=(x_{ij}-m_j)/s_j.$$(3)

Step 2: apply univariate outlier detection to all variables. After the columnwise standardization in (3), we define a new matrix $\mathbf{U}$ with entries (4) $$u_{ij}=\left\{\begin{array}{cc}{z}_{ij}\hfill & \text{if}|z_{ij}|⩽c\hfill \\ \text{NA}\hfill & \text{if}|z_{ij}|>c.\hfill \end{array}\right.$$(4) Due to the standardization in (3), formula (4) is a columnwise outlier detector. The cutoff value c is taken as (5) $$c=\sqrt{{\chi }_{1,p}^2},$$(5) where χ²_{1, p} is the pth quantile of the chi-squared distribution with 1 degree of freedom, where the probability p is 99% by default so that under ideal circumstances only 1% of the entries gets flagged.

Step 3: bivariate relations. For any two variables h ≠ j, we compute their correlation as (6) $${\text{cor}}_{jh}={\text{robCorr}}_i(u_{ij},u_{ih}),$$(6) where robCorr is a robust correlation measure given in the Appendix (the computation is over all i for which neither u_ij or u_ih is NA). From here onward we will only use the relation between variables j and h when (7) $$|{\text{cor}}_{jh}|⩾\text{corrlim}$$(7) in which corrlim is set to 0.5 by default. Variables j that satisfy (7) for some h ≠ j will be called connected, and contain useful information about each other. The others are called standalone variables. For the pairs (j, h) satisfying (7) we also compute (8) $$b_{jh}={\text{robSlope}}_i\left(u_{ij}\right|u_{ih}),$$(8) where robSlope computes the slope of a robust regression line without intercept that predicts variable j from variable h (see the Appendix). These slopes will be used in the next step to provide predictions for connected variables.

Step 4: predicted values. Next we compute predicted values ${\widehat{z}}_{ij}$ for all cells. For each variable j we consider the set H_j consisting of all variables h satisfying (7), including j itself. For all i = 1, …, n, we then set (9) $${\widehat{z}}_{ij}=G\left(\{b_{jh}{u}_{ih};h\text{in}{H}_j\}\right),$$(9) where G is a combination rule applied to these numbers, which omits the NA values and is zero when no values remain. (The latter corresponds to predicting the unstandardized cell by the estimated location of its column, which is a fail-safe.) Our current preference for G is a weighted mean with weights w_jh = |cor_jh| but other choices are possible, such as a weighted median. The advantage of (9) is that the contribution of an outlying cell z_ih to ${\widehat{z}}_{ij}$ is limited because |u_ih| ⩽ c and it can only affect a single term, which would not be true for a prediction from a least-square multiple regression of z_.j on the d − 1 remaining variables z_.h together. We will illustrate the accuracy of the predicted cells by simulation in Section 5.

Step 5: deshrinkage. Note that a prediction such as (9) tends to shrink the scale of the entries, which is undesirable. We could try to shrink less in the individual terms b_jhu_ih but this would not suffice because these terms can have different signs for different h. Therefore, we propose to deal with the shrinkage after applying the combination rule (9). For this purpose, we replace ${\widehat{z}}_{ij}$ by $a_j{\widehat{z}}_{ij}$ for all i and j, where (10) $$a_j:={\text{robSlope}}_{i^{\text{'}}}\left(z_{i^{\text{'}}j}\right|{\widehat{z}}_{i^{\text{'}}j})$$(10) comes from regressing the observed z_.j on the (shrunk) predicted ${\widehat{z}}_{.j}$.

Step 6: flagging cellwise outliers. In Steps 4 and 5, we have computed predicted values ${\widehat{z}}_{ij}$ for all cells. Next we compute the standardized cell residuals (11) $$r_{ij}=\frac{z_{ij}-{\widehat{z}}_{ij}}{{\text{robScale}}_{i^{\text{'}}}(z_{i^{\text{'}}j}-{\widehat{z}}_{i^{\text{'}}j})}.$$(11) In each column j we then flag all cells with |r_ij| > c as anomalous, where c was given in (5). We also assemble the “imputed” matrix ${\mathbf{Z}}_{\text{imp}}$, which equals $\mathbf{Z}$ except that it replaces deviating cells z_ij and NA’s by their predicted values ${\widehat{z}}_{ij}$. The unflagged cells remain as they were.

We have the option of setting the flagged cells to NA and repeating steps 4 to 6 to improve the accuracy of the estimates. This can be iterated.

Step 7: flagging rowwise outliers. To decide whether to flag row i, we could just count the number of cells whose |r_ij| exceeds a cutoff a, but this would miss rows with many fairly large |r_ij| < a. The other extreme would be to compare ${\text{ave}}_j\left(r_{ij}^2\right)$ to a cutoff, but then a row with one very outlying cell would be flagged as outlying, which would defeat the purpose. We do something in between. Under the null hypothesis of multivariate Gaussian data without any outliers, the distribution of the r_ij is close to standard Gaussian, so the cdf of r²_ij is approximately the cdf F of χ²₁. This leads us to the criterion (12) $$T_i=\underset{j=1}{\overset{d}{\text{ave}}}F\left(r_{ij}^2\right).$$(12) We then standardize the T_i as in (3) and flag the rows i for which the standardized T_i exceed the cutoff c of (5).

When we find that row i has an unusually large T_i this does not necessarily “prove” that it corresponds to a member of a different population, but at least it is worth looking into.

Even though the T_i can flag some rowwise outliers, there are types of rowwise outliers that it may not detect, for instance, in the barrow wheel configuration of Stahel and Mächler (2009) where a rowwise outlier may not have any large r²_ij. Therefore, it is recommended to use rowwise robust methods in subsequent analyses of the data, after dropping the rows that were flagged.

Step 8: destandardize. Next, we turn the imputed matrix ${\mathbf{Z}}_{\text{imp}}$ into an imputed matrix ${\mathbf{X}}_{\text{imp}}$ by undoing the standardization in (3). The main output of DDC is ${\mathbf{X}}_{\text{imp}}$ together with the list of cells flagged as outlying and, optionally, the list of rowwise outliers that were found.

Predicting the value of a cell in Steps 4 and 5 above evokes the imputation of a missing value, as in the EM algorithm of Dempster, Laird, and Rubin (1977). However, there are two reasons why the EM method may not work in the present situation. First of all, we would need estimators of location and scatter $\widehat{\mu }$ and $\widehat{\Sigma }$ that are robust against cellwise and rowwise outliers, which we only know how to obtain after detecting those outliers. And second, the conditional expectation of x_ij is a linear combination of all the available x_ih with h ≠ j but some of those cells may be outlying themselves, thus spoiling the sum.

The DDC method employs measures of bivariate correlation, and regression for bivariate data. This may seem naive compared to estimating scatter matrices and/or performing multiple regression on sets of q > 2 variables. However, the latter would imply that we must compute scatter matrices between q variables in an environment where a fraction ϵ of cells is contaminated, so that according to (1) a fraction 1 − (1 − ϵ)^q of the rows of length q is expected to be contaminated. The latter fraction grows very quickly with q, and once it exceeds 50% there is no hope of estimating the scatter matrix by a rowwise robust method. So the larger q is, the less this approach can work. Therefore, we chose the smallest possible value q = 2, which has the additional advantage that the computational effort to robustly estimate bivariate scatter is much lower than for higher q.

To illustrate this, Table 1 shows the effect of fitting substructures in dimension q < d for various values of q. The total time complexity is the product of the time needed for a 50%-rowwise-breakdown fit (regression or scatter) to a q-variate dataset (in which elemental sets are formed by q − 1 data points and the origin) and the number of combinations of q out of d variables (for large d). The complexity grows very fast with q. The next column is the expected breakdown value, that is, what fraction of cellwise outliers at random positions these fits can typically withstand, that is, 1 − 2^{− 1/q}. The expected breakdown value drops quickly with q. The final column is the probability that a subrow of q entries contains no outliers when the probability of a cell being outlying is 10%. Table 1 is in fact overly optimistic because it does not account for the problem that the predictions from such q-variate fits are affected whenever at least one of the cells in it is outlying.

This should not be confused with the approach of Kriegel et al. (2009) who look for rowwise outliers under the assumption that outlyingness of a row is due to only a few variables whereas the other variables are deemed irrelevant. That situation is different from both the general rowwise outlier setting and the cellwise outlier model, in each of which all variables may be relevant.

As^{Figure 7}described, the computation time (number of operations) of DDC is on the order of O(nd²) due to computing all d(d − 1)/2 correlations between variables, where each correlation takes O(n) time. The required memory (storage) space is then also O(nd²). However, if the dimension d is over (say) 1000, we switch to predicting each column by the k columns that are most correlated with it, thereby bringing down the space requirement to O(nd), which cannot be reduced further as it is proportional to the size of the data matrix itself. On a parallel architecture, we can reduce the computation time by letting each processor compute the k correlations for a subset of columns to be predicted. On a nonparallel system, we could switch to a fast approximate k-nearest neighbor method such as the celebrated algorithm of Arya et al. (1999), which would lower the computation time from O(nd²) to $O\left(nd\text{log}\right(d\left)\right)$. Compared to analyzing each column separately, this algorithm only spends an additional time factor $\text{log}\left(d\right),$ which grows very slowly with d. A related speedup is used in Google Correlate (Vanderkam et al. 2013).

5. Comparison to Other Detection Methods

We now compare DDC to a method that detects outlying cells per column, the univariate GY filter of Gervini and Yohai (2002) as described in (Agostinelli et al. 2015). Applying the univariate GY filter to the columns of the cars data actually gives the same result as the simpler columnwise detection in Figure 3. The key difference between these methods is that the GY filter uses an adaptive cutoff rather than the fixed cutoff (5).

We will only report a small part of our simulations, the results of other settings being similar. First, clean multivariate Gaussian data were generated with $\mu =\mathbf{0}$ and two types of covariance matrices $\Sigma$ with unit diagonal. The ALYZ (Agostinelli et al. 2015) random correlation matrices yield relatively low correlations between the variables, whereas the A09 correlation matrix given by ρ_jh = ( − 0.9)^{|h − j|} yields both high and low correlations.

Next, these clean data were contaminated. Outlying cells were generated by replacing a random subset (say, 10%) of the n × d cells by a value γ, which was varied to see its effect.

Outlying rows were generated in the hardest direction, that of the last eigenvector $\mathbf{v}$ of the true covariance matrix $\Sigma$. Next we rescale $\mathbf{v}$ to the typical size of a data point, by making ${\mathbf{v}}^{\text{'}}{\Sigma }^{-1}\mathbf{v}=E\left[Y^2\right]=d$ where  Y² ∼ χ²_d . We then replaced a random set of rows of $\mathbf{X}$ by $\gamma \mathbf{v}$.

Figure 7 shows the outlier misclassification rate (for n = 200 and d = 20). This is the number of cells that were generated as outlying and detected by the method, divided by the total number of outlying cells generated. This fraction depends on γ: it is high for small γ (but those cells or rows are not really outlying) and then goes down. In the top panels, we see that all methods for detecting cells work about equally well when the correlations between variables are fairly small, whereas the new method does better when there are at least some higher correlations. The bottom panels show a similar effect for detecting rows, and indicate that using the row filter, that is, Step 7 in Section 4, is an important component of DDC.

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 7. Outlier misclassification rates of different detection methods: (top) under cellwise contamination, (bottom) under rowwise contamination. Here GY is the Gervini–Yohai filter.

Display full size

The same simulation yields Figure 8, which shows the mean squared error of the imputed cells relative to the original cell values (before any contamination took place), in which flagged rows were taken out. Again DDC does best when there are some higher correlations. This confirms the usefulness of computing predicted cells in Step 4 of the DDC algorithm.

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 8. Mean squared error of the imputed cells in the same simulation.

Display full size

We have repeated the simulations leading to Figures 7 and 8 for dimensions d ranging up to 200, n between 10 and 2000, and contamination fractions up to 20%, each time with 50 replications. The method does not suffer from high dimensions, and even for n = 20 cases in d = 200 dimensions the comparisons between these methods yield qualitatively the same conclusions.

6. Two-Step Methods for Location and Scatter

Now assume that we are interested in estimating the parameters of the multivariate Gaussian distribution $N(\mu ,\Sigma )$ from which the data were generated (before cells and/or rows were contaminated).

If the data were clean, one would estimate $\mu$ and $\Sigma$ by the mean of the data and the classical covariance matrix. However, when the data may contain both cellwise and rowwise outliers, the problem becomes much more difficult. For this Agostinelli et al. (2015) proposed a two-step procedure called 2SGS, which is the current best method. The first step applies the univariate GY filter (from the previous section) to each column of the data matrix $\mathbf{X}$ and sets the flagged cells to NA. The second step then applies the sophisticated generalized S-estimator (GSE) of Danilov, Yohai, and Zamar (2012) to this incomplete dataset, yielding $\widehat{\mu }$ and $\widehat{\Sigma }$. The GSE is a rowwise robust estimator of $\mu$ and $\Sigma$ that was designed to work with data containing missing values, following earlier work by Little (1988) and Cheng and Victoria-Feser (2002).

Our version of this is to replace GY in the first step by DDC, followed by the same second step. When the first step flags a row, we take it out of the subsequent computations. We also include the Huberized Stahel-Donoho (HSD) estimator of Van Aelst, Vandervieren, and Willems (2012), as well as the minimum covariance determinant (MCD) estimator (Rousseeuw 1985; Rousseeuw and Van Driessen 1999), which is robust to rowwise outliers but not to cellwise outliers. For each method, we measure how far $\widehat{\Sigma }$ is from the true $\Sigma$ by the likelihood ratio type deviation $$\begin{array}{c}\hfill \text{LRT}=\text{trace}\left(\widehat{\Sigma }{\Sigma }^{-1}\right)-\text{log}\left(\text{det}\left(\widehat{\Sigma }{\Sigma }^{-1}\right)\right)-d\end{array}$$(which is a special case of the Kullback–Leibler divergence) and average this quantity over all replications.

Figure 9 compares these methods in the same simulation settings as Figures 7 and 8. In the top panels, we see that the new method performs about as well as 2SGS when the correlations between the variables are fairly low, but does much better when there are some higher correlations. For rowwise outliers their performance is quite similar.

Detecting Deviating Data Cells

All authors

Peter J. Rousseeuw & Wannes Van Den Bossche

https://doi.org/10.1080/00401706.2017.1340909

Published online:

06 December 2017

Figure 9. LRT deviation of three estimates from the true scatter matrix.

Display full size

Agostinelli et al. (2015) showed that 2SGS is consistent, that is, it gets the right answer for data generated from the model without any contamination when n goes to infinity. The proof follows from the fact that the fraction of cells flagged by the univariate GY filter goes to zero in that setting. This is not true for DDC because the cutoff value (5) is fixed, so some cells will still be flagged. Nevertheless the above simulations indicate that with actual contamination, using DDC in the first step often does better than GY.

A limitation of the GSE in the second step is that it requires n > d and in fact the dimension d should not be above 20 or 30, whereas the raw DDC method can deal with higher dimensions as we saw in the glass example with d = 750.

7. Conclusions and Outlook

The proposed method detects outliers at the level of the individual cells (entries) of the data matrix, unlike existing outlier detection methods, which consider the rows (cases) of that matrix as the basic units. Its main construct is the prediction of each cell. This turns a high dimensionality into a blessing instead of a curse, as having more variables (columns) available increases the amount of information and may improve the accuracy of the predicted cells.

The new method requires more computation than considering each variable in isolation, but on the other hand is able to detect outliers that would otherwise remain hidden as we saw in the first example. In simulations, we saw that DetectDeviatingCells performs as well as columnwise outlier detection when there is little correlation between the columns, whereas it does much better when there are higher correlations. Also, using it as the first step in the method of Agostinelli et al. (2015) outperforms a columnwise start.

A topic for further research is to extend this work to nonnumeric variables, such as binary and nominal. For the interaction between numeric and nonnumeric variables, and between nonnumeric variables, this necessitates replacing correlation by other measures of bivariate association. Also, for predicting cells the linear regression would need to be replaced by other techniques such as logistic regression.

8. Software Availability

The R code of the DetectDeviatingCells algorithm is available on CRAN in the package cellWise, which also contains functions for drawing cell maps and allows to reproduce all the examples in this article. Equivalent MATLAB code is available from the website http://wis.kuleuven.be/stat/robust/software.