Principal Subsimplex Analysis

Compositional data, which are also referred to as simplicial data, are multivariate observations consisting of vectors of proportions. Such data are prevalent in domains where the relative magnitude between variables is the primary concern, including fields such as geochemistry (proportions of constituent elements), microbiology (proportions of species), and economics (proportions of portfolio components). Compositional data vectors are characterized by the constraints that the entries are nonnegative and sum to one. This constraint complicates the application of traditional multivariate analysis techniques, which are designed for Euclidean data, to compositional datasets.

A major challenge with compositional data is identifying meaningful lower dimensional approximations and the corresponding modes of variation (see Section 3.1.4. of Marron and Dryden, (2021).) A mode of variation is a set of data objects (e.g. composition vectors) that is one-dimensional in some sense which describes a component of variation. A prototypical example is found in classical Principal Components Analysis (PCA), where modes of variations are the straight lines through the mean spanned by the loading vectors.

While classical PCA is popular, it often fails to provide effective approximations and modes of variation for non-Euclidean data, including compositional data. Figure 1 demonstrates this issue using a compositional dataset with three variables. The distribution of this 2-dimensional compositional data set is effectively visualized by a ternary plot, a rotation of the 2-simplex, as shown in Figure 1. The blue dots represent data points, while the green arrows depict the first and second principal component directions. The red dots show the one-dimensional approximations of the two selected blue data points, as indicated by the dashed lines. Notably, these red dots fall outside the original compositional space, meaning that the lower-dimensional representations have negative proportions and no longer represents compositions. This example motivates development of more interpretable and effective methods for lower-dimensional approximation in compositional data analysis.

A popular strategy to address such challenges is applying transformations such as log-ratio transformations (Aitchison, 1982, 1983, 1986) and power transformations (Aitchison, 1986) followed by the usual PCA. Log-ratio transformations bijectively map the interior of the simplex to a Euclidean space of the same dimension, thus enabling application of techniques for Euclidean data. Log-ratio transformations well capture curvature present in compositional data and are the basis of the well-known logistic normal distribution (Aitchison, 1986). On the other hand, power transforms provide a natural continuum between no transformation and log-ratio transformations. Power transforms are effective in normalizing marginal distributions of data, especially remedying high skewness which is common in compositional data. However, these transformations can strongly distort the original compositional space and lead to undesirable results in some applications. Log-ratio transformations magnify variation within variables with low average proportions (minor variables), obscuring important signals present in variables with high average proportions (major variables). In addition, log-ratio transformations do not naturally accommodate data points with zero proportions, requiring either a separate treatment of those data points or imputation for zeros. Scealy et al., (2015) proposed a different set of power transformations that map compositional data vectors onto the surface of various manifolds. They then applied classical and robust PCA on the tangent space to the manifold. All power transformations inherit the problem of direct PCA in the spirit of Figure 1, that the lower dimensional approximations may leave the object space, and the corresponding modes of variation are not interpretable in a compositional sense.

The above discussion motivates the development of an alternative PCA method that aims to describe variation in major variables while avoiding the limitations of transformations. In other contexts there has also been a recent growing shift away from transformation based approaches, with a focus on analyzing compositional data on its original scale of measurement (Xiong et al., 2015, Weistuch et al., 2022, Fiksel et al., 2022, Scealy and Wood, 2023, Firth and Sammut, 2023, Lundborg and Pfister, 2023, Scealy et al., 2024). This article aligns with this view.

This paper proposes applying the backwards Principal Component Analysis approach to compositional data. Backwards PCA was motivated by Jung et al., (2012) then discussed in detail by Damon and Marron, (2014) and Marron and Dryden, (2021, Section 8.6). Given an object space that is a rank $d$ manifold, backwards PCA sequentially finds the best-fitting rank $r$ submanifold of the previous rank $r+1$ submanifold. A major advantage of the backwards PCA is that its approximations of any rank naturally remain in the object space, and the nested rank $r$ approximations have natural relationships with the rank $r+1$ approximations which lead to easily defined $r$th scores. The nested nature of backwards PCA naturally leads to modes of variation, which distinguishes backwards PCA from general dimensionality reduction methods.

