1 Introduction

Real data often contain outliers, which can create serious problems when analyzing it. Many methods have been developed to deal with outliers, often by constructing a fit that is robust to them and then detecting the outliers by their large deviation (distance, residual) from that fit. For a brief overview of this approach see Rousseeuw and Hubert (2018). Unfortunately, most robust methods cannot handle data with missing values, some rare exceptions being Cheng and Victoria-Feser (2002) and Danilov, Yohai, and Zamar (2012). Moreover, they are typically restricted to casewise outliers, which are cases that deviate from the majority. We call these rowwise outliers because multivariate data are typically represented by a rectangular matrix in which the rows are the cases and the columns are the variables (measurements). In general, robust methods require that fewer than half of the rows are outlying, see, for example, Lopuhaä and Rousseeuw (1991). However, recently a different type of outliers, called cellwise outliers, have received much attention (Alqallaf et al. 2009; Van Aelst, Vandervieren, and Willems 2012; Agostinelli et al. 2015). These are suspicious cells (entries) that can occur anywhere in the data matrix. Figure 1 illustrates the difference between these types of outliers. The regular cells are shown in gray, whereas black means outlying. Rows 3 and 7 are rowwise outliers, and the other rows contain a fairly small percentage of cellwise outliers. As in this example, a small proportion of outlying cells can contaminate over half the rows, which causes most methods to break down. This effect is at its worst when the dimension (the number of columns) is high.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 1 Data matrix with missing values and cellwise and rowwise contamination.

Display full size

In high-dimensional situations, which are becoming increasingly common, one often applies principal component analysis (PCA) to reduce the dimension. However, the classical PCA (CPCA) method is not robust to either rowwise or cellwise outliers. Robust PCA methods that can deal with rowwise outliers include Croux and Ruiz-Gazen (2005), Hubert, Rousseeuw, and Verboven (2002), Locantore et al. (1999), Maronna (2005), and the ROBPCA method (Hubert, Rousseeuw, and Vanden Branden 2005). The latter method combines projection pursuit ideas with robust covariance estimation.

To deal with missing values, Nelson, Taylor, and MacGregor (1996) and Kiers (1997) developed the iterative classical PCA algorithm (ICPCA), see Walczak and Massart (2001) for a tutorial. The ICPCA follows the spirit of the EM algorithm. It starts by replacing the missing values by initial estimates such as the columnwise means. Then it iteratively fits a CPCA, yielding scores that are transformed back to the original space resulting in new estimates for the missing values, until convergence.

Serneels and Verdonck (2008) proposed a rowwise robust PCA method that can also cope with missing values. We will call this method MROBPCA (ROBPCA for missing values) as its key idea is to combine the ICPCA and ROBPCA methods. MROBPCA starts by imputing the NAs by robust initial estimates. The main difference with the ICPCA algorithm is that in each iteration the PCA model is fit by ROBPCA, which yields different imputations and flags rowwise outliers.

As of yet there are no PCA methods that can deal with cellwise outliers in combination with rowwise outliers and NAs. This article aims to fill that gap by constructing a new method called MacroPCA, where “Macro” stands for Missingness And Cellwise and Rowwise Outliers. It starts by applying a multivariate method called DetectDeviatingCells (Rousseeuw and Van den Bossche 2018) for detecting cellwise outliers, which provides initial imputations for the outlying cells and the NAs as well as an initial measure of rowwise outlyingness. In the next steps, MacroPCA combines ICPCA and ROBPCA to protect against rowwise outliers and to create improved imputations of the outlying cells and missing values. MacroPCA also provides graphical displays to visualize the different types of outliers. R code for MacroPCA is publicly available (Section 8).

2 The MacroPCA Algorithm

2.2 Dealing With Missing Values and Cellwise and Rowwise Outliers

We propose the MacroPCA algorithm for analyzing data that may contain one or more of the following issues: missing values, cellwise outliers, and rowwise outliers. Throughout the algorithm we will use the following two notations:

• the NA-imputed matrix ${\stackrel{°}{\mathit{X}}}_{n,d}$ only imputes the missing values of ${\mathit{X}}_{n,d}$;

• the cell-imputed matrix ${\stackrel{•}{\mathit{X}}}_{n,d}$ has imputed values for the outlying cells that do not belong to outlying rows, and for all missing values.

Both of these matrices still have n rows. Neither is intended to simply replace the true data matrix ${\mathit{X}}_{n,d}$. Note that ${\stackrel{•}{\mathit{X}}}_{n,d}$ does not try to impute outlying cells inside outlying rows, which would mask these rows in subsequent computations.

Since we do not know in advance which cells and rows are outlying, the set of flagged cellwise and rowwise outliers (and hence ${\stackrel{°}{\mathit{X}}}_{n,d}$ and ${\stackrel{•}{\mathit{X}}}_{n,d}$) will be updated in the course of the algorithm.

The first part of MacroPCA is the DetectDeviatingCells (DDC) algorithm. The description of this method can be found in Rousseeuw and Van den Bossche (2018) and in Section 1 of the supplementary material. The main purpose of the DDC method is to detect cellwise outliers. DDC outputs their positions $I_{c,\text{DDC}}$ as well as imputations for these outlying cells and any missing values. It also yields an initial outlyingness measure on the rows, which is, however, not guaranteed to flag all outlying rows. The set of flagged rows $I_{r,\text{DDC}}$ will be improved in later steps.

