Sparse robust PCA algorithm based on the ROBPCA algorithm of Hubert et al. (2005).
rospca(X, k, kmax = 10, alpha = 0.75, h = NULL, ndir = "all", grid = TRUE,
lambda = 10^(-6), sparse = "varnum", para, stand = TRUE, skew = FALSE)A list with components:
Loadings matrix containing the sparse robust loadings (eigenvectors), a numeric matrix of size \(p\) by \(k\).
Numeric vector of length \(k\) containing the robust eigenvalues.
Scores matrix (computed as \((X-center) \cdot loadings)\), a numeric matrix of size \(n\) by \(k\).
Numeric vector of length \(k\) containing the centre of the data.
Matrix used to standardise the data before applying sparse PCA (identity matrix if stand=FALSE), a numeric matrix of size \(p\) by \(p\).
Number of (chosen) principal components.
Logical vector of size \(n\) indicating if an observation is in the initial h-subset.
Logical vector of size \(n\) indicating if an observation is kept in the non-sparse reweighting step (in robust part).
Loadings matrix before applying sparse reweighting step, a numeric matrix of size \(p\) by \(k\).
Numeric vector containing the indices of the variables that are used in the sparse reweighting step.
Logical vector of size \(n\) indicating if an observation is kept in the sparse reweighting step.
Loadings matrix before estimating eigenvalues, a numeric matrix of size \(p\) by \(k\).
Logical vector of size \(n\) indicating if an observation is kept in the final SD reweighting step.
The robustness parameter \(\alpha\) used throughout the algorithm.
The \(h\)-parameter used throughout the algorithm.
Numeric vector of size \(n\) containing the robust score distances within the robust PCA subspace.
Numeric vector of size \(n\) containing the orthogonal distances to the robust PCA subspace.
Cut-off value for the robust score distances.
Cut-off value for the orthogonal distances.
Numeric vector of size \(n\) containing the SD-flags of the observations. The observations whose score distance is larger than cutoff.sd receive an SD-flag equal to zero. The other observations receive an SD-flag equal to 1.
Numeric vector of size \(n\) containing the OD-flags of the observations. The observations whose orthogonal distance is larger than cutoff.od receive an OD-flag equal to zero. The other observations receive an OD-flag equal to 1.
Numeric vector of size \(n\) containing the flags of the observations. The observations whose score distance is larger than cutoff.sd or whose orthogonal distance is
larger than cutoff.od can be considered as outliers and receive a flag equal to zero.
The regular observations receive flag 1.
An \(n\) by \(p\) matrix or data matrix with observations in the rows and variables in the columns.
Number of principal components that will be used.
Maximal number of principal components that will be computed, default is 10.
Robustness parameter, default is 0.75.
The number of outliers the algorithm should resist is given by \(n-h\). Any value for h between \(n/2\) and \(n\) may be specified. Default is NULL which uses h=ceiling(alpha*n)+1. Do not specify alpha and h at the same time.
Number of directions used when computing the outlyingness (or the adjusted outlyingness when skew=TRUE), see outlyingness and adjOutl for more details.
Logical indicating if the grid version of sparse PCA should be used (SPcaGrid with method="sd" from rrcovHD). Otherwise, the version of Zou et al. (2006) is used (spca from elasticnet). Default is TRUE.
Sparsity parameter of SPcaGrid (when grid=TRUE) or ridge parameter of spca (when
grid=FALSE), default is \(10^{-6}\).
Parameter for spca (only used when grid=FALSE), see spca for more details.
Parameter for spca (only used when grid=FALSE), see spca for more details.
If TRUE, the data are standardised robustly in the beginning and classically before applying sparse
PCA. If FALSE, the data are only mean-centred before applying sparse PCA. Default is TRUE.
Logical indicating if the version for skewed data should be applied, default is FALSE.
Tom Reynkens, based on R code from Valentin Todorov for PcaHubert in rrcov (released under GPL-3) and Matlab code from Katrien Van Driessen (for the univariate MCD).
The ROSPCA algorithm consists of an outlier detection part (step 1), and a sparsification part (steps 2 and 3). We give an overview of these steps here and refer to Hubert et al. (2016) for more details.
Step 1: This is a robustness step similar to ROBPCA. When a standardisation is appropriate, the variables are first robustly standardised by means of the componentwise median and the \(Q_n\). Using the singular value decomposition (SVD) of the resulting data matrix, the \(p\)-dimensional data space is reduced to the affine subspace spanned by the \(n\) observations. Then, the subset of the \(h\) observations with smallest outlyingness is selected (\(H_0\)). Thereafter, a reweighting step is applied: given the orthogonal distances to the preliminary PCA subspace determined by the observations in \(H_0\), all observations with orthogonal distances (ODs) smaller than the corresponding cut-off are kept (\(H_1\)).
Step 2: First, the data points with indices in \(H_1\) are standardised using the componentwise median and the \(Q_n\) and sparse PCA is applied to them. Then, an additional reweighting step is performed which incorporates information about the sparse structure of the data. Variables with zero loadings on all \(k\) PCs are discarded and then the orthogonal distances to the estimated sparse PCA subspace are computed. This yields an index set \(H_2\) of observations with orthogonal distance smaller than the cut-off corresponding to these new orthogonal distances. Thereafter, the subset of observations with indices in \(H_2\) is standardised using the componentwise median and the \(Q_n\) of the observations in \(H_1\) (the same standardisation as in the first time sparse PCA is applied) and sparse PCA is applied to them which gives sparse loadings. Adding the discarded zero loadings again gives the loadings matrix \(P_2\).
Step 3: In the last step, the eigenvalues are estimated robustly by applying the \(Q_n^2\) estimator on the scores of the observations with indices in \(H_2\). In order to robustly estimate the centre, the score distances are computed and all observations of \(H_2\) with a score distance smaller than the corresponding cut-off are considered, this is the set \(H_3\). Then, the centre is estimated by the mean of these observations. Finally, the estimates of the eigenvalues are recomputed as the sample variance of the (new) scores of the observations with indices in \(H_3\). The eigenvalues are sorted in descending order, so the order of the PCs may change. The columns of the loadings and scores matrices are changed accordingly.
Note that when it is not necessary to standardise the data, they are only centred as in the scheme above, but not scaled.
In contrast to Hubert et al. (2016), we allow for SPCA (Zou et al., 2006) to be used as the sparse PCA method inside ROSPCA (grid=FALSE). Moreover, we also include a skew-adjusted version of ROSPCA (skew=TRUE) similar to the skew-adjusted version of ROBPCA (Hubert et al., 2009). This adjusted version is not detailed in Hubert et al. (2016).
Hubert, M., Reynkens, T., Schmitt, E. and Verdonck, T. (2016). ``Sparse PCA for High-Dimensional Data with Outliers,'' Technometrics, 58, 424--434.
Hubert, M., Rousseeuw, P. J., and Vanden Branden, K. (2005), ``ROBPCA: A New Approach to Robust Principal Component Analysis,'' Technometrics, 47, 64--79.
Hubert, M., Rousseeuw, P. J., and Verdonck, T. (2009), ``Robust PCA for Skewed Data and Its Outlier Map," Computational Statistics & Data Analysis, 53, 2264--2274.
Croux, C., Filzmoser, P., and Fritz, H. (2013), ``Robust Sparse Principal Component Analysis,'' Technometrics, 55, 202--214.
Zou, H., Hastie, T., and Tibshirani, R. (2006), ``Sparse Principal Component Analysis,'' Journal of Computational and Graphical Statistics, 15, 265--286.
PcaHubert, robpca, outlyingness, adjOutl, SPcaGrid, spca
X <- dataGen(m=1, n=100, p=10, eps=0.2, bLength=4)$data[[1]]
resRS <- rospca(X, k=2, lambda=0.4, stand=TRUE)
diagPlot(resRS)
Run the code above in your browser using DataLab