Backwards PCA is a unifying framework and requires special attention for each class of object spaces. Backwards PCA has been implemented for a variety of object spaces, including spherical data (Jung et al., 2012), skeletal representations (Pizer et al., 2013), polyspheres (Eltzner et al., 2015), nonnegative data (Zhang et al., 2015), and high-dimensional tori (Eltzner et al., 2018, Zoubouloglou et al., 2023). More discussion on the applications of backwards PCA can be found in Section S1 of the supplementary material.

Using the backwards PCA framework, this paper proposes Principal Subsimplex Analysis (PSA). Given a $d$ dimensional compositional data set, PSA effectively identifies a nested sequence of simplices of dimension $r$ for $r=d-1,\cdots,1,0$. The vertices of each subsimplex form a partition of parts, thus PSA is closely related to amalgamation in the sense defined in Aitchison, (1986). We introduce two versions: PSA via Simplices (PSA-S) and PSA via Orthants (PSA-O), which use the Euclidean metric and the spherical metric (arc length on the unit sphere), respectively. It will become clear that PSA-S approximates data through amalgamation, while PSA-O combines parts in a related but distinct manner.

Our proposed method has the following four key features.

1. 1.

   Targets Variation among Major Variables: Because PSA uses Euclidean or spherical metrics, it has better ability than transformation-based approaches to capture important variation that occurs among major variables. In addition, PSA is more robust to the existence of pure noise variables.

2. 2.

   Natural Treatment of Zeros: By construction of the method, PSA naturally handles data sets with zeros.

3. 3.

   Compositional Lower Dimensional Approximations: Each lower dimensional approximation is a simplex whose vertices represent groups of variables, which offers interpretable structure. This approach preserves the important nature of compositional data, such as nonnegativity and potential negative correlations between variables. In particular, the 2-dimensional representation provides useful visualization.

4. 4.

   Linear Modes of Variation: PSA produces interpretable modes of variation. In addition, these modes of variation are linear, which serve as a more interpretable alternative to the non-linear modes of variation produced by log-ratio PCA.

Quinn and Erb, (2020) proposed a method Amalgams which uses amalgamation for dimensionality reduction of compositional data. For $d$-dimensional compositional data and a prescribed number $k<d$, Amalgams searches for a $k$-dimensional amalgamation that optimizes a given objective function, for example, either the average inter-sample distance or classification accuracy. The most important difference between our method and Amalgams is that our method, based on backwards PCA, produces modes of variation. Our method may also be computationally faster as Amalgams employed a slow genetic algorithm for estimation.

The rest of the paper is organized as follows. Section 2 describes the geometry of simplices and proposes PSA-S. Motivated by scaling challenges, Section 3 proposes PSA-O. These two versions of PSA, together with Euclidean PCA and the two transformation-based PCAs, are applied to simulated data sets in Section 5, illustrating distinctive properties of our new PSA methods. The same methods are compared on relative abundances of diatom species during the late Pliocene in Section 6. The analysis for this article was performed using R Statistical Software (R Core Team,, 2024). The R package PSA and the authors’ code for this article are available at https://github.com/haneone33/Principal-Subsimplex-Analysis.

In a Euclidean space $\mathbb{R}^{d}$, a set $\boldsymbol{V}=\{\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{r+1}\}$ of $r+1$ vectors is said to be affinely independent if $(\boldsymbol{v}_{1}-\boldsymbol{v}_{r+1}),\cdots,(\boldsymbol{v}_{r}-\boldsymbol{v}_{r+1})$ are linearly independent. An $r$-dimensional simplex is the convex hull $S_{r}$ of an affinely independent set of $r+1$ vectors, that is,

In this case, the vectors $\{\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{r+1}\}$ are called the vertices of $S_{r}$. Given a simplex $S$, a subsimplex of $S$ is a simplex of any dimension that is also a subset of $S$.

The $d$-dimensional unit simplex, denoted by $\Delta_{d}$, is a special case of a simplex in $\mathbb{R}^{d+1}$ whose vertices are the unit vectors $\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{d+1}$, where $\boldsymbol{e}_{j}$ is a vector of zeros with a unique one in the $j$th coordinate. A vector of compositions of $d+1$ elements naturally belongs to the $d$-dimensional unit simplex. Given a compositional data set in $\Delta_{d}$, PSA reveals interpretable lower dimensional representations taking the form of subsimplices of $\Delta_{d}$ of varying dimensions.

The sequence of nested subsimplices (1) has the following explicit representation using vertices,