The second part of MacroPCA constructs principal components along the lines of the ICPCA algorithm but employing a version of ROBPCA (Hubert, Rousseeuw, and Vanden Branden 2005) to fit subspaces. It consists of the following steps, with all notations listed in Section 2 of the supplementary material.

1. Projection pursuit. The goal of this step is to provide an initial indication of which rows are the least outlying. For this ROBPCA starts by identifying the h < n least outlying rows by a projection pursuit procedure. We write $0.5 ⩽ \alpha =h/n<1$. This means that we can withstand up to a fraction $1-\alpha$ of outlying rows. To be on the safe side the default is $\alpha =0.5$.

   However, due to cellwise outliers there may be far fewer than h uncontaminated rows, so we cannot apply this step to the original data ${\mathit{X}}_{n,d}$. We also cannot use the entire imputed matrix ${\tilde{\mathit{X}}}_{n,d}$ obtained from DDC in which all outlying cells are imputed, even those in potentially outlying rows, as this could mask outlying rows. Instead we use the cell-imputed matrix ${\stackrel{•}{\mathit{X}}}_{n,d}^{(0)}$ defined as follows:

   1. In all rows we replace the missing values by the values $\stackrel{°}{{\mathit{x}}_i^{(0)}}$ as imputed by DDC.

   2. Next, in the h rows with the fewest cells flagged by DDC but not in $I_{r,\text{DDC}}$ we impute those flagged cells, that is, $\stackrel{•}{{\mathit{x}}_i^{(0)}}={\tilde{\mathit{x}}}_i$.

      As in ROBPCA the outlyingness of a point $\stackrel{•}{{\mathit{x}}_i^{(0)}}$ is then computed as(2) $$\text{outl}\left(\stackrel{•}{{\mathit{x}}_i^{\left(0\right)}}\right)=\underset{\mathit{v}\in B}{\max }\frac{|\mathit{v}\prime \stackrel{•}{{\mathit{x}}_i^{\left(0\right)}}-m_{\text{MCD}}(\mathit{v}\prime \stackrel{•}{{\mathit{x}}_j^{\left(0\right)}}\left)\right|}{s_{\text{MCD}}(\mathit{v}\prime \stackrel{•}{{\mathit{x}}_j^{\left(0\right)}})},$$(2)

      where $m_{\text{MCD}}(\mathit{v}\prime \stackrel{•}{{\mathit{x}}_j^{(0)}})$ and $s_{\text{MCD}}(\mathit{v}\prime \stackrel{•}{{\mathit{x}}_j^{(0)}})$ are the univariate MCD location and scale estimates (Rousseeuw and Leroy 1987) of $\{\mathit{v}\prime \stackrel{•}{{\mathit{x}}_1^{(0)}},\dots ,\mathit{v}\prime \stackrel{•}{{\mathit{x}}_n^{(0)}}\}$. The set B contains 250 directions through two data points (or all of them if there are fewer than 250). Finally, the indices of the h rows $\stackrel{•}{{\mathit{x}}_i^{(0)}}$ with the lowest outlyingness and not belonging to $I_{r,\text{DDC}}$ are stored in the set H₀.

2. Subspace dimension. Here, we choose the number of principal components. For this we build a new cell-imputed matrix ${\stackrel{•}{\mathit{X}}}_{n,d}^{(1)}$ which imputes the outlying cells in the rows of H₀ and imputes the NAs in all rows. This means that $\stackrel{•}{{\mathit{x}}_i^{(1)}}={\tilde{\mathit{x}}}_i$ for $i\in H_0$, and $\stackrel{•}{{\mathit{x}}_i^{(1)}}=\stackrel{°}{{\mathit{x}}_i^{(0)}}$ if $i\notin H_0$. Then we apply classical PCA to the $\stackrel{•}{{\mathit{x}}_i^{(1)}}$ with $i\in H_0$. Their mean ${\mathit{m}}_d^{(1)}$ is an estimate of the center, whereas the spectral decomposition of their covariance matrix yields a loading matrix ${\mathit{P}}_{d,d}^{(1)}$ and a diagonal matrix ${\mathit{L}}_{d,d}^{(1)}$ with the eigenvalues sorted from largest to smallest. These eigenvalues can be used to construct a screeplot from which an appropriate dimension k of the subspace can be derived. Alternatively, one can retain a certain cumulative proportion of explained variance, such as 80%. The maximal number of principal components that MacroPCA will consider is the tuning constant $k_{\text{max}}$, which is set to 10 by default.

