Robust Principal Component Analysis for Compositional Tables

TL;DR

This paper introduces a robust PCA method for compositional tables that uses a special coordinate system linked to clr coefficients, facilitating interpretation and analysis of relationships between factors in compositional data.

Contribution

It proposes a new coordinate choice related to clr coefficients, enabling robust PCA for compositional tables with improved interpretability and handling of data singularities.

Findings

Enables dimension reduction in compositional tables

Facilitates interpretation of factor relationships

Handles data singularity issues effectively

Abstract

A data table which is arranged according to two factors can often be considered as a compositional table. An example is the number of unemployed people, split according to gender and age classes. Analyzed as compositions, the relevant information would consist of ratios between different cells of such a table. This is particularly useful when analyzing several compositional tables jointly, where the absolute numbers are in very different ranges, e.g. if unemployment data are considered from different countries. Within the framework of the logratio methodology, compositional tables can be decomposed into independent and interactive parts, and orthonormal coordinates can be assigned to these parts. However, these coordinates usually require some prior knowledge about the data, and they are not easy to handle for exploring the relationships between the given factors. Here we propose a special choice of coordinates with a direct relation to centered logratio (clr) coefficients, which are particularly useful for an interpretation in terms of the original cells of the tables. With these coordinates, robust principal component analysis (PCA) is performed for dimension reduction, allowing to investigate the relationships between the factors. The link between orthonormal coordinates and clr coefficients enables to apply robust PCA, which would otherwise suffer from the singularity of clr coefficients.

Tables

Table 1. Table 1: Unemployed people in thousands partitioned according to their gender and age groups [ 25 ] .

	Australia		Austria		Belgium		Canada
Age	Men	Women	Men	Women	Men	Women	Men	Women
15-24	129	111	29	25	53	43	250	178
25-39	90	85	40	35	86	83	241	192
40-54	66	68	36	27	65	52	242	188
55+	37	19	7	3	13	11	121	75

Equations

$$\mathbf{x}=\left(\begin{array}{ccc}{x}_{11}& \cdots & x_{1J}\\ ⋮& \ddots & ⋮\\ x_{I1}& \cdots & x_{IJ}\end{array}\right),x_{ij}>0,i=1,...,I,j=1,...,J,$$

$${\mathcal{S}}^D=\left\{\mathbf{x}=(x_1,...,x_D)|x_i>0,i=1,...,D,\sum _{i=1}^D{x}_i=\kappa \right\}$$

$$\mathcal{C}(\mathbf{x})=\left(\frac{\kappa \cdot x_1}{{\sum }_{i=1}^D{x}_i},\frac{\kappa \cdot x_2}{{\sum }_{i=1}^D{x}_i},...,\frac{\kappa \cdot x_n}{{\sum }_{i=1}^D{x}_i}\right)$$

$$\mathbf{x}\oplus \mathbf{y}=(x_1{y}_1,...,x_D{y}_D),\alpha \odot \mathbf{x}=(x_1^{\alpha },...,x_D^{\alpha }),$$

$${⟨\mathbf{x},\mathbf{y}⟩}_A=\frac{1}{2D}\sum _{i=1}^D\sum _{j=1}^D\ln \frac{x_i}{x_j}\ln \frac{y_i}{y_j},{\Vert \mathbf{x}\Vert }_A=\sqrt{{⟨\mathbf{x},\mathbf{x}⟩}_A},$$

$$d_A(\mathbf{x},\mathbf{y})={\Vert \mathbf{x}\ominus \mathbf{y}\Vert }_A,$$

$$\mathrm{clr}(\mathbf{x})=\left(\ln \frac{x_1}{g(\mathbf{x})},\ln \frac{x_2}{g(\mathbf{x})},...,\ln \frac{x_D}{g(\mathbf{x})}\right),$$

$$\mathbf{z}=\mathrm{ilr}(\mathbf{x})=\left({⟨\mathbf{x},{\mathbf{e}}^1⟩}_A,{⟨\mathbf{x},{\mathbf{e}}^2⟩}_A,...,{⟨\mathbf{x},{\mathbf{e}}^{D-1}⟩}_A\right),$$

$$\mathrm{clr}(\mathbf{x})=\mathbf{Vz}=[\mathrm{clr}({\mathbf{e}}_1{)}^T,\mathrm{clr}({\mathbf{e}}_2{)}^T,...,\mathrm{clr}({\mathbf{e}}_{D-1}{)}^T]\cdot \mathrm{ilr}(\mathbf{x}).$$

$$z_i^{(l)}=\sqrt{\frac{D-i}{D-i+1}}\ln \frac{x_i^{(l)}}{\sqrt[D-i]{{\prod }_{j=i+1}^D{x}_j^{(l)}}},$$

$$\mathrm{clr}(x_l^{(l)})=\sqrt{\frac{D}{D-1}}{z}_1^{(l)}.$$

$$\mathbf{C}={\mathbf{GLG}}^T,$$

$${\mathbf{X}}^{\ast }=(\mathbf{X}-{1\mathbf{t}}^T)\mathbf{G},$$

$${\mathbf{z}}_i^{\ast }={\mathbf{G}}^T({\mathbf{z}}_i-\mathbf{t}).$$

$$\mathcal{C}(\mathbf{x})=\left(\begin{array}{ccc}\frac{\kappa x_{11}}{{\sum }_{i,j}{x}_{ij}}& \cdots & \frac{\kappa x_{1J}}{{\sum }_{i,j}{x}_{ij}}\\ ⋮& \ddots & ⋮\\ \frac{\kappa x_{I1}}{{\sum }_{i,j}{x}_{ij}}& \cdots & \frac{\kappa x_{IJ}}{{\sum }_{i,j}{x}_{ij}}\end{array}\right)$$