where $\left\{\boldsymbol{v}_{1}^{(r-1)},\cdots,\boldsymbol{v}_{r}^{(r-1)}\right\}$ is an affinely independent subset of $S_{r+1}\left(\boldsymbol{v}_{1}^{(r)},\cdots,\boldsymbol{v}_{r+1}^{(r)}\right)$. For notational consistency, we set $\boldsymbol{v}_{j}^{(d)}=\boldsymbol{e}_{j},j=1,\cdots,d+1$ so that

Given an $r$-simplex $S_{r}=S_{r}\left(\boldsymbol{v}_{1}^{(r)},\cdots,\boldsymbol{v}_{r+1}^{(r)}\right)$, an easy and interpretable way to construct $S_{r-1}\left(\boldsymbol{v}_{1}^{(r-1)},\cdots,\boldsymbol{v}_{r}^{(r-1)}\right)$ is to form the vertex set by aggregating two of the vertices of $S_{r}$ at a ratio $\alpha_{r}\in[0,1]$, while keeping the other $r-1$ vertices. Without loss of generality, suppose the two merging vertices are $\boldsymbol{v}^{(r)}_{1}$ and $\boldsymbol{v}^{(r)}_{2}$. The new vertex $\boldsymbol{v}_{1}^{(r-1)}$ is given as the convex combination of these two vertices at an appropriate ratio $\alpha_{r}\in[0,1]$,

This new vertex together with the existing vertices after the reindexing

are affinely independent give an $(r-1)$ subsimplex $S_{r-1}\left(\boldsymbol{v}_{1}^{(r-1)},\cdots,\boldsymbol{v}_{r}^{(r-1)}\right)\subset S_{r}$.

Given an $r$-dimensional simplex $S_{r}$, the size of an $(r-1)$-dimensional subsimplex varies depending on its vertex set, in terms of its $(r-1)$-dimensional volume. In particular, a $(d-1)$-dimensional subsimplex of the unit $d$-simplex $\Delta_{d}$ is smaller than the unit $(d-1)$-simplex unless the vertices of the subsimplex are chosen from the vertices of $\Delta_{d}$. This raises a challenge in terms of variation explained by the Residual Sum of Squared scores (RSS). This scale problem is handled by correctly identifying an $r$-simplex with the unit $r$-simplex.

An $r$-simplex $S_{r}$ is naturally identified with the unit $r$-simplex by matching the vertices of $S_{r}$ to those of the unit $r$-simplex. We define the scaling function $\phi_{S_{r}}$ by

Given a sample $\boldsymbol{x}=\sum_{j=1}^{r+1}\alpha_{j}\boldsymbol{v}_{j}$, the proportion associated with vertex $\boldsymbol{v}_{j}$ is given by $\alpha_{j}$. In this sense, $\phi_{S_{r}}$ can be interpreted as writing $\boldsymbol{x}$ in a new coordinate system with basis $\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{r+1}$.

To describe the PSA-S procedure, assume we have a $d$-dimensional compositional data set $\left\{\boldsymbol{x}_{1},\cdots,\boldsymbol{x}_{n}\right\}$ in $\Delta_{d}$. We denote by $\hat{\boldsymbol{x}}^{(r)}_{i}$ the rank $r$ approximation of $\boldsymbol{x}_{i}$ which is a member of $S_{r}\left(\boldsymbol{v}^{(r)}_{1},\cdots,\boldsymbol{v}^{(r)}_{r+1}\right)$, for $r=d-1,\cdots,0$. For notational consistency, we set $\hat{\boldsymbol{x}}^{(d)}_{i}=\boldsymbol{x}_{i},i=1,\cdots,n$.