3. Iterative subspace estimation. This step aims to estimate the k-dimensional subspace fitting the data. As in ICPCA this requires iteration, for $s ⩾ 2$:

   1. The scores matrix in Equation (1) based on the cell-imputed cases is computed as ${\stackrel{•}{\mathit{T}}}_{n,k}^{(s-1)}=({\stackrel{•}{\mathit{X}}}_{n,d}^{(s-1)}-1_n({\mathit{m}}_d^{(s-1)})\prime ){\mathit{P}}_{d,k}^{(s-1)}$. The predicted data values are set to ${\widehat{\mathit{X}}}_{n,d}^{(s)}=1_n({\mathit{m}}_d^{(s-1)})\prime +{\stackrel{•}{\mathit{T}}}_{n,k}^{(s-1)}({\mathit{P}}_{d,k}^{(s-1)})\prime$. We then update the imputed matrices to ${\stackrel{°}{\mathit{X}}}_{n,d}^{(s)}$ and ${\stackrel{•}{\mathit{X}}}_{n,d}^{(s)}$ by replacing the appropriate cells by the corresponding cells of ${\widehat{\mathit{X}}}_{n,d}^{(s)}$. That is, for ${\stackrel{°}{\mathit{X}}}_{n,d}^{(s)}$ we update all the imputations of missing cells, whereas for ${\stackrel{•}{\mathit{X}}}_{n,d}^{(s)}$ we update the imputations of the outlying cells in rows of H₀ as well as the NAs in all rows.

   2. The PCA model is reestimated by applying classical PCA to the $\stackrel{•}{{\mathit{x}}_i^{(s)}}$ with $i\in H_0$. This yields a new estimate ${\mathit{m}}_d^{(s)}$ as well as an updated loading matrix ${\mathit{P}}_{d,k}^{(s)}$.

      The iterations are repeated until s = 20 or until convergence is reached, that is, when the maximal angle between a vector in the new subspace and the vector most parallel to it in the previous subspace is below some tolerance (by default 0.005). Following Krzanowski (1979) this angle is computed as $\mathit{arccos}(\sqrt{{\delta }_k})$ where δ_k is the smallest eigenvalue of $({\mathit{P}}_{d,k}^{(s)})\prime {\mathit{P}}_{d,k}^{(s-1)}({\mathit{P}}_{d,k}^{(s-1)})\prime {\mathit{P}}_{d,k}^{(s)}$.

      After all iterations we have the NA-imputed matrix ${\stackrel{°}{\mathit{X}}}_{n,d}^{(s)}$, the cell-imputed matrix ${\stackrel{•}{\mathit{X}}}_{n,d}^{(s)}$ as well as the estimated center ${\mathit{m}}_d^{(s)}$ and the loading matrix ${\mathit{P}}_{d,k}^{(s)}$.

4. Reweighting. In robust statistics one often follows an initial estimate by a reweighting step to improve the statistical efficiency at a low computational cost, see, for example, Rousseeuw and Leroy (1987) and Engelen, Hubert, and Vanden Branden (2005). Here, we use the orthogonal distance of each $\stackrel{•}{{\mathit{x}}_i^{(s)}}$ to the current PCA subspace:$$\stackrel{•}{{\text{OD}}_i}=\left|\right|\stackrel{•}{{\mathit{x}}_i^{\left(s\right)}}-{\widehat{\stackrel{•}{\mathit{x}}}}_i^{\left(s\right)}\left|\right|=\left|\right|\stackrel{•}{{\mathit{x}}_i^{\left(s\right)}}-({\mathit{m}}_d^{\left(s\right)}+(\stackrel{•}{{\mathit{x}}_i^{\left(s\right)}}-{\mathit{m}}_d^{\left(s\right)}\left){\mathit{P}}_{d,k}^{\left(s\right)}\right({\mathit{P}}_{d,k}^{\left(s\right)})\prime )\left|\right|.$$

   The orthogonal distances to the power 2/3 are roughly Gaussian except for the outliers (Hubert, Rousseeuw, and Vanden Branden 2005), so we compute the cutoff value(3) $$c_{\text{OD}}:={\left(m_{\text{MCD}}\left(\right\{\stackrel{•}{{\text{OD}}_j^{2/3}}\left\}\right)+s_{\text{MCD}}\left(\right\{\stackrel{•}{{\text{OD}}_j^{2/3}}\left\}\right){\Phi }^{-1}\left(0.99\right)\right)}^{3/2}.$$(3) All cases for which $\stackrel{•}{{\text{OD}}_i} ⩽ c_{\text{OD}}$ are considered non-outlying with respect to the PCA subspace, and their indices are stored in a set $H^{*}$. As before, any $i\in I_{r,\mathit{DDC}}$ is removed from $H^{*}$. The final cell-imputed matrix ${\stackrel{•}{\mathit{X}}}_{n,d}$ is given by ${\stackrel{•}{x}}_{i,j}={\widehat{\stackrel{•}{x}}}_{i,j}^{(s)}$ if $i\in H^{*}$ and $j\in I_{c,\mathit{DDC}}$ and ${\stackrel{•}{x}}_{i,j}={\stackrel{°}{x}}_{i,j}^{(s)}$ otherwise. Applying classical PCA to the $n^{*}$ rows $\stackrel{•}{{\mathit{x}}_i}$ in $H^{*}$ yields a new center ${\mathit{m}}_d^{*}$ and a new loading matrix ${\mathit{P}}_{d,k}^{*}$.

