cellPcaCoDa: Cellwise-robust PCA for compositional data
In robCompositions: Compositional Data Analysis

cellPcaCoDa

R Documentation

Cellwise-robust PCA for compositional data

Description

Performs principal component analysis for compositional data that is robust against cellwise contamination in raw parts. The method accounts for the dependent contamination structure induced by closure on the simplex.

Usage

cellPcaCoDa(
  x,
  k = NULL,
  method = "adaptive",
  rho = "tukey",
  alpha = NULL,
  maxiter = 100,
  tol = 1e-06,
  reg = 0,
  w0 = 0.1,
  alpha_detect = 0.01,
  trace = FALSE,
  init = c("mcd", "classical")
)

Arguments

`x`	compositional data with at least three parts (a data.frame or matrix of strictly positive values whose rows are compositions, D >= 3)
`k`	number of principal components to retain. Cellwise detection needs a non-empty residual space, so use `k < D - 1`; with `k = D - 1` the reconstruction is exact and no cells are flagged. If `NULL` (default), `k` is chosen by the Kaiser criterion (eigenvalues above the average) — an implementation convenience, not part of Templ (2026); set `k` explicitly (e.g. from a scree plot) for publication results.
`method`	contamination detection method: currently only `"adaptive"` is implemented. A `"fixed"` method is planned for a future release.
`rho`	loss function: `"tukey"` (default, Tukey biweight, as in Templ 2026) or `"huber"`.
`alpha`	tuning constant for the loss. If `NULL`, 4.685 (Tukey) or 1.345 (Huber) is used, each giving 95% Gaussian efficiency.
`maxiter`	maximum number of alternating iterations (default: 100)
`tol`	convergence tolerance on the relative change in the objective (default: 1e-6)
`reg`	ridge parameter added to the diagonal of the weighted covariance for ill-conditioned high-dimensional compositions (default: 0, i.e. the published algorithm). An implementation extension; see Details.
`w0`	floor weight assigned to contaminated cells (default: 0.1). Must be in (0, 1).
`alpha_detect`	Bonferroni significance level for the cellwise contamination detection test (default: 0.01).
`trace`	logical; if `TRUE`, iteration progress is printed
`init`	initialisation of the alternating scheme: `"mcd"` (default) uses a robust MCD fit of the ilr coordinates (the default in Templ 2026, underpinning its breakdown results); `"classical"` uses a fast classical SVD start. Falls back to `"classical"` (with a warning) if `covMcd` fails.

Details

Compositional data are first expressed in isometric logratio (pivot) coordinates via pivotCoord. An alternating optimisation scheme then jointly estimates cell-level contamination flags in the original raw-part space, maps them to coordinate-level weights through the contrast matrix, and computes a weighted PCA in the ilr space. Loadings and the robust centre are back-transformed to the clr space for interpretation, following the same convention as pcaCoDa.

By default the alternating scheme is initialised from a robust MCD fit of the ilr coordinates (init = "mcd", the default in Templ 2026, underpinning its breakdown results); a fast classical SVD start is available via init = "classical". Robustness is controlled by a bounded loss rho, which defaults to the Tukey biweight (95% Gaussian efficiency), as in the paper.

The detection step uses a projection-based test derived from the propagation theorem: if raw part j alone is contaminated, the perturbation in ilr space lies along the direction d_j = V^\top (e_j - 1/D) (normalised to unit length). MAD-standardised ilr residuals are projected onto each d_j, and the projection score is standardised by its correlation-adjusted scale \sigma_{p,j} = \sqrt{d_j^\top R\, d_j}, where R is the correlation matrix of the standardised residuals, so the standardised score is approximately N(0,1) under the null. A cell is flagged when this score exceeds the Bonferroni-adjusted threshold z_{1 - \alpha/(2D)}.

The weight mapping from raw-part flags to ilr weights is fully vectorised via a matrix multiply (no triple-nested loop). For each observation i and ilr coordinate l the weight is

w_{il} = \max\Bigl(w_0,\; \exp\bigl(\sum_j c_{ij}\,\log\bigl(1-(1-w_0)\,|V_{jl}|/\max_j |V_{jl}|\bigr)\bigr)\Bigr)

where c_{ij} are the binary contamination flags (1 = flagged).