The subsimplices in (3) are constructed inductively through repeated aggregation of vertices. Suppose we have obtained the rank $r$ subsimplex $S_{r}=S_{r}\left(\boldsymbol{v}^{(r)}_{1},\cdots,\boldsymbol{v}^{(r)}_{r+1}\right)$ and the approximations $\hat{\boldsymbol{x}}^{(r)}_{i}$. Further suppose that the approximating subsimplex $S_{r-1}=S_{r-1}\left(\boldsymbol{v}^{(r-1)}_{1},\cdots,\boldsymbol{v}^{(r-1)}_{r}\right)$ is formed by combining $\boldsymbol{v}^{(r)}_{1}$ and $\boldsymbol{v}^{(r)}_{2}$ at a ratio $\alpha_{r}$. (The first two vertices are chosen for notational convenience.) PSA-S approximates a data point $\hat{\boldsymbol{x}}^{(r)}=\sum_{j=1}^{r+1}\hat{x}^{(r)}_{j}\boldsymbol{v}^{(r)}_{j}$ in $S_{r}$ in a mass preserving way, as illustrated in Figure 2. The point $\hat{\boldsymbol{x}}^{(r)}$ in $S_{r}$ is mapped onto $S_{r-1}$ along the line segment parallel to the edge connecting $\boldsymbol{v}_{1}^{(r)}$ and $\boldsymbol{v}_{2}^{(r)}$. This mapping preserves the weights $x_{3},\cdots,x_{r+1}$ and only re-distributes the first two weights according to the ratio $\alpha_{r}$. In other words, the rank $r-1$ approximation is given by the map $\pi^{PSA-S}_{r}:S_{r}\left(\boldsymbol{v}_{1}^{(r)},\cdots,\boldsymbol{v}_{r+1}^{(r)}\right)\to S_{r-1}\left(\boldsymbol{v}_{1}^{(r-1)},\cdots,\boldsymbol{v}_{r}^{(r-1)}\right)$,

Thus the rank $(r-1)$ approximation is $\hat{\boldsymbol{x}}^{(r-1)}_{i}=\pi^{PSA-S}_{r}\left(\hat{\boldsymbol{x}}^{(r)}_{i}\right)$. Notably, the mapping is well defined for data points with zero proportions.

The $r$th score of a point $\boldsymbol{x}_{i}=\sum_{j=1}^{r+1}\hat{x}^{(r)}_{i,j}\boldsymbol{v}^{(r)}_{j}$ is the signed distance from $\hat{\boldsymbol{x}}^{(r)}_{i}$ to $\hat{\boldsymbol{x}}^{(r-1)}_{i}$ where $S_{r}$ is rescaled to have unit size through $\phi_{S_{r}}$:

To arrive at an analogue of modes of variation in PCA, define the $r$th loading vector $\boldsymbol{l}_{r}$ to be the difference of the two merged vertices, i.e. $\boldsymbol{l}_{r}=\boldsymbol{v}^{(r)}_{2}-\boldsymbol{v}^{(r)}_{1}$. This loading vector indicates the direction pointing from each $\hat{\boldsymbol{x}}^{(r)}_{i}$ to $\hat{\boldsymbol{x}}^{(r-1)}_{i}$. The score $s^{(r)}_{i}$ represents the amount of the movement in that direction because

As described in Section 3 of the supplementary material, the corresponding mode of variation is $\left\{\boldsymbol{v}^{(0)}+\frac{s^{(r)}_{i}}{\sqrt{2}}\boldsymbol{l}_{r}:i=1,...,n\right\}$ and thus $\boldsymbol{l}_{r}$ is indeed the direction of the mode of variation. More detailed investigation of modes of variation for PSA and backwards PCA are found in the same section of the supplementary material.

The merged vertices and the ratio $\alpha_{r}$ are chosen so that they minimize the Residual Sum of Squared scores (RSS) $\sum_{i=1}^{n}\left(s^{(r)}_{i}\right)^{2}$. Given the vertices $\boldsymbol{v}_{1}$ and $\boldsymbol{v}_{2}$ to be merged, the optimal $\alpha_{r}$ has a closed-form solution

The pair of vertices are chosen so that RSS is minimized at the optimal $\alpha_{r}$.

The procedure of PSA-S is summarized in Algorithm 1 of the Supplementary material.

Let $\boldsymbol{V}=\left\{\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{r+1}\right\}$ be a set of $r+1$ orthonormal vectors in $\mathbb{R}^{d}$. The $r$-dimensional nonnegative orthant with the vertex set $\boldsymbol{V}$ is the set $O_{r}$ of linear combinations of $\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{r+1}\in\mathbb{R}^{r}$ with nonnegative coefficients whose squared sum is one. In other words,

Note that $O_{r}$ is a subset of the $(d-1)$-dimensional unit sphere. Important subsets of an $r$-nonnegative orthant $O_{r}$ are suborthants, which are subsets of $O_{r}$ that are $r^{\prime}$-nonnegative orthants for some $r^{\prime}<r$. The unit $d$-nonnegative orthant $\mathcal{O}_{d}$ is a special case of a nonnegative orthant in $\mathbb{R}^{d+1}$ whose vertices are $\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{d+1}$.