5. DetMCD. Now we want to estimate a robust basis of the estimated subspace. The columns of ${\mathit{P}}_{d,k}^{*}$ from Step 4 need not be robust, because some of the $n^{*}$ rows in $H^{*}$ might be outlying inside the subspace. These so-called good leverage points do not harm the estimation of the PCA subspace but they can still affect the estimated eigenvectors and eigenvalues, as illustrated by a toy example in Section A.1 of the Appendix. In this step, we first project the $n^{*}$ points of $H^{*}$ onto the subspace, yielding$${\stackrel{•}{\mathit{T}}}_{n^{*},k}=\left({\stackrel{•}{\mathit{X}}}_{n^{*},d}-1_{n^{*}}{\mathit{m}}_d^{*\prime }\right){\mathit{P}}_{d,k}^{*}.$$

   Next, the center and scatter matrix of the scores ${\stackrel{•}{\mathit{T}}}_{n^{*},k}$ are estimated by the DetMCD method of Hubert, Rousseeuw, and Verdonck (2012). This is a fast, robust and deterministic algorithm for multivariate location and scatter, yielding ${\mathit{m}}_k^{\text{MCD}}$ and ${\mathit{S}}_{k,k}^{\text{MCD}}$. Its computation is feasible because the dimension k of the subspace is quite low. The spectral decomposition of ${\mathit{S}}_{k,k}^{\text{MCD}}$ yields a loading matrix ${\mathit{P}}_{k,k}^{\text{MCD}}$ and eigenvalues ${\widehat{\lambda }}_j$ for $j=1,\dots ,k$. We set the final center to ${\mathit{m}}_d={\mathit{m}}_d^{*}+{\mathit{P}}_{d,k}^{*}{\mathit{m}}_k^{\text{MCD}}$ and the final loadings to ${\mathit{P}}_{d,k}={\mathit{P}}_{d,k}^{*}{\mathit{P}}_{k,k}^{\text{MCD}}$.

6. Scores, predicted values, and residuals. We now provide the final output. We compute the scores of ${\stackrel{°}{\mathit{X}}}_{n,d}$ as ${\stackrel{°}{\mathit{T}}}_{n,k}=({\stackrel{°}{\mathit{X}}}_{n,d}-1_n{\mathit{m}}_d^{\prime }){\mathit{P}}_{d,k}$ and the predictions of ${\stackrel{°}{\mathit{X}}}_{n,d}$ as ${\widehat{\stackrel{°}{\mathit{X}}}}_{n,d}=1_n{\mathit{m}}_d^{\prime }+{\stackrel{°}{\mathit{T}}}_{n,k}({\mathit{P}}_{d,k})\prime$. (The formulas for ${\stackrel{•}{\mathit{T}}}_{n,k}$ and ${\widehat{\stackrel{•}{\mathit{X}}}}_{n,d}$ are analogous.) This yields the difference matrix ${\stackrel{°}{\mathit{X}}}_{n,d}-{\widehat{\stackrel{°}{\mathit{X}}}}_{n,d}$, which we then robustly scale by column, yielding the final standardized residual matrix ${\mathit{R}}_{n,d}$. The orthogonal distance of $\stackrel{°}{{\mathit{x}}_i}$ to the PCA subspace is given by(4) $$\stackrel{°}{{\text{OD}}_i}=\left|\right|\stackrel{°}{{\mathit{x}}_i}-{\widehat{\stackrel{°}{\mathit{x}}}}_i\left|\right|.$$(4)

See Section 8 for the R code carrying out MacroPCA.

MacroPCA can be carried out in $O(\mathit{nd}(\min (n,d)+\hspace{0.17em}\text{log}\hspace{0.17em}(n)+\hspace{0.17em}\text{log}\hspace{0.17em}(d)))$ time (see Section A.2 of the Appendix) which is not much more than the complexity $O(\mathit{nd}\min (n,d))$ of classical PCA. Figure 2 shows times as a function of n and d indicating that MacroPCA is quite fast. The fraction of NAs in the data had no substantial effect on the computation time, as seen in Figure A.3 in Section A.2, Appendix.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 2 Computation times of MacroPCA in seconds on Intel i7-4800MQ at 2.70 GHz, as a function of the number of cases n (left) and of the dimension d (right).

Display full size

Fig. 2 Computation times of MacroPCA in seconds on Intel i7-4800MQ at 2.70 GHz, as a function of the number of cases n (left) and of the dimension d (right).

Note that PCA loadings are highly influenced by the variables with the largest variability. For this the MacroPCA code provides the option to divide each variable by a robust scale. This does not increase the computational complexity.

3 Outlier Detection

MacroPCA provides several tools for outlier detection. We illustrate them on a dataset collected by Alfons (2016) from the website of the Top Gear car magazine. It contains data on 297 cars, with 11 continuous variables. Five of these variables (price, displacement, BHP, torque, top speed) are highly skewed, and were logarithmically transformed. The dataset contains 95 missing cells, which is only 2.9% of the 297 × 11 = 3267 cells. We retained two principal components (k = 2).

The right-hand panel of Figure 3 shows the results of MacroPCA by a modification of the cell map introduced by Rousseeuw and Van den Bossche (2018). The computations were performed on all 297 cars, but to make the map fit on a page it only shows 24 cars, including some of the more eventful cases. The color of the cells stems from the standardized residual matrix ${\mathit{R}}_{n,d}$ obtained by MacroPCA. Cells with $\left|r_{\mathit{ij}}\right| ⩽ \sqrt{{\chi }_{1,0.99}^2}=2.57$ are considered regular and colored yellow in the residual map, whereas the missing values are white. Outlying residuals receive a color which ranges from light orange to red when r_ij > 2.57 and from light purple to dark blue when $r_{\mathit{ij}}<-2.57$. So a dark red cell indicates that its observed value is much higher than its fitted value, while a dark blue cell means the opposite.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 3 Residual map of selected rows from Top Gear data: (left) when using ICPCA; (right) when using MacroPCA. The numbers shown in the cells are the original data values (with price in units of 1000 UK Pounds).