The weighted covariance uses per-coordinate normalisation: entry (l, m) of the covariance is normalised by \sqrt{\sum_i w_{il} \cdot \sum_i w_{im}} so that coordinates with different total weights are treated correctly.

Implementation choices beyond Templ (2026). Two conveniences are not part of the published algorithm and exist only for ease of use: (i) when k is NULL it is chosen by the Kaiser criterion (retain eigenvalues above the average) — the paper leaves k to the analyst, so set it explicitly (e.g. from a scree plot or a clear eigenvalue gap) for publication-quality results; and (ii) the ridge parameter reg (default 0, i.e. the paper's algorithm) adds a small diagonal term for ill-conditioned high-dimensional compositions, a regime the paper lists as future work. Convergence is declared when the objective's relative change falls below tol, or when a period-2 limit cycle is detected (a few borderline cells oscillating under the hard-threshold reweighting); in the latter case it stops at the stabilised cycle (the solution is stable up to those few cells) and cycled is TRUE.

Value

An object of class "cellPcaCoDa" with components:

`scores`	n x k score matrix in ilr space
`loadings`	D x k loadings in clr space
`eigenvalues`	the first k eigenvalues of the weighted covariance
`eigenvalues_all`	all eigenvalues of the final weighted covariance (length D-1), for explained-variance reporting
`cellflags`	n x D logical matrix of contamination flags in raw-part space
`cellweights`	n x (D-1) matrix of cell weights in ilr space
`center`	robust centre in ilr space (length D-1)
`center_clr`	robust centre in clr space (length D)
`converged`	logical: did the algorithm converge within `maxiter` iterations?
`iterations`	number of iterations actually used
`method`	the detection method used
`k`	number of retained components
`rho`	the loss function name
`init`	the initialisation actually used (`"mcd"` or `"classical"`)
`cycled`	logical: `TRUE` if the algorithm stopped at a period-2 limit cycle (the loadings are stable up to a few oscillating borderline cells) rather than at a fixed point; `converged` is then also `TRUE`

Note

cellPcaCoDa targets cellwise contamination (individual corrupted parts); use pcaCoDa when whole observations are outlying (rowwise/casewise contamination). Inspect colMeans(res$cellflags) for the per-part contamination rate and rowSums(res$cellflags) for the per-observation count. The cell-weight floor w0 is stable on [0.05, 0.2]; for heavier contamination (rate \ge 0.15) use [0.05, 0.08]. The detection level alpha_detect (default 0.01) is Bonferroni-adjusted across the D parts. The input must be strictly positive (no zeros); replace rounded or structural zeros first, e.g. with imputeBDLs. At least three parts are required (D \ge 3), and cellwise detection needs k < D-1 (at k = D-1 the residual space is empty and no cells are flagged). The method and its defaults assume a fixed number of parts with growing sample size.

Author(s)

Matthias Templ

References

Templ, M. (2026). cellPcaCoDa: cellwise-robust principal components for compositional data. Advances in Data Analysis and Classification. In press.

Templ, M. (2026). Log-ratio propagation on the simplex: a theory of cellwise contamination for compositional data. arXiv 2605.31345.

Examples

## --- a 5-part composition: detection enabled (k = 2 < D-1) -------------
data(expenditures)
set.seed(123)
x <- expenditures
x[1, 2] <- x[1, 2] * 10      # contaminate one cell
x <- constSum(x)             # re-close to the simplex
res <- cellPcaCoDa(x, k = 2)
res                          # print method
head(res$cellflags)          # logical n x D contamination flags


plot(res)                    # contamination heatmap
plot(res, which = 2)         # clr biplot

## --- paper application: GEMAS geochemistry (Templ 2026, Sec. 5) --------
data(gemas)
elements <- intersect(c("Al", "Ca", "Fe", "K", "Mg", "Mn", "Na"),
                      colnames(gemas))
xg <- gemas[, elements]
xg <- xg[complete.cases(xg), ]
xg <- constSum(xg, const = 100)               # close to 100%
resg <- cellPcaCoDa(xg, k = 3)                 # 7 parts, k = 3 < D-1
## per-element cellwise contamination rate
round(sort(colMeans(resg$cellflags), decreasing = TRUE), 3)
plot(resg)                   # cellwise flags / clr biplot

robCompositions documentation built on Nov. 1, 2026, 1:06 a.m.