For each $r$-nonnegative orthant $O_{r}\left(\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{r+1}\right)$ there exists a corresponding $r$-simplex $S_{r}\left(\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{r+1}\right)$ with the same vertex set. Each nonnegative orthant and its corresponding simplex are intimately related through the invertible projection map $P$,

The inverse correspondence is

This projection is equivalent to transporting $\boldsymbol{x}\in S_{r}$ along the line from the origin through $\boldsymbol{x}$ until it intersects with the unit sphere.

The (geodesic) distance along the surface of the sphere $S^{d-1}$ between two points $\boldsymbol{x},\boldsymbol{y}$ in an $r$-nonnegative orthant $O_{r}\left(\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{r+1}\right)$ is the length of the shorter arc that connects two points, which is a part of a great circle. The distance function is given by $d(\boldsymbol{x,y})=\cos^{-1}(\boldsymbol{x^{\top}y})$. Equipped with this geodesic distance, any two $r$-nonnegative orthants are isometric.

PSA-O proceeds in a manner similar to that of PSA-S but the fundamental low-rank approximation happens in the unit nonnegative orthant $\mathcal{O}_{d}=O_{d}(\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{d+1})$. PSA-O has three steps: (i) project data points $\boldsymbol{x}_{i}$ in $\Delta_{d}(\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{d+1})$ to $\tilde{\boldsymbol{x}}_{i}$ in $O_{d}(\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{d+1})$ through the projection $\boldsymbol{x}\mapsto\frac{\boldsymbol{x}}{\|\boldsymbol{x}\|_{2}}$, (ii) find a nested sequence of nonnegative orthants $O_{r}(\tilde{\boldsymbol{v}}^{(r)}_{1},\cdots,\tilde{\boldsymbol{v}}^{(r)}_{r+1})$ and lower dimensional approximations $\tilde{\boldsymbol{x}}_{r}$ for $r=d,d-1,\cdots,0$, and (iii) project $O_{r}(\tilde{\boldsymbol{v}}^{(r)}_{1},\cdots,\tilde{\boldsymbol{v}}^{(r)}_{r+1})$ and $\tilde{\boldsymbol{x}}_{r}$ in $O_{d}(\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{d+1})$ back onto $\Delta_{d}(\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{d+1})$ through the inverse of the projection map, $\tilde{\boldsymbol{x}}\mapsto\frac{\tilde{\boldsymbol{x}}}{\|\tilde{\boldsymbol{x}}\|_{1}}$. As mentioned earlier, the use of nonnegative orthants naturally avoids the scaling problem because all nonnegative orthants of the same dimension have the same size. We will use $\tilde{\boldsymbol{x}}$ to indicate points on nonnegative orthants, and continue to use $\boldsymbol{x}$ and $\hat{\boldsymbol{x}}$ for points on simplices.

The above three steps are detailed as follows. Firstly, we project $\boldsymbol{x}_{i}$ on $\Delta_{d}$ to the corresponding point $\tilde{\boldsymbol{x}}_{i}^{(d)}=\frac{\hat{\boldsymbol{x}}_{i}^{(d)}}{\|\hat{\boldsymbol{x}}_{i}^{(d)}\|_{2}}$ on $O_{r}$ and initialize $\tilde{\boldsymbol{v}}_{j}^{(d)}=\boldsymbol{e}_{j}$ for $j=1,\cdots,d+1$ so that

Next, similar to the PSA-S algorithm, PSA-O iteratively combines two vertices at some ratio $\alpha_{r}\in[0,1]$. Suppose the rank $r$ approximating orthant $O_{r}\left(\tilde{\boldsymbol{v}}_{1}^{(r)},\cdots\tilde{\boldsymbol{v}}_{r+1}^{(r)}\right)$ has been computed and $\tilde{\boldsymbol{v}}_{1}^{(r)}$ and $\tilde{\boldsymbol{v}}_{2}^{(r)}$ are merged to give the new vertex

This new vertex and the remaining vertices $\tilde{\boldsymbol{v}}_{j}^{(r-1)}=\tilde{\boldsymbol{v}}_{j+1}^{(r)}$ for $j=2,\cdots,d$ form an orthonormal set and thus defines an $(r-1)$-nonnegative suborthant