Display full size

To the right of each row in the map is a circle whose color varies from white to black according to the orthogonal distance $\stackrel{°}{{\text{OD}}_i}$ given by Equation (4) compared to the cutoff Equation (3). Cases with $\stackrel{°}{{\text{OD}}_i} ⩽ c_{\text{OD}}$ lie close to the PCA subspace and receive a white circle. The others are given darker shades of gray up to black according to their $\stackrel{°}{{\text{OD}}_i}$.

On these data, we also ran the ICPCA method, which handles missing values in classical PCA. It differs from MacroPCA in some important ways: the initial imputations are by nonrobust column means, the iterations carry out CPCA and do not exclude outlying rows, and the residuals are standardized by the nonrobust standard deviation. By itself ICPCA does not provide a residual map, but we can construct one anyway by plotting the nonrobust standardized residuals with the same color scheme, yielding the left panel of Figure 3.

The ICPCA algorithm finds high orthogonal distances (dark circles) for the BMW i3, the Chevrolet Volt, the Renault Twizzy, and the Vauxhall Ampera. These are hybrid or purely electrical cars with a high or missing MPG (miles per gallon). Note that the Ssangyong Rodius and Renault Twizzy get blue cells for their acceleration time of zero seconds, which is physically impossible.

On this dataset the ICPCA algorithm provides decent results because the total number of outliers is small compared to the size of the data, and indeed the residual map of all 297 cars was mostly yellow. But MacroPCA (right panel) detects more deviating behavior. The orthogonal distance of the hybrid Citroen DS5 and the electrical Mitsubishi i-MiEV are now on the high side, and the method flags the Bugatti Veyron and Pagani Huayra supercars as well as the Land Rover Defender and Mercedes-Benz G all-terrain vehicles. It also flags more cells, giving a more complete picture of the special characteristics of some cars.

We can also compute the score distance of each case, which is the robustified Mahalanobis distance of its projection on the PCA subspace among all such projected points. It is easily computed as(5) $$\stackrel{°}{{\text{SD}}_i}=\sqrt{\sum _{j=1}^k{\left(\stackrel{°}{t_{\mathit{ij}}}\right)}^2/{\widehat{\lambda }}_j},$$(5) where $\stackrel{°}{t_{\mathit{ij}}}$ are the scores and ${\widehat{\lambda }}_j$ the eigenvalues obtained by MacroPCA. This allows us to construct a PCA outlier map of cases as introduced in Hubert, Rousseeuw, and Vanden Branden (2005), which plots the orthogonal distances $\stackrel{°}{{\text{OD}}_i}$ on the vertical axis versus the score distances $\stackrel{°}{{\text{SD}}_i}$. The MacroPCA outlier map of these data is the right panel of Figure 4. The vertical line indicates the cutoff $c_{\text{SD}}=\sqrt{{\chi }_{k,0.99}^2}$ and the horizontal line is the cutoff $c_{\text{OD}}$.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 4 Outlier map of Top Gear data: (left) when using ICPCA; (right) when using MacroPCA.

Display full size

Regular cases are those with a small $\stackrel{°}{{\text{SD}}_i} ⩽ c_{\text{SD}}$ and a small $\stackrel{°}{{\text{OD}}_i} ⩽ c_{\text{OD}}$. Cases with large $\stackrel{°}{{\text{SD}}_i}$ and small $\stackrel{°}{{\text{OD}}_i}$ are called good leverage points. The cases with large $\stackrel{°}{{\text{OD}}_i}$ can be divided into orthogonal outliers (when their $\stackrel{°}{{\text{SD}}_i}$ is small) and bad leverage points (when their $\stackrel{°}{{\text{SD}}_i}$ is large too). We see several orthogonal outliers such as the Vauxhall Ampera as well as some bad leverage points, especially the BMW i3. There are also some good leverage points.

The left panel displays the outlier map for ICPCA. It flags the BMW i3 as an orthogonal outlier. This behavior is typical because a bad leverage point will attract the fit of classical methods, making it appear less special. For the same reason ICPCA considers some of the good leverage points as regular cases. That the ICPCA outlier map is still able to flag some outliers is due to the fact that this dataset only has a small percentage of outlying rows.

4 Online Data Analysis

Applying MacroPCA to a dataset ${\mathit{X}}_{n,d}$ yields a PCA fit. Now suppose that a new case (row) x comes in, and we would like to impute its missing values, detect its outlying cells and impute them, estimate its scores, and find out whether or not it has an outlying score and/or a high orthogonal distance. We could of course append x to ${\mathit{X}}_{n,d}$ and rerun MacroPCA, but that would be very inefficient.

Instead we propose a method to analyze x using only the output of MacroPCA on the initial set ${\mathit{X}}_{n,d}$. This can be done quite fast, which makes the procedure suitable for online process control. For outlier-free data with NAs this was studied by Nelson, Taylor, and MacGregor (1996) and Walczak and Massart (2001). Folch-Fortuny, Arteaga, and Ferrer (2015) call this model exploitation, as opposed to model building (fitting a PCA model). Our procedure consists of two stages, along the lines of MacroPCA.