$$\mathbf{x}\oplus \mathbf{y}=\left(\begin{array}{ccc}{x}_{11}{y}_{11}& \cdots & x_{1J}{y}_{1J}\\ ⋮& \ddots & ⋮\\ x_{I1}{y}_{I1}& \cdots & x_{IJ}{y}_{IJ}\end{array}\right),\alpha \odot \mathbf{x}=\left(\begin{array}{ccc}{x}_{11}^{\alpha }& \cdots & x_{1J}^{\alpha }\\ ⋮& \ddots & ⋮\\ x_{I1}^{\alpha }& \cdots & x_{IJ}^{\alpha }\end{array}\right),$$

$${⟨\mathbf{x},\mathbf{y}⟩}_A=\frac{1}{2IJ}\sum _{i,j}\sum _{k,l}\ln \frac{x_{ij}}{x_{kl}}\ln \frac{y_{ij}}{y_{kl}}.$$

$$\mathbf{x}={\mathbf{x}}_{ind}\oplus {\mathbf{x}}_{int}.$$

$${\mathrm{row}}^{\perp }(\mathbf{x})=\left(\begin{array}{ccc}g(x_{11},...,x_{1J})& \cdots & g(x_{11},...,x_{1J})\\ \cdots & \cdots & \cdots \\ g(x_{I1},...,x_{IJ})& \cdots & g(x_{I1},...,x_{IJ})\end{array}\right),$$

$${\mathrm{col}}^{\perp }(\mathbf{x})=\left(\begin{array}{ccc}g(x_{11},...,x_{I1})& ⋮& g(x_{1J},...,x_{IJ})\\ ⋮& ⋮& ⋮\\ g(x_{11},...,x_{I1})& ⋮& g(x_{1J},...,x_{IJ})\end{array}\right),$$

$$x_{ij}^{ind}={\left(\prod _{k=1}^I\prod _{l=1}^J{x}_{kj}{x}_{il}\right)}^{\frac{1}{IJ}}\propto {\left(\prod _{k=1}^I{x}_{kj}\right)}^{\frac{1}{I}}{\left(\prod _{l=1}^J{x}_{il}\right)}^{\frac{1}{J}},$$

$$x_{ij}^{int}={\left(\prod _{k=1}^I\prod _{l=1}^J\frac{x_{ij}}{x_{kj}{x}_{il}}\right)}^{\frac{1}{IJ}}.$$

$$z_i^r=\sqrt{\frac{(I-i)J}{1+I-i}}\ln \frac{g({\mathbf{x}}_{i\bullet })}{[g({\mathbf{x}}_{i+1\bullet }),...,g({\mathbf{x}}_{I\bullet }){]}^{1/(I-i)}},i=1,...,I-1,$$

$$z_j^c=\sqrt{\frac{I(J-j)}{1+J-j}}\ln \frac{g({\mathbf{x}}_{\bullet j})}{[g({\mathbf{x}}_{\bullet j+1}),...,g({\mathbf{x}}_{\bullet J}){]}^{1/(J-j)}},j=1,...,J-1,$$

$$z_{rs}^{OR}=\sqrt{\frac{1}{(I-r)(J-s)(I-r+1)(J-s+1)}}\ln \prod _{i=r+1}^I\prod _{j=s+1}^J\frac{x_{ij}{x}_{rs}}{x_{is}{x}_{rj}}.$$

$$\mathrm{clr}({\mathbf{x}}_{ind}{)}_{kl}=\sqrt{\frac{I-1}{IJ}}{z}_1^{r(k)}+\sqrt{\frac{J-1}{IJ}}{z}_1^{c(l)},$$

$$\mathrm{clr}({\mathbf{x}}_{int}{)}_{kl}=\sqrt{\frac{(I-1)(J-1)}{IJ}}{z}_{11}^{OR(kl)},$$

$$\mathrm{clr}({\mathbf{x}}_{ind}{)}_{ij}=\ln \frac{x_{ij}^{ind}}{g({\mathbf{x}}_{\bullet \bullet }^{ind})},\mathrm{clr}({\mathbf{x}}_{int}{)}_{ij}=\ln \frac{x_{ij}^{int}}{g({\mathbf{x}}_{\bullet \bullet }^{int})},$$

$$\mathrm{clr}({\mathbf{x}}_{ind}{)}_{ij}=\ln \frac{g({\mathbf{x}}_{i\bullet })g({\mathbf{x}}_{\bullet j})}{g({\mathbf{x}}_{\bullet \bullet }{)}^2},\mathrm{clr}({\mathbf{x}}_{int}{)}_{ij}=\ln \frac{x_{ij}g({\mathbf{x}}_{\bullet \bullet })}{g({\mathbf{x}}_{i\bullet })g({\mathbf{x}}_{\bullet j})},$$