The lower dimensional approximation $\tilde{\boldsymbol{x}}_{i}^{(r-1)}$ of $\tilde{\boldsymbol{x}}_{i}^{(r)}\in O_{r}\left(\tilde{\boldsymbol{v}}_{1}^{(r)},\cdots\tilde{\boldsymbol{v}}_{r+1}^{(r)}\right)$ is the point in $O_{r-1}\left(\tilde{\boldsymbol{v}}_{1}^{(r-1)},\cdots,\tilde{\boldsymbol{v}}_{r}^{(r-1)}\right)$ that is closest to $\tilde{\boldsymbol{x}}_{i}^{(r)}$ with respect to the geodesic distance (shown as an arc in Figure 3(b)) on $\mathbb{S}^{d}$ and is given by

This projection is equivalent to transporting $\tilde{\boldsymbol{x}}_{i}^{(r)}$ onto the suborthant $O_{r-1}\left(\tilde{\boldsymbol{v}}_{1}^{(r-1)},\cdots,\tilde{\boldsymbol{v}}_{r}^{(r-1)}\right)$ along the great circle that is perpendicular to the suborthant (Figure 3). The score is the signed geodesic distance,

Again, modes of variation are based on the $r$th loading vector defined as $\boldsymbol{l}_{r}=\boldsymbol{v}^{(r)}_{2}-\boldsymbol{v}^{(r)}_{1}$. Note that although the difference $\hat{\boldsymbol{x}}^{(r)}_{i}-\hat{\boldsymbol{x}}^{(r-1)}_{i}$ is not exactly parallel to the loading vector as in the case of PSA-S, the loading vectors still serve as effective representatives of the directions.

Lastly, approximating subsets $O_{r}(\tilde{\boldsymbol{v}}^{(r)}_{1},\cdots,\tilde{\boldsymbol{v}}^{(r)}_{r+1})$ and lower dimensional approximations $\hat{\boldsymbol{x}}^{(r)}_{i}$ for $r=d,d-1,\cdots,0$ are mapped onto $\Delta_{r}$ through the inverse of the projection map.

The search for the optimal pair of vertices and the ratio $\alpha_{r}$ is implemented by multiple grid searches for $\alpha_{r}$, for each of $r(r+1)/2$ pairs of vertices. Note that the grid search is fast and effective because for each pair of vertices, $\alpha_{r}$ is a one-dimensional parameter in the bounded region [0,1].

The procedure of PSA-O is summarized in Algorithm 2 of the Supplementary material.

The log transform is defined for positive values, so we denote by $\Delta_{d}^{\circ}$ the interior of $\Delta_{d}$ or the open unit $d$-simplex, that is,

Aitchison, (1986) proposed the additive log-ratio transformation, $\boldsymbol{w}_{\text{alr}}:\Delta_{d}^{\circ}\to\mathbb{R}^{d}$ defined by

Division by the last component gives the additive log-ratio transformation the nice property that the component-wise log transformation $\log\boldsymbol{x}$ does not have: it is a one-to-one mapping from the open unit $d$-simplex to the entire Euclidean space $\mathbb{R}^{d}$.

This correspondence is particularly useful in parametric analysis of compositional data with no zero entries. A random variable supported on $\Delta_{d}^{\circ}$ is said to follow logistic normal distribution with parameters $\boldsymbol{\mu}\in\mathbb{R}^{d}$ and $\boldsymbol{\Sigma}\in\mathbb{R}^{d^{2}}$, denoted by $\mathcal{L}^{d}(\boldsymbol{\mu,\Sigma})$, if $\boldsymbol{w}_{alr}(\boldsymbol{x})$ follows $N_{d}(\boldsymbol{\mu,\Sigma})$. The class of logistic normal distributions is more flexible than the traditional class of Dirichlet distributions in the sense that it can model correlation between components in a way that Dirichlet distributions cannot.

Because the additive log-ratio transformation depends on the choice of ordering of components and it may greatly affect subsequent analyses, Aitchison, (1986) further proposed the centered log-ratio transformation, $\boldsymbol{w}_{\text{clr}}:\Delta_{d}^{\circ}\to\mathbb{R}^{d+1}$, which centers each point $\boldsymbol{x}$ by its geometric mean $g(\boldsymbol{x})=\sqrt[d+1]{x_{1}\cdots x_{d+1}}$.