1. DDCpredict is a new function which only uses x and the output of DDC on the initial data ${\mathit{X}}_{n,d}$. First the entries of x are standardized using the robust location and scale estimates from DDC. Then all x_j with $\left|x_j\right|>\sqrt{{\chi }_{1,0.99}^2}=2.57$ are replaced by NAs. Next, all NAs are estimated as in DDC making use of the pre-built coefficients b_jh and weights w_jh. Also the deshrinkage step uses the original robust slopes. The DDCpredict stage yields the imputed vector ${\tilde{\mathit{x}}}^{(0)}$ and the standardized residual of each cell x_j.

2. MacroPCApredict improves on the initial imputation ${\tilde{\mathit{x}}}^{(0)}$. The improvements are based solely on the ${\mathit{m}}_d$ and ${\mathit{P}}_{d,k}$ that were obtained by MacroPCA on the original data ${\mathit{X}}_{n,d}$. Step $s ⩾ 1$ is of the following form:

   1. Project the imputed case ${\tilde{\mathit{x}}}^{(s-1)}$ on the MacroPCA subspace to obtain its scores vector ${\mathit{t}}^{(s)}=({\mathit{P}}_{d,k})\prime ({\tilde{\mathit{x}}}^{(s-1)}-{\mathit{m}}_d)$;

   2. transform the scores to the original space, yielding ${\widehat{\mathit{x}}}^{(s)}={\mathit{m}}_d+{\mathit{P}}_{d,k}{\mathit{t}}^{(s)}$;

   3. reimpute the outlying cells and missing values of x by the corresponding values of ${\widehat{\mathit{x}}}^{(s)}$, yielding ${\tilde{\mathit{x}}}^{(s)}$.

These steps are iterated until convergence (when the new imputed values are within a tolerance of the old ones) or the maximal number of steps (by default 20) is reached. We denote the final ${\tilde{\mathit{x}}}^{(s)}$ as $\tilde{\mathit{x}}$.

Next, we create $\stackrel{°}{\mathit{x}}$ by replacing the missing values in x by the corresponding cells in $\tilde{\mathit{x}}$. We then compute the orthogonal distance $\text{OD}(\stackrel{°}{\mathit{x}})$ and the score distance $\text{SD}(\stackrel{°}{\mathit{x}})$, and classify it as a regular case, a leverage point or an orthogonal outlier. Finally the cell residuals $\stackrel{°}{x_j}-{\widehat{\stackrel{°}{x}}}_j$ are standardized as in the last step of MacroPCA, and used to flag outlying cells in x.

To illustrate this prediction procedure we reanalyze the Top Gear dataset. We exclude the 24 cars shown in the residual map of Figure 3 and build the MacroPCA model on the remaining data. This model was then provided to analyze the 24 selected cars as “new” data. Figure 5 shows the result. As before the cells are colored according to their standardized residual, and the circles on the right are filled according to their $\stackrel{°}{\text{OD}}$. The left panel is the MacroPCA residual map shown in Figure 3, which was obtained by applying MacroPCA to the entire dataset. The right panel shows the result of analyzing these 24 cases using the fit obtained without them. The residual maps are quite similar. Note that each cell now shows its standardized residual (instead of its data value as in Figure 3), making it easier to see the differences.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 5 Top Gear dataset: residual maps obtained by (left) including and (right) excluding these 24 cars when fitting the PCA model.

Display full size

5 Simulations

We have compared the performance of ICPCA, MROBPCA, and MacroPCA in an extensive simulation study. Several contamination models were considered with missing values, cellwise outliers, rowwise outliers, and combinations of them. Only a few of the results are reported here since the others yielded similar conclusions.

The clean data ${\mathit{X}}_{n,d}^0$ are generated from a multivariate Gaussian with $\mathit{\mu }=0$ and two types of covariance matrices ${\mathbf{\Sigma }}_{d,d}$. The first one is based on the structured correlation matrix called A09 where each off-diagonal entry is ${\rho }_{i,j}={\left(-0.9\right)}^{|i-j|}$. The second type of covariance matrix is based on the random correlation matrices of Agostinelli et al. (2015) and will be called ALYZ. These correlation matrices are turned into covariance matrices with other eigenvalues. More specifically, the diagonal elements of the matrix ${\mathit{L}}_{d,d}$ from the spectral decomposition ${\mathbf{\Sigma }}_{d,d}={\mathit{P}}_{d,d}{\mathit{L}}_{d,d}\mathit{P}{\prime }_{d,d}$ are replaced by the desired values listed below. The specifications of the clean data are n = 100, d = 200, ${\mathit{L}}_{d,d}=\text{diag}(30,25,20,15,10,5,0.098,0.0975,\dots ,0.0020,0.0015)$, and k = 6 (since ${\sum }_{j=1}^6{\lambda }_j/{\sum }_{j=1}^{200}{\lambda }_j=91.5\%$). MacroPCA takes less than a second for n = 100, d = 200 as seen in Figure 2.

In a first simulation setting, the clean data ${\mathit{X}}_{n,d}^0$ are modified by replacing a random subset of 5%, 10%,… up to 30% of the n × d cells with NAs. The second simulation setting generates NAs and outlying cells by randomly replacing 20% of the cells x_ij by missing values and 20% by the value $\gamma {\sigma }_j$ where ${\sigma }_j^2$ is the jth diagonal element of ${\mathbf{\Sigma }}_{d,d}$ and γ ranges from 0 to 20. The third simulation setting generates NAs and outlying rows. Here 20% of random cells are replaced by NAs and a random subset of 20% of the rows is replaced by rows generated from $N(\gamma {\mathit{v}}_{k+1},{\mathbf{\Sigma }}_{d,d})$, where γ varies from 0 to 50 and ${\mathit{v}}_{k+1}$ is the $(k+1)$th eigenvector of ${\mathbf{\Sigma }}_{d,d}$. The last simulation setting generates 20% of NAs, together with 10% of cellwise outliers and 10% of rowwise outliers in the same way.