$$\mathrm{clr}({\mathbf{e}}^r)=\{\begin{array}{ll}\sqrt{\frac{I-i}{(I-i+1)J}}& \text{for the elements in pivot row }i,\\ -\sqrt{\frac{1}{(I-i+1)J(I-i)}}& \text{for the elements in rows }i+1,...,I,\\ 0& \text{otherwise},\end{array}$$

Full text

Robust Principal Component Analysis

for Compositional Tables

J. Rendlová1, K. Hron1, K. Fačevicová1 and P. Filzmoser2

1 Department of Mathematical Analysis and Applications of Mathematics,

Faculty of Science, Palacký University, Olomouc, Czech Republic;

2 Institute of Statistics and Mathematical Methods in Economics,

Vienna University of Technology, Vienna, Austria.

Abstract

A data table which is arranged according to two factors can often be considered as a compositional table. An example is the number of unemployed people, split according to gender and age classes. Analyzed as compositions, the relevant information would consist of ratios between different cells of such a table. This is particularly useful when analyzing several compositional tables jointly, where the absolute numbers are in very different ranges, e.g. if unemployment data are considered from different countries. Within the framework of the logratio methodology, compositional tables can be decomposed into independent and interactive parts, and orthonormal coordinates can be assigned to these parts. However, these coordinates usually require some prior knowledge about the data, and they are not easy to handle for exploring the relationships between the given factors.

Here we propose a special choice of coordinates with a direct relation to centered logratio (clr) coefficients, which are particularly useful for an interpretation in terms of the original cells of the tables. With these coordinates, robust principal component analysis (PCA) is performed for dimension reduction, allowing to investigate the relationships between the factors. The link between orthonormal coordinates and clr coefficients enables to apply robust PCA, which would otherwise suffer from the singularity of clr coefficients.

Keywords: Compositional data; Compositional table; Robust principal component analysis; Independence table; Interaction table; Pivot coordinates.

1 Introduction

Many practical data sets (e.g. in economics [11, 12]; biology [18, 5]; or sociology [7, 27]) consist of observations which contain inherently relative information about the distribution according to two factors. From a mathematical perspective, this leads to a two-factorial extension of vector compositional data [1, 30] carrying information about a relationship between and within these (row and column) factors. In [7], it was shown that such a structure, called a compositional table $\boldsymbol{x}$,

[TABLE]

can be decomposed into its independent and interactive parts where their separate analysis can be advantageous for further interpretation concerning both factors and their relationships. Similar as compositional data, compositional tables are driven by the Aitchison geometry [28] which inhibits a direct application of multivariate statistical tools. Instead, an appropriate coordinate representation needs to be established first to map the information from the Aitchison geometry into the real space. Such a coordinate system and a decomposition into independent and interactive parts were proposed in [12].

A common primary task in multivariate statistics is to reduce the dimensionality of the data at hand, done using principal component analysis (PCA). As stated above, in case of compositional data, and consequently also compositional tables, this needs to be done in a proper coordinate representation that maps the Aitchison geometry of compositions to a standard Euclidean geometry [28]. To eliminate the influence of outlying observations in PCA, in [15] it was proposed to estimate the covariance matrix for robust PCA by the Minimum Covariance Determinant (MCD) estimator [23]. Since centered logratio (clr) coefficients [1], that aggregate relative information on single compositional parts, lead to singularity and are not appropriate for most robust methods including the MCD estimator, loadings and scores of PCA need to be computed from isometric logratio (ilr) coordinates [9] of the compositional data and then transformed back to clr coefficients for a better interpretation of the resulting compositional biplot.

Accordingly, the aim of this paper is to generalize the previous considerations on dimension reduction of vector compositional data and to propose a robust approach to principal component analysis of compositional tables. In Sections 2 and 3, the logratio methodology of vector compositional data and their dimension reduction using principal component analysis, respectively, are briefly reviewed. Compositional tables and their interpretable processing using robust PCA are introduced in Section 4. The new methodology is illustrated in Sections 5 and 6 on real data sets from OECD Statistics using the statistical software R, namely the robCompositions package. Data from several different countries containing unemployment information with gender distribution and age structure are processed as a set of $2\times 4$ compositional tables. Therefore, a robust compositional biplot is a possible tool to analyze the distribution of unemployment rates in these countries as well as gender and age differences. Data from the area of education, carrying relative information about fields of study and the resulting degree in given countries, are approached as larger $3\times 8$ compositional tables, and results for men and women are compared.

2 Logratio methodology of compositional data

In many applications, the relative structure of the observations is more interesting than the absolute values of their components. For example, when considering numbers of students attending bachelor, master and doctoral studies at different universities, the ratios among these three groups might be more relevant for a statistical analysis than just the empirical values, which might not be comparable because of different total student numbers. In other words, the actual total number of students (sufficiently high so that the impact of a measurement error with small sample sizes can be neglected) might be considered as not informative for the purpose of the analysis. However, to work with quantitatively described contributions on a given whole in a concise and meaningful manner, some concepts need to be introduced first.

A positive (row) vector $\boldsymbol{x}=(x_{1},x_{2},...,x_{D})$ is defined to be a D-part composition if it carries relative information, i.e. the ratios between the components are informative [1, 30]. Any compositional vectors with equal number of parts are considered to be representatives of the same equivalence class if one vector is obtained from another by a positive scalar multiplication [30]. Accordingly, equivalence classes of compositional data are represented without loss of information in a D-part simplex,

[TABLE]

for any $\kappa>0$. The choice of $\kappa$ (being 1 for proportions and 100 for percentages) is irrelevant for the analysis and can also vary throughout the compositional data set. Formally, the closure operation

[TABLE]

can be applied to rescale the data to a given constant sum ($\kappa$) representation. Accordingly, the $D$-part simplex is a sample space of (representatives of equivalence classes of) compositions. The constant sum representation is useful, e.g. for a comparison of several compositions from a sample. As an interesting consequence, possibly deviating (outlying) observations from the main data cloud are characterized by aberrant ratios rather than by significantly high or low absolute values of components [16].

As mentioned above, relevant information in data with compositional nature is contained only in ratios. Consequently, results of statistical processing should not depend on the sum $\kappa$ of compositional parts and instead of Euclidean distances, relative differences are used to express distances between observations. This principle called scale invariance is the first of three basic compositional principles [30]. Often the original data contain some non-informative part(s) in the compositional vector that are not of interest. Hence, we do not expect any change of results concerning the respective subcomposition when these parts are removed from the data. Subcompositional coherence is a principle stating that results obtained from a d-part subcomposition, $d<D$, are not in contradiction with results obtained by an analysis of the original D-part composition. Finally, permutation invariance states that the results are independent from a chosen order of parts within the composition, an expectable assumption for any reasonable statistical processing.

Due to the above principles and the relative scale of compositional data, the Euclidean geometry needs to be replaced with the Aitchison geometry [9]. Operations of perturbation and power transformation, being analogous to a sum of two vectors and a multiplication of a vector by a scalar in the real Euclidean geometry, are defined as

[TABLE]

where $\boldsymbol{x}$ and $\boldsymbol{y}$ are D-part compositions, and $\alpha$ is a real constant. Accordingly, operations of perturbation and power transformation form a $(D-1)-$dimensional vector space $(\mathcal{S}^{D},\oplus,\odot)$ [30].

To obtain Euclidean vector space structure, the Aitchison inner product, norm and distance are defined for D-part compositions $\boldsymbol{x}$ and $\boldsymbol{y}$ as

[TABLE]

[TABLE]

respectively, where $\boldsymbol{x}\ominus\boldsymbol{y}=\boldsymbol{x}\oplus\left[(-1)\odot\boldsymbol{y}\right]$.

Given the introduced specifics of compositional data endowed with the Aitchison geometry, standard multivariate statistical methods cannot be applied directly on raw data. It is advisable to employ the working on coordinates principle; that is to represent the compositional data in real coordinates before starting with a statistical analysis [29, 24]. There are two types of isometric coordinate representations with respect to the Aitchison geometry. Centered logratio coefficients [1] are defined as

[TABLE]

where $g(\boldsymbol{x})$ stands for the geometrical mean of the whole composition. Each clr coefficient aggregates all pairwise logratios with a given compositional part, thus enabling for a simple and meaningful interpretation in terms of dominance of that part with respect to the other parts on average. Consequently, clr coefficients are useful for a graphical interpretation of compositional data including compositional biplots as a result of a dimension reduction using PCA [2]. On the other hand, clr coefficients sum up to zero which leads to a singular covariance matrix, being inappropriate for processing using common robust statistical methods [15, 16].

Fortunately, isometric logratio coordinates $\boldsymbol{z}\in\boldsymbol{R}^{D-1}$ representing orthonormal coordinates with respect to the Aitchison geometry can be derived as

[TABLE]

where $D$-part compositions $\boldsymbol{e}^{i}=\mathcal{C}(e_{1}^{i},e_{2}^{i},...,e_{D}^{i})$, $i=1,...,D-1$, form an orthonormal basis on the simplex.

Obviously, the interpretation of ilr coordinates might be more tricky than in the case of clr coefficients as there are infinitely many possibilities of constructing ilr coordinates depending on the choice of basis vectors $\boldsymbol{e}^{i},\ i=1,\dots,D-1$. Sequential binary partitioning (SBP) of compositional parts is one possibility for providing a meaningful choice of $\boldsymbol{e}^{i}$ for the practitioner which is corresponding to the prior knowledge about compositions and resulting in coordinates called balances [8]. Moreover, there is a linear transformation between ilr coordinates and clr coefficients, done through a $D\times(D-1)$ matrix $\boldsymbol{V}$ of clr representations of the ilr basis vectors,

[TABLE]

Recently, pivot coordinates [17, 19]

[TABLE]

were introduced as a special case of ilr coordinates. They are appropriate especially in situations where no prior knowledge about how to perform SBP is available, e.g. in [3, 4, 6, 21]. Here, $x_{i}^{(l)}$ refers to the i-th part of the re-ordered composition $(x_{l},x_{1},\ldots,$ $x_{l-1},x_{l+1},\ldots,x_{D})$ which can be rewritten as $(x_{1}^{(l)},x_{2}^{(l)},\ldots,x_{l}^{(l)},x_{l+1}^{(l)},\ldots,x_{D}^{(l)})$. In each of the $D$ coordinate systems, a permutation of compositional parts needs to be performed, so that the l-th part of $\boldsymbol{x}$ ($l=1,\dots,D$) stands at the first position. Accordingly, the first pivot coordinate in each system, $z_{1}^{(l)}$, then clearly explains all relative information about part $x_{l}$ and, additionally, it is proportional to the respective clr coefficient

[TABLE]

The expression (8) holds also for pivot coordinates which is particularly useful in the context of this paper.

3 Robust principal component analysis for compositional data

When dealing with large-scale data sets, dimension reduction is often of primary interest. PCA is one of the widely used methods for this purpose also in a compositional approach, converting possibly correlated original variables from the data at hand into a smaller set of linearly uncorrelated variables called principal components (PCs). Additionally, the first component accounts for the largest variance of the given data, the second one for a maximum of the remaining variance, etc., under the constraint of being orthogonal to all the previous principal components [20].

The covariance matrix $\boldsymbol{C}$ estimated from a real data matrix $\boldsymbol{X}$ can be spectrally decomposed into

[TABLE]

where $\boldsymbol{G}$ is a matrix of eigenvectors and $\boldsymbol{L}$ represents a diagonal matrix of eigenvalues of $\boldsymbol{C}$. It is then possible to define the PCA transformation as

[TABLE]

where $\boldsymbol{t}$ is the location estimator and $\boldsymbol{1}$ is a vector of ones with length $n$. The columns of the matrix $\boldsymbol{X}^{\ast}$, the coordinates of the principal components, are called scores and the columns of $\boldsymbol{G}$, containing the respective basis vectors, are called loadings. Typically, only the first few principal components are considered for further analysis. Taking into account only two PCs, a graphical outcome called biplot can depict both loadings as arrows and scores as points in one plot, where associations can be revealed.

It is common to take $\boldsymbol{t}$ as the arithmetic mean and $\boldsymbol{C}$ as the sample covariance matrix, however, both are very sensitive to outlying observations. Robust alternatives can be obtained by using the Minimum Covariance Determinant (MCD) estimators of location and covariance [23]. However, this approach inquires working in ilr coordinates to obtain full rank data in order to get the MCD estimate of the covariance matrix and the respective matrix of eigenvectors $\boldsymbol{G}$. In addition, ilr coordinates ensure subcompositional coherence and enable to keep affine equivariance of the results to the change of basis.

Accordingly, robust principal component analysis of compositional data based on the MCD estimator requires ilr coordinates $\boldsymbol{z}_{i}$ as an input, and the scores $\boldsymbol{z}_{i}^{\ast}$ are given by

[TABLE]

Once PCA is performed, the loadings can be transformed back to clr coefficients as

$\boldsymbol{G}_{\mathrm{clr}}=\boldsymbol{VG}$,

accounting for compositional biplot construction with meaningful interpretation, whereas the scores remain identical and only a column of zeros is added to the end. Clr coefficients are also worth as such for their simple construction as an amalgamation of pairwise logratios of a given part. Due to the zero-sum constraint of clr coefficients, their covariance structure is distorted, thus the interpretation of the biplot in terms of the correlation between coefficients (through angles between arrows) might be misleading. Instead, the focus is on links between vertices of arrows as they stand for a proportionality between the original compositional parts [2]. On the other hand, due to the link with pivot coordinates, the single clr variables (or the respective loadings) can be used to identify observations with a high dominance of the respective parts in a compositional vector [22].

4 Compositional tables

Compositional tables describe quantitatively relative contributions to a given whole that is distributed according to two (row and column) factors. Mathematically, this leads to a two-factorial extension of vector compositional data, carrying information about a relationship between these factors. A compositional table (1) thus can be represented, e.g. either as a contingency table (with sufficiently high numbers of counts in the cells) or as a table of the same order with ML estimates of the respective probabilities – due to scale invariance, the relative information (contained in the ratios between the cells) is the same in both cases [7, 10]. Hence, the concept of compositional tables (1) covers both the discrete case of contingency tables and its continuous counterpart (e.g. input-output tables, see [11]). Nevertheless, in the compositional context, a concrete table represents usually just one realization in a sample from a multivariate continuous distribution. Due to the decision to treat the data from such a distribution compositionally, the possible order of the factor categories (e.g. age levels) is ignored, making this a relevant subject for future research.

Since compositional tables form a direct extension of vector compositional data, all the principles and operations introduced in Section 2 apply, up to some minor modifications due to the two-factorial structure of the tables.

Accordingly, the closure operation

[TABLE]

is used to represent a compositional table $\boldsymbol{x}$ in an $IJ$-part simplex $\mathcal{S}^{IJ}$ of vectorized tables $\mathrm{vec}(\boldsymbol{x})=(x_{11},\ldots,x_{I1},\ldots,x_{IJ})$. Perturbation $\boldsymbol{x}\oplus\boldsymbol{y}$, powering $\alpha\odot\boldsymbol{x}$, and the Aitchison inner product $\left\langle\boldsymbol{x},\boldsymbol{y}\right\rangle_{A}$ of two tables $\boldsymbol{x},\,\boldsymbol{y}$ and a real number $\alpha$ can be defined analogously [7, 10]:

[TABLE]

[TABLE]

It is straightforward to derive that the dimension of the simplex $\mathcal{S}^{IJ}$ is $IJ-1$, corresponding to the dimensionality of $(I\times J)$-compositional tables.

Permutation invariance and subcompositional coherence are valid with respect to the two factors of the compositional tables, allowing to permute and discard entire rows or columns only.

To analyze compositional tables, it is beneficial to work also with the so-called independence and interaction tables. These can be obtained through an orthogonal decomposition [7]

[TABLE]

Here, the independence table is constructed to extract all the relative information about row and column factors under the assumption that the original compositional table is a product of its row and column geometric marginals, and the interaction table contains information about the relationships between the row and column factors. Therefore, in case of actual independence in the data at hand (in the above sense), all the entries of the interaction table are the same, since there is no remaining information left in the data after the decomposition; the interaction table thus forms a neutral element with respect to the Aitchison geometry of compositional tables. Otherwise, the interactive part describes the nature of the deviation from an independent situation. Needless to say, analyzing each of these new tables separately allows for a deeper insight into the original data.

It turns out that the introduced decomposition can be easily derived from row and column projections of the compositional table onto marginal subspaces (for further details, see [7]),

[TABLE]

[TABLE]

where $g(.)$ denotes the geometric mean of the cells in the argument and $\bot$ stands for orthogonality of the projections.

Recalling the case of independence in probability tables, it is instant to get the independence table simply by perturbing both these projections, $\boldsymbol{x}_{ind}=\mathrm{row}^{\bot}(\boldsymbol{x})\oplus\mathrm{col}^{\bot}(\boldsymbol{x})$. From (16) it follows that the interaction table is just a decomposition remainder in $\boldsymbol{x}_{int}=\boldsymbol{x}\ominus\boldsymbol{x}_{ind}$. For practical calculations, the following formulas are used to obtain the single entries of these tables,

[TABLE]

[TABLE]

It is crucial to realize that the dimensions of $\boldsymbol{x}_{ind}$ and $\boldsymbol{x}_{int}$ lower to $I+J-2$ for the independence tables, which follows immediately from the dimensions of the row and column projections being, respectively, $I-1$ and $J-1$, and to $(I-1)(J-1)$ for the interaction tables, which is easily obtained from the orthogonality of the decomposition.

As stated in the previous sections, a coordinate representation which respects the sample space dimensionality is needed to perform robust PCA of compositional data. In case of compositional tables (and particularly their decomposed parts), a generalization of balance coordinates needs to consider two SBPs according to each factor [13]. Even for moderate numbers of rows and columns, the interpretation of this coordinate representation gets rather complex without a deeper expert knowledge. Therefore, only a two-factorial alternative to pivot coordinates is usually used in practice [12, 13]. Interestingly, the coordinates of the entire compositional table can be divided into two groups according to the dimensionality of the independence and interaction tables, respectively. This becomes the main advantage also when using pivot coordinates for robust PCA since it allows for a comparison of the results from the whole table and its decomposed parts.

Generally, there are three types of pivot coordinates corresponding to the row, column and “odds ratio” partitioning of the compositional table [12]. The first two types together form a coordinate representation of the independence table, the third one is used for the interaction table. Altogether, they form a coordinate representation of the original compositional table. In case of row and column types of coordinates, the entire first row or column, respectively, is taken as the pivot element and separated from the rest. In the next step, this pivot is not considered anymore and the following row or column is taken as the new pivot element, and so on, until the following $I+J-2$ coordinates are obtained,

[TABLE]

[TABLE]

where $g(\boldsymbol{x}_{i\bullet})$ and $g(\boldsymbol{x}_{\bullet j})$ stand for the geometric mean of the $i$-th row and $j$-th column, respectively.

The process of obtaining the remaining $(I-1)(J-1)$ coordinates is based on a division of the original compositional table into four blocks, say upper left A, upper right B, lower left C and lower right D, where A contains always just one (pivot) cell indexed by $rs$. The odds ratio interpretation should be now easily seen from the following formula, where the elements of blocks A and D are in the numerator, and the elements of blocks B a C in the denominator of the logratio;

[TABLE]

To obtain all coordinates of the odds ratio type in a proper order corresponding to the $z^{r}$ and $z^{c}$ coordinates, the position of the pivot cell is moving firstly by rows with fixed first column, $r=1,\ldots,I-1$, then by columns with fixed last row, $s=1,\ldots,J-1$, and afterward the row position is always leveled back down by one and the column position moves again from $1$ to $J-1$ for the given row until all sizes of the $r\times s$ table are covered.

For the sake of completeness, permutations of the entire rows or columns following the same principle as stated in Section 2 could be performed. Hereby for all combinations of rows and columns, different coordinate systems consisting of $z_{i}^{r(k)}$, $z_{j}^{c(l)}$ and $z_{rs}^{OR(kl)}$, where $(kl),k=1,\ldots,I,l=1,\ldots,J,$ defines row and column permuted to the pivot position within the whole table, would be gained [12]. Following (10), also the first coordinates of the three types from each system can then be expressed as proportional (up to a constant) to respective clr coefficients,

[TABLE]

[TABLE]

which is an important fact for the interpretation of the analysis. The resulting clr coefficients, computed originally from the elements of the independence and interaction tables,

[TABLE]

can thus be expressed also in terms of cells of the input compositional table as

[TABLE]

where $g(\boldsymbol{x}_{i\bullet})$, $g(\boldsymbol{x}_{\bullet j})$ and $g(\boldsymbol{x}_{\bullet\bullet})$ stand for the geometric mean of the $i$-th row, the $j$-th column and the whole compositional table (and its independent and interactive counterparts), respectively. As a consequence, each $\mathrm{clr}(\boldsymbol{x}_{ind})_{ij}$ expresses a dominance of a given combination of factor values in case of independence. This dominance is then either amplified or weakened according to the interaction table which depends on whether the interaction is shifted in a positive or a negative direction. The interaction table refers also to sources of departures from independence, nevertheless, the information obtained only from $\mathrm{clr}(\boldsymbol{x}_{int})_{ij}$ does not provide a complete picture about the dominance of the respective cell to all other averaged cells.

Furthermore, note that each coordinate $\mathrm{clr}(\boldsymbol{x}_{ind})_{ij}$ is formed by the sum of clr coefficients of the respective row and column marginals, $1/J\sum_{j}{\mathrm{clr}(\boldsymbol{x}_{ind})_{ij}}\\ =\ln g_{i\bullet}/g_{\bullet\bullet}$ and $1/I\sum_{i}{\mathrm{clr}(\boldsymbol{x}_{ind})_{ij}}=\ln g_{\bullet j}/g_{\bullet\bullet}$, which amount to zero. Thus, there are only $I+J-2$ linearly independent clr coefficients, reflecting the dimensionality of the sample space of independence tables again. A similar feature holds also for clr coefficients of interaction tables that sum up to zero across each row or column. Consequently, in the case of an interaction table, the number of linearly independent clr coefficients reduces to $(I-1)(J-1)$. Since this dependency makes it impossible to use the clr coefficients for the robust PCA of independence and interaction tables, the strategy to perform robust PCA for compositional tables is the same as in case of vector compositional data: PCA loadings and scores are computed in ilr coordinates and then back-transformed using relation (8) to the clr space, where the loadings can be interpreted in terms of dominance of single cells. Here, clr coefficients of basis vectors for rows $\boldsymbol{e}^{r}$, columns $\boldsymbol{e}^{c}$ and interactions $\boldsymbol{e}^{OR}$, forming the columns of the matrix $\boldsymbol{V}$, are defined as follows,

[TABLE]

[TABLE]

and

[TABLE]

reinterpreting the expressions from [12]. As a result of (8), row-wise clr coefficients of the whole table are obtained for the $IJ-1$ columns of the matrix $\boldsymbol{V}$. Alternatively, if the matrix $\boldsymbol{V}$ has just $I+J-2$ columns formed by clr coefficients of basis vectors corresponding to the ilr representation of the independence table (20), (21), its respective clr coefficients are derived (and similarly for the interaction table with its coordinates (22)). Finally, the transformed loadings and scores can be used to construct a biplot in order to reveal the multivariate structure of the sample of compositional tables and relations between both factors.

5 Unemployment data analysis

In the following, the methodological results are applied to two real-world data sets in order to illustrate the main features and possible limitations of the approach. In the first case, data from OECD Statistics about more than 150 million unemployed people from 42 different countries in 2010 [25] are analyzed using the statistical software environment R [31].

The data set contains the numbers of unemployed people together with their gender and age category for the following countries: Australia, Austria, Belgium, Canada, Chile, Czech Republic, Denmark, Estonia, Finland, France, Germany, Greece, Hungary, Iceland, Ireland, Israel, Italy, Japan, Korea, Mexico, Netherlands, New Zealand, Norway, Poland, Portugal, Slovakia, Slovenia, Spain, Sweden, Switzerland, Turkey, United Kingdom, United States, Colombia, Costa Rica, Latvia, Lithuania, China, India, Indonesia, Russian Federation and South Africa. An example of (transposed) raw data from the first four countries is shown in Table 1. The numbers in the tables are basically counts of unemployed people according to two factors. As the population size varies among the countries, the interest here is not in the absolute values of the counts in the single countries, but rather the relative structure of unemployment. Particularly, ratios of men and women and ratios among age groups 15-24, 25-39, 40-54 and 55+, as well as proportionality among countries will be analyzed. Since outliers can be anticipated due to completely different economics, education levels, gender balance and also traditions of the listed countries, the analysis will be carried out in a robust manner.

All compositional tables in this example have two rows and four columns, i.e., gender is the row factor and age structure is the column factor. The sample space of tables thus has dimension $7$ out of which independence tables account for a dimension of $4$ with pivot coordinates $z_{1}^{r}$, $z_{1}^{c}$, $z_{2}^{c}$ and $z_{3}^{c}$, while the remaining coordinates $z_{11}^{OR}$, $z_{12}^{OR}$ and $z_{13}^{OR}$ correspond to the interaction tables with a dimension of $3$.

To point out the differences between the classical and robust PCA, both are performed and compared through the resulting covariance compositional biplots. Recall that classical PCA can directly be applied on clr coefficients. Nevertheless, since for this data set robust PCA may be more relevant because of potential outlying tables, ilr coordinates are used in both cases, and the results are transformed to clr for the biplot construction. This can only result in a different rotation of the classical biplot, however, it obviously does not alter the results.

In order to perform PCA in ilr coordinates, the standard function princomp in R can be used, where the parameter covmat is set to covMcd (MCD estimator of covariance) in case of robust PCA. Thereafter, loadings need to be transformed to clr coefficients as described in Section 2 using the matrix $\boldsymbol{V}$ with columns defined by (27) - (29) for the entire compositional table, and by (27) and (28), or by (29) for its independent and interactive part, respectively. The resulting classical biplots are depicted on the right-hand side of Figure 1, while the robust PCA output is on the left.

Assessing Figure 1, it can be noticed that in all three cases, robust PCA performs better in terms of explained variability by the first two principal components. As mentioned above, some outliers might be present in the data, and even from the classical biplot of the whole compositional tables (upper right corner of Figure 1), at least two outlying tables (Turkey-TUR and South Africa-ZAF) could be expected, and thus the robust approach should provide more meaningful results. However, an outlier detection is performed additionally in order to confirm these expectations. In the R package robCompositions [32], there is a function outCoDa defined for this purpose, based on robust Mahalanobis distances computed from ilr transformed data [14]. Using the pivot coordinates and applying the $0.975$ quantile of the chi-squared distribution as the common cut-off value, 15 out of all 42 countries are identified as outlying observations, clearly supporting the choice of robust analysis. Note that similarly, 10 observations from the set of independence tables and 6 from the interaction tables were detected as potential outliers.

Additionally, from the same part of the figure, it is easy to identify from the direction of the arrows which countries tend to have relatively higher unemployment among younger people and which ones have a rather higher rate in the opposite situation. Although no clear compact clusters are visible, it is obvious that most European countries together with the USA and Canada tend to have most likely problems with employing older people, say 40+, while for Central and South America together with China, India, and Indonesia the unemployment depending on age structure has rather opposite tendencies. Moreover, also some gender differences can be observed, except for the youngest generation. The structure in the classical biplot (upper right plot) is similar but driven by the identified outlying observations.

The left plot in the middle part of Figure 1 shows the “ideal” situation in case the relationships between gender and age factors would be filtered out. While the positions of the countries are not apparently changed compared to the previously discussed covariance biplot (upper right corner), the general relationships between the factors are remarkably illustrative. In case of independence, nearly gender equity would be achieved, while on the contrary, relationships among the age levels would be disproportionally weaker. Also, a bigger difference between results provided by robust and classical PCA is present here. One can easily understand how the classical approach does not handle outliers and how those can affect the output; the biplot on the right side is quite far away from picturing the same ideal situation.

As demonstrated, the independence table captures the hypothetical balanced state with each clr interpreted in terms of dominance of a given combination of factors in case of independence. However, this dominance is then either amplified, or weakened according to the interaction table; in terms of clr coefficients, it depends on whether the logratio dominance is shifted in a positive, or in a negative direction. Note that information obtained from the interaction table only does not provide a complete picture about the dominance of single cells in the table. For example, in the lower left biplot (robust biplot of interaction tables), Costa Rica (CRI) is placed towards the loading “male 55+”, but this does not necessarily lead to a conclusion that unemployment in this group is higher in general in this country; it simply marks the cell whose dominance causes imbalance for Costa Rica, although the actual proportion of unemployment for this age group can be lower than its average dominance. Therefore, the conclusion about the higher dominance of unemployed men in the oldest group than expected in the hypothetical case of independence can be stated only after looking at the biplot of the independence tables. For the compositional tables with dimension $2\times J$ (or alternatively $I\times 2$), this feature is nicely illustrated by the depicted loadings themselves, placed along a line corresponding to increasing dominance of one factor value at the expense of the latter value of the same factor. The opposite relation between the respective clr coefficients is clearly visible from both the biplots and the form of $\mathrm{clr}(\boldsymbol{x}_{int})_{ij}$ in (26): as it was already discussed in Section 4, clr coefficients of the interaction table sum up to zero across each row and column which results in the identity $\mathrm{clr}(\boldsymbol{x}_{int})_{1j}=-\mathrm{clr}(\boldsymbol{x}_{int})_{2j}$, holding for each $j$ when $I=2$ (and similarly for $J=2$). While in case of higher data dimension the property is no longer visible in the graphs, in this example it can be seen that the two possible values of the gender factor lead to precisely contradictory loadings for any chosen value of the age factor. Thus they might only carry the information about the origin of the dominance shift, but no longer about the direction of the shift for which the difference from the independence table has to be consulted.

6 Education data analysis

It was illustrated in the previous example how outliers can affect results of classical PCA. Especially the gender equity achieved in the robust biplot of the independence tables would not be present in the classical one. Hence, in this second example, only robust analysis outputs are discussed. The data set contains information about more than 7 million female and nearly 6 million male students, divided according to 8 different fields of study, being Education, Humanities and arts, Social sciences, business and law, Science, mathematics and computing, Engineering, manufacturing and construction, Agriculture and veterinary, Health and welfare, and Services [26]. The information about the achieved degree (bachelor, master, and doctoral, respectively) is recorded as well for about 30 different countries.

Compositional tables are analyzed for both genders separately in order to allow for a comparison of possible differences between them later on. Biplots as graphical PCA outcomes of the whole compositional table as well as independence and interaction tables are collected in Figure 2. Due to a larger dimensionality of compositional tables than in previous case ($3\times 8$), the biplots contain three times more variables and an objective interpretation becomes more difficult. An additional aspect is that since it is necessary to go many dimensions down to achieve the PCA projection using the first two principal components, it is expected to obtain more approximative picture of the multivariate data structure in the biplot. However, for data of similar or even bigger size, the proposed methods still offer an extremely useful rank-two approximation capturing the relationships between both factors. The performance of robust PCA is still good enough (at least $54.08\%$) for this concrete case in terms of cumulative variability explained by the first two principal components.

Despite the previous interpretational doubts, it can be seen that the effect of the chosen final degree is possibly stronger than the effect of the study field since the loadings tend to create a quite clear division of bachelors, masters, and doctors for most of the biplots. This property is more obvious for men while for women the difference between bachelors and masters is partially wiped out, maybe also due to a less compact data structure. Finally, employing the outCoDa function [32, 14] again, some outliers for the data set of the whole compositional tables are detected: Austria, Norway and Spain for men, and United Kingdom, Turkey and United States for women.

From the second row of the figures, some overall idea about the hypothetical state of independence between degree and study field factors might be obtained. A stronger effect of the chosen degree and a weaker effect of the study field would still be apparent for men, and one new feature could be observed: there would be a strong relation between Education and Humanities and arts study fields for each degree. For women, e.g. a similarity of educational systems in Sweden, Finland, Belgium and Switzerland is reflected. Also, in the case of independence, higher occurrence of outliers would be present for both men and women.

When looking at the biplot for the interaction tables for women (lower right figure), two of the mentioned countries are shifted away from the fields of study that would dominate if independence was achieved, being especially Agriculture and veterinary, Science, mathematics and computing, and Services. These countries are Sweden and Belgium, while Finland would correspond approximately to the independence between the factors, and Switzerland actually accounts for even stronger dominance of those fields (particularly master and doctoral studies in Services). For both men and women, it could be stated that the actual relationships between the factors are quite distant from the relative dominance given in the independence state. Stronger patterns concerning both factors are generally observed for men.

7 Conclusions

Vector compositional data, or their generalization to compositional tables, are characterized by their relative scale and scale invariance properties captured by the Aitchison geometry. Accordingly, a representation of the compositional tables isometrically in clr coefficients or ilr coordinates is essential for the further statistical analysis using popular multivariate methods. Compositional tables can be decomposed onto their independence and interaction parts; a statistical analysis of both is recommended to get insight into the ideal situation when relationships between factors are filtered away, as well as into interactions between factors forming the compositional table. As most real-world datasets contain outlying observations, robust methods requiring an orthonormal coordinate representation have been considered. To reduce the dimension of data at hand, a robust PCA using the MCD estimator can be applied to pivot coordinates of compositional, independence and interaction tables (but also any other ilr coordinates which respect the decomposition of compositional tables could be used). The necessity of respecting dimensionality of independence and interaction tables forms the main difference to (vector) compositional data where such feature does not occur. In case of pivot coordinates, it is also possible to find a simple relationship between them and the respective clr coefficients as well as to enhance interpretation of the latter variables. Therefore, PCA loadings obtained in pivot coordinates are transformed back to clr coefficients where they are used for the construction of biplots and their meaningful analysis. In case of $(2\times J)$ table dimensions, an additional feature can be observed in the graphical output of interaction tables, which can be traced back to the interpretation of the clr coefficients. As a next step, three factors of a compositional cube $(I\times J\times K)$ could be analyzed by robust PCA and appropriate 3D graphical tools employed. Consequently, the educational datasets for men and women could be merged into one and approached accordingly. The case of compositional cubes and also its further extension to generally $N$ factors belong to primary research interests of the authors.

Funding

This work was supported by COST Action CRoNoS under Grant IC1408; The Czech Science Foundation under Grant 18-05432S; and the grant IGA_PrF_2018_024 Mathematical Models of the Internal Grant Agency of the Palacký University in Olomouc.