The centered log-ratio transformation is a one-to-one mapping from $\Delta_{d}^{\circ}$ to the $d$-dimensional subspace $H=\{\boldsymbol{w}\in\mathbb{R}^{d+1}:\boldsymbol{1}_{d+1}^{\top}\boldsymbol{w}=0\}$. One can remove the redundant dimension by isometrically identifying the subspace with $\mathbb{R}^{d}$, which leads to the isometric log-ratio transformation, $\boldsymbol{w}_{\text{ilr}}:\Delta_{d}^{\circ}\to\mathbb{R}^{d}$. The identification can be computed by

where $\boldsymbol{H}_{d+1}$ is the lower $d\times(d+1)$ submatrix of the Helmert matrix of order $(d+1)$. This refinement is required when a statistical method requires the data to be of full rank.

For principal component analysis purpose, however, there is no real difference between using clr and ilr, but clr transformation is intuitively and computationally simpler, for clr maps all data points into the hyperplane $H$ spanned by the simplex which essentially passes through the mean. Consequently, the last principal component of the clr-transformed data is the normal vector $\boldsymbol{1}_{d+1}$ and variation explained by the principal component is zero. In addition, all lower dimensional approximating subspaces and approximation points of clr-PCA lie on $H$ and these approximations can be effectively mapped back onto the open unit simplex through the inverse the of clr map.

Figure S4 illustrates an example of 2-dimensional compositional data in which log-ratio PCA is applied. In Figure S4 (a), two outlying data points are colored blue and red, while the other points are colored black. Figure S4 (b) shows the distribution of the same data after ilr transformation, along with the magenta and green line segments indicating the first and the second principal directions. These line segments are mapped back onto the simplex in Figure S4 (c). Note that significant curvature are introduced through the log-ratio transformation.

The log-ratio transformations are not applicable when a data set contains zeros. Using terminologies in Chapter 11 of Aitchison, (1986), zeros in compositional data are called essential zeros (or structural zeros) if the corresponding parts were truly not there, and rounded zeros (or count zeros) if the zeros occur because the true values were below detection level. Essential zeros may suggest existence of distinct subpopulations that can be modeled by hierarchical models, for example, zero inflated Poisson models. Poisson component is modeling count zeros, and the point mass at the zero is modeling essential zeros. Rounded zeros are typically replaced by small values through imputation.

Box-Cox power transformations are also frequently used to handle high skewness, and they generalize the log transformation in the sense that $\lim_{\alpha\to 0}\frac{\boldsymbol{x}^{\alpha}-1}{\alpha}=\log\boldsymbol{x}$. Aitchison, (1986) proposed to apply Box-Cox power transformation to the ratio of the components, which we denote by $w_{\alpha}^{(1)}:\Delta_{d}^{\circ}\to\mathbb{R}^{d}$,

for some $\alpha>0$. On the other hand, one can consider Box-Cox power transformation on the components but not on the ratio, $\boldsymbol{w}_{\alpha}^{(2)}:\Delta_{d}^{\circ}\to\mathbb{R}^{d+1}$,

This approach was extended regarding robustness in Scealy et al., (2015). In contrast to the case of the centered log-ratio transformation, the data points resulting from power transformations lie on a curved space but not a linear space.

For PCA purposes, $\alpha$ is often chosen so that it maximizes the profile log-likelihood assuming normality after power transformation. Another approach is simply using $\alpha=1/2$ and consider $\boldsymbol{x}\mapsto\boldsymbol{x}^{1/2}$. Note that this is a translation and a dilation of the Box-Cox power transformation. This transformation effectively maps a simplex onto the nonnegative orthant of the same dimension, which proves particularly advantageous for parametric analysis of compositional data. This enables the utilization of models defined on spheres (Scealy and Welsh, 2011, Zhu and Müller, 2024).

In Section 6 of the main manuscript, we investigated the distribution of the compositions of the Diatom data. Figure S5 (a) shows the depth distribution with an equally spaced grid of the unit interval for the y-axis. This shows that the depths where the samples were measured are not uniformly distributed but appear in clusters. Figure S5(b) shows the percent of nonzero observations for each of the diatom species. The bar plot suggests that the percent of nonzero observations is not proportional to the variance, and the most observations of the last 10-15 species are zero.