In each setting, we consider the set C consisting of the rows i that were not replaced by rowwise outliers, with $c:=\#C$, and the data matrix ${\mathit{X}}_{c,d}^0$ consisting of those rows of the clean data ${\mathit{X}}_{n,d}^0$. As a baseline for the simulation we apply classical PCA to ${\mathit{X}}_{c,d}^0$ and denote the resulting predictions by ${\widehat{x}}_{\mathit{ij}}^C$ for i in C. We then measure the mean squared error (MSE) from the baseline:$$\text{MSE}=\frac{1}{\mathit{cd}}\sum _{i\in C}\sum _{j=1}^d{\left({\widehat{x}}_{\mathit{ij}}-{\widehat{x}}_{\mathit{ij}}^C\right)}^2,$$where ${\widehat{x}}_{\mathit{ij}}$ is the predicted value for x_ij obtained by applying the different methods to the contaminated data. The MSE is then averaged over 100 replications.

Figure 6 shows the performance of ICPCA, MROBPCA, and MacroPCA when some data become missing. As CPCA and ROBPCA cannot deal with NAs, they are not included in this comparison. Since there are no outliers the classical ICPCA performs best, followed by MROBPCA and MacroPCA which perform similarly to each other, and only slightly worse than ICPCA considering the scale of the vertical axis which is much smaller than in the other three simulation settings.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 6 Average MSE as a function of the fraction ε of missing values. The data were generated using A09 (left) and ALYZ (right).

Display full size

Now we set 20% of the data cells to missing and add 20% of cellwise contamination given by γ. Figure 7 shows the performance of ICPCA, MROBPCA, and MacroPCA in this situation. The MSE of both ICPCA and MROBPCA grows very fast with γ which indicates that these methods are not at all robust to cellwise outliers. Note that d = 200 so on average $1-{(1-0.2)}^{200}\approx 100\%$ of the rows are contaminated, whereas no purely rowwise method can handle more than 50%. MacroPCA is the only method that can withstand cellwise outliers here. When γ is smaller than five the MSE goes up, but this is not surprising as in that case the values in the contaminated cells are still close to the clean ones. As soon as the contamination is sufficiently far away, the MSE drops to a very low value.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 7 Average MSE for data with 20% of missing values and 20% of cellwise outliers, as a function of γ which determines the distance of the cellwise outliers.

Display full size

Figure 8 presents the results of ICPCA, MROBPCA, and MacroPCA when there are 20% of missing values combined with 20% of rowwise contamination. As expected, the ICPCA algorithm breaks down while MROBPCA and MacroPCA provide very good results. MROBPCA and MacroPCA are affected the most (but not much) by nearby outliers, and very little by far contamination.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 8 Average MSE for data with 20% of missing values and 20% of rowwise outliers, as a function of γ, which determines the distance of the rowwise outliers.

Display full size

Finally, Figure 9 presents the results in the situation of 20% of missing values combined with 10% of cellwise and 10% of rowwise contamination. In this scenario the ICPCA and MROBPCA algorithms break down whereas MacroPCA still provides reasonable results.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 9 Average MSE for data with 20% of missing values, 10% of cellwise outliers and 10% of rowwise outliers, as a function of γ which determines the distance of both the cellwise and the rowwise outliers.

Display full size

In this section, the missing values were generated in a rather simple way. In Section A.3 they are generated in a more challenging way but still MAR, with qualitatively similar results.

6 Real Data Examples

6.1 Glass Data

The glass dataset (Lemberge et al. 2000) contains spectra with d = 750 wavelengths of n = 180 archeological glass samples. It is available in the R package cellWise (Raymaekers et al. 2019). The MacroPCA method selects four principal components and yields a 180 × 750 matrix of standardized residuals. There is not enough resolution on a page to show so many individual cells in a residual map. Therefore, we created a map (the top panel of Figure 10) which combines the residuals into blocks of 5 × 5 cells. The color of each block now depends on the most frequent type of outlying cell in it, the resulting color being an average. For example, an orange block indicates that quite a few cells in the block were red and most of the others were yellow. The more red cells in the block, the darker red the block will be. We see that MacroPCA has flagged a lot of cells, that happen to be concentrated in a minority of the rows where they show patterns. In fact, the colors indicate that some of the glass samples (between 22 and 30) have a higher concentration of phosphor, whereas rows 57–63 and 74–76 had an unusually high concentration of calcium. The bottom part of the residual map looks very different, due to the fact that the measuring instrument was cleaned before recording the last 38 spectra. One could say that those outlying rows belong to a different population.

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 10 Residual maps of the glass dataset when fitting the PCA model by MacroPCA (top) and ROBPCA (bottom).

Display full size

Fig. 10 Residual maps of the glass dataset when fitting the PCA model by MacroPCA (top) and ROBPCA (bottom).

Since the dataset has no NAs and we found that fewer than half of the rows are outlying, it can also be analyzed by the original ROBPCA method as was done by Hubert, Rousseeuw, and Vanden Branden (2005), also for k = 4. This detects the same rowwise outliers. In principle ROBPCA is a purely rowwise method that does not flag cells. Even though ROBPCA does not produce a residual map, we can construct one analogously to that of MacroPCA. First we construct the residual matrix of ROBPCA, the rows of which are given by ${\mathit{x}}_i-{\widehat{\mathit{x}}}_i$, where ${\widehat{\mathit{x}}}_i$ is the projection of x on the ROBPCA subspace. We then standardize the residuals in each column by dividing them by a robust one-step scale M-estimate. This yields the bottom panel of Figure 10. We see that the two residual maps look quite similar.

This example illustrates that purely rowwise robust methods can be useful to detect cellwise outliers when these cells occur in fewer than 50% of the rows. But if the cellwise outliers contaminate more rows, this approach is insufficient.

6.2 DPOSS Data

In our last example, we analyze data from the Digitized Palomar Sky Survey (DPOSS) described by Odewahn et al. (1998). This is a huge database of celestial objects, from which we have drawn 20,000 stars at random. Each star has been observed in the color bands J, F, and N. Each band has seven variables. Three of them measure light intensity: for the J band they are MAperJ, MTotJ, and MCoreJ, where the last letter indicates the band. The variable AreaJ is the size of the star based on its number of pixels. The remaining variables IR2J, csfJ, and EllipJ combine size and shape. (There were two more variables in the original data, but these measured the background rather than the star itself.) There are substantial correlations between these 21 variables.

In this dataset 84.6% of the rows contain NAs (in all there are 50.2% missing entries). Often an entire color band is missing, and sometimes two. We applied MacroPCA to these data, choosing k = 4 components according to the screeplot. The left panel of Figure 11 shows the loadings of the first and second component. It appears that the first component captures the overall negative correlation between two groups of variables: those measuring light intensity (the first three variables in each band) and the others (variables 4–7 in each band). The right panel is the corresponding scores plot, in which the 150 stars with the highest orthogonal distance $\stackrel{°}{\text{OD}}$ are shown in red. Most of these stand out in the space of PC1 and PC2 (bad leverage points), whereas some only have a high $\stackrel{°}{\text{OD}}$ (orthogonal outliers).

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 11 DPOSS stars data: (left) loadings of the first (black full line) and the second (blue dashed line) component of MacroPCA, with vertical lines separating the three color bands; (right) plot of the first two scores, with filled red circles for stars with high orthogonal distance $\stackrel{°}{\text{OD}}$ and open black circles for the others.

Display full size

Fig. 11 DPOSS stars data: (left) loadings of the first (black full line) and the second (blue dashed line) component of MacroPCA, with vertical lines separating the three color bands; (right) plot of the first two scores, with filled red circles for stars with high orthogonal distance $\stackrel{°}{\text{OD}}$ and open black circles for the others.

Figure 12 shows the residual map of MacroPCA, in which each row block combines 25 stars. The six rows at the top correspond to the 150 stars with highest $\stackrel{°}{\text{OD}}$. We note that the outliers tend to be more luminous (MTot) than expected and have a larger Area, which suggests giant stars. The analogous residual map of ICPCA (not shown) did not reveal much. Note that the non-outlying rows in the bottom part of the residual map are yellow, and the missing color bands show up as blocks in lighter yellow (a combination of yellow and white cells).

MacroPCA: An All-in-One PCA Method Allowing for Missing Values as Well as Cellwise and Rowwise Outliers

All authors

Mia Hubert , Peter J. Rousseeuw & Wannes Van den Bossche

https://doi.org/10.1080/00401706.2018.1562989

Published online:

13 May 2019

Fig. 12 MacroPCA residual map of stars in the DPOSS data, with 25 stars per row block. The six row blocks at the top correspond to the stars with highest $\stackrel{°}{\text{OD}}$.

Display full size

Fig. 12 MacroPCA residual map of stars in the DPOSS data, with 25 stars per row block. The six row blocks at the top correspond to the stars with highest $\stackrel{°}{\text{OD}}$.

7 Conclusions

The MacroPCA method is able to handle missing values, cellwise outliers, and rowwise outliers. This makes it well-suited for the analysis of possibly messy real data. Simulation showed that its performance is similar to a classical method in the case of outlier-free data with missing values, and to an existing robust method when the data only has rowwise outliers. The algorithm is fast enough to deal with many variables, and has recently been sped up by recoding parts in C$++$, making it much faster than the pure R code in Figure 2 and Figure A.3.

MacroPCA can analyze new data as they come in, only making use of its existing output obtained from the initial dataset. It imputes missing values in the new data, flags and imputes outlying cells, and flags outlying rows. This computation is fast, so it can be used to screen new data in quality control or even online process control. (One can update the initial fit offline from time to time.) The advantage of MacroPCA is that it not only tells us when the process goes out of control, but also which variables are responsible.

Potential extensions of MacroPCA include methods of PCA regression and partial least squares able to deal with rowwise and cellwise outliers and missing values.

8 Software Availability

The methods described in this article as well as the DPOSS data have been made available in the R package cellWise (Raymaekers et al. 2019) on CRAN, with a vignette reproducing the examples shown.

Supplementary Materials

This is a text with more details about the first part of the MacroPCA algorithm and a list of notations.