Figure S6(a)-(e) are the score plot matrices of the five methods. The same color scheme as in Figure 8 of the main manuscript was used and the outlier is marked by a red circle as before. Note that the diagonal panels are different from the density plots that are usual for scatter plot matrices. Each of the diagonal panels has the x-axis of the score and the y-axis of the depth. Figure S11 provides the corresponding loading vectors.

The detailed discussion on the scores and the loading vectors of PSA-O given in Section 6 is compared to those of PSA-S. Recall that the first, second, and fourth modes of variation are driven by the first phase (Depth $>$ 67), the second phase (Depth $<$ 67), and the outlier, respectively. A similar investigation into the diagonal panels of Figure S6(a) suggests that each of the modes of variation captures variation that is driven by either different time periods or the outlier. The first four modes of variation of PSA-S are associated with the transition of species during the first phase, the second phase, the earlier period of the first phase, and the later period of the second phase, respectively. In addition, the second mode of variation is also associated with the outlier. This observation can be attributed to that the first loading vector of PSA-S is similar to that of PSA-O, as represented by positive weight of F. barronii and R. antarctica and negative weight of T. inura and A. karstenii. The second loading vector of PSA-S is similar to a mixture of the second and the fourth loading vectors of PSA-O, as represented by positive weight of T. insigna and Denticulopsis. spp., and negative weight of T. inura and A. karstenii. However, importantly, the second mode of PSA-S is not exactly same as the mixture as A. karstenii indicates the opposite directions, negative in the second mode of PSA-S and positive in the mixture of two modes of PSA-O. The two methods are pointing out different yet important aspects of the outlier. PSA-S pulls out A. ingense which is the most particular species in the outlier, while PSA-O extracts most of the species that are particular in the outlier with altered weights.

The diagonal panels of Figure S6(c) suggests that the first, second, and fourth modes of variation of PCA are associated with the transition of diatom species through the entire period, the second phase, and the later part of the first phase, respectively, while the third mode is associated with the outlier. The second and the third rows of Figure S11 indicate that the first mode of PCA is a mixture of the negative of the first mode and the positive of the second mode of PSA-O, as represented by postivie weight of T. insigna and A. karstenii and negative weight of F. barronii and R. antarctica. The second mode of PCA is a mixture of the second and the third mode of PSA-O as represented by extreme weights of T. instigna and F. bohatyi. The third mode of PCA and the fourth mode of PSA-O share some species but also have some distinct species, indicating that they see different aspects of the outlier.

The first two scores of the power transform PCA show similar trends to those of PCA except that the outlier begins to stick out in the second mode of variation. Accordingly, A. ingens is introduced in the second loading vector of the power transform PCA while the first two loading vectors of the power transform PCA are mostly the same as those of PCA. On contrary, the third mode of the power transform PCA is more noisy and does not manifest clear association with any time period. The fourth mode is driven by the outlier.

Similarly, the first two scores and the loading vectors of the log-ratio PCA are close to those of PCA, while the outlier is more clearly pulled out in the second mode of variation as a result of higher weight of A. ingens in the loading vector.

An alternative way of representing loading vector is using parallel coordinate plots. Figure S12 shows parallel coordinate plots for PSA-S and PSA-O loading vectors. An advantage of using parallel coordinate plots over using bar plots is that the same species appears at the same x-axis so one can easily compare the loading vectors across different methods. For example, by comparing Figure S12(a) and S12(b) one can quickly notice that the first loading vectors of PSA-S and PSA-O are similar and the second loading vector of PSA-S is a mixture of the second and the fourth loading vectors of PSA-O. On contrary, an advantage of using bar plots is that one can easily recognize the important players in the loading vectors, particularly when there are more than a handful of variables.

We saw in Figure 10 that the rank $2$ approximations of PSA-O can be effectively visualized through a ternary plot. The same display for PSA-S is shown in Figure S14. It is clear from the ternary plot that the observations in the first phase (yellow to green) have similarly low levels of $V3$, while the the observations in the second phase (green to navy) and the outlier have varied proportion of $V3$. In addition, $V3$ consists of T. insigna together with five more variables, effectively pin-pointing the variables which drive the variation in the second phase. Similarly to the vertices of PSA-O, $V1$ is more related to the sea-ice affinity diatoms and $V2$ is more related to open-ocean diatoms.