Variable Selection in Compositional Data Analysis Using Pairwise Logratios

Abstract

In the approach to compositional data analysis originated by John Aitchison, a set of linearly independent logratios (i.e., ratios of compositional parts, logarithmically transformed) explains all the variability in a compositional data set. Such a set of ratios can be represented by an acyclic connected graph of all the parts, with edges one less than the number of parts. There are many such candidate sets of ratios, each of which explains 100% of the compositional logratio variance. A simple choice consists in using additive logratios, and it is demonstrated how to identify one set that can serve as a substitute for the original data set in the sense of best approximating the essential multivariate structure. When all pairwise ratios of parts are candidates for selection, a smaller set of ratios can be determined by automatic selection, but preferably assisted by expert knowledge, which explains as much variability as required to reveal the underlying structure of the data. Conventional univariate statistical summary measures as well as multivariate methods can be applied to these ratios. Such a selection of a small set of ratios also implies the choice of a subset of parts, that is, a subcomposition, which explains a maximum percentage of variance. This approach of ratio selection, designed to simplify the task of the practitioner, is illustrated on an archaeometric data set as well as three further data sets in an “Appendix”. Comparisons are also made with existing proposals for selecting variables in compositional data analysis.

Appendix

1.1 A.1 Three Additional Data Sets

Three more data sets are analyzed, to demonstrate the benefit of using ALRs as a substitute for the full compositional data set. Two of these compositional data sets are taken from Aitchison (2005) and the third one is considered by Greenacre (2016) in the context of CA. For each data set, the sets of ALRs are computed, using each part in turn as the reference in the denominator. The set of ALRs that lead to inter-case distances that best match the logratio distances, using the Procrustes correlation as the criterion, is identified.

• Data Set 1 (Aitchison 2005)

• Minerals compositions: 21 samples, 8 minerals

• qu: Quartz or: orthoclase al: albite an: anorthite

• en: Enstatite ma: magnetite il: ilmenite ap: apatite

The ALRs with respect to quartz (qu) give the best agreement—the Procrustes correlation (between full space configurations) is equal to 0.995. Figure 10 shows the two-dimensional LRA based on all 28 logratios alongside the PCA of the 7 ALRs, showing the almost identical configurations of sample points.

Fig. 10

Two-dimensional biplots of mineral compositions: a LRA biplot; b PCA biplot of optimal set of ALRs

• Data Set 2 (Aitchison 2005)

• Activity pattern of a statistician: 20 days, 6 activities

• te = Teaching; co = consultation; ad = administration;

• re = Research; ot = other wakeful activities; sl = sleep

The ALRs with respect to sleep (sl) give the best agreement—the Procrustes correlation (between full space configurations) is equal to 0.960. Figure 11 shows the two-dimensional LRA based on all 15 logratios alongside the PCA of the 5 ALRs, showing the highly similar configurations of sample points. The first dimension of the ALR analysis accounts for a much higher percentage of variance, similar to the glass cup example in the main text, suggesting that there is only one relevant dimension and that the LRA analysis is inflated with redundant variance.

Fig. 11

Two-dimensional biplots of activity patterns of statisticians: a LRA biplot; b PCA biplot of optimal set of ALRs

• Data set 3 (see Greenacre 2016, Appendix E)

• Fatty acid data: 42 samples, 25 fatty acids with nonzero values

This data set consists of groups of marine organisms collected in three different seasons. The ALRs with respect to fatty acid 16:0 give the best agreement to the multivariate structure—the Procrustes correlation (between full space configurations) is equal to 0.989. Figure 12 shows the two-dimensional LRA based on all 300 logratios alongside the PCA of the 24 ALRs, showing the similar groupings of the three seasonal subsets of data, separated by the ALR analysis just as well as by the LRA. The four ratios that stand out in the contribution biplot on the right are made up of the four parts prominently radiating out from the centre in the LRA on the left, expressed relative to the more centrally located fatty acid 16:0 (Fig. 12).

Fig. 12

Two-dimensional biplots of fatty acids: a LRA biplot; b PCA biplot of optimal set of ALRs. In respective biplots, the labels of fatty acids and fatty acid ratios close to the center (i.e., with low contributions to the solution) have been omitted to improve legibility

Fig. 13

Graph of set of logratios in first column of Table 4

1.2 A.2 Procrustes Analysis and Procrustes Correlation

The following matrix formulation summarizes the computations required:

Suppose F₁ (n₁ × p) and F₂ (n₂ × p) are two matrices of coordinates defining two configurations of the same labelled points in separate p-dimensional spaces. Both matrices are column-centered (i.e., column means are zero). Then the following steps lead to the Procrustes correlation.

• 1. Normalize both matrices: ${\mathbf{F}}_{1}^{*} = {\mathbf{F}}_{1} /\sqrt {{\text{trace(}}{\mathbf{F}}_{1}^{\text{T}} {\mathbf{F}}_{1} )} , \, {\mathbf{F}}_{2}^{*} = {\mathbf{F}}_{2} /\sqrt {{\text{trace(}}{\mathbf{F}}_{2}^{\text{T}} {\mathbf{F}}_{2} )}$

• 2. Compute cross-product matrix: ${\mathbf{S}} = {\mathbf{F}}_{1}^{{* \, \text{T}}} {\mathbf{F}}_{2}^{*}$

• 3. Perform singular value decomposition (SVD): ${\mathbf{S}} = {\mathbf{UD}}_{\alpha } {\mathbf{V}}^{\text{T}}$

• 4. Procrustes rotation matrix: ${\mathbf{Q}} = {\mathbf{VU}}^{\text{T}}$

• 5. Sum of squared errors between normalized coordinates after rotation of the second matrix:

  $$ E = {\text{trace[(}}{\mathbf{F}}_{1}^{* \, } - {\mathbf{F}}_{2}^{*} {\mathbf{Q}})^{\text{T}} ({\mathbf{F}}_{1}^{* \, } - {\mathbf{F}}_{2}^{*} {\mathbf{Q}})] $$

• 6. Procrustes correlation: $r = \sqrt {1 - E}$

1.3 A.3 Comparison of the Present Logratio Approach with the Principal Balances of Martín-Fernández et al. (2018)

Martín-Fernández et al. (2018) developed an algorithm for a stepwise selection of ILR balances, by successively partitioning the parts using an exhaustive search at each step of this divisive algorithm. They apply their method to the ten-part Aar Massif geochemical data set from the book by Van den Boogart and Tolosana-Delgado (2013), and their approach uses unweighted parts, which is the present practice of the CODA school. A major difference between their approach and the one in the present article is they do not use variance explained in the sense used here, but rather “variance contained” in, or “variance contributed” to the logratio variance (although they sometimes do use the term “variance explained”, but they mean “variance contained”). This is a weaker criterion than the variance explained one that is proposed in the present study, because a part of variance contributed by a logratio or a balance is a measure in isolation from the remainder of the variability in the rest of the data set (see Sect. 3.5 of the article for more explanation). Thus, in order to compare our results with those of Martín-Fernández et al. (2018), the explained variances have had to be computed for the sequence of ILR balances published in that paper (Table 4, columns 3 and 4). In addition, the simpler approach of selecting logratios proposed in the present study was executed (Table 4, columns 1 and 2, see Fig. 13 for a graph of these ratios). As a yet further comparison, the simple logratios of amalgamated parts, using the same partitioning sequence as the ILR balances, were also computed and their explained variances computed—these can be termed “amalgamation balances” (Table 4, fifth column). Finally, the variances explained by the principal component axes (i.e., dimensions of the unweighted LRA of the data), which are the optimal explained variances, are reproduced (Table 4, columns 6 and 7). Note that these last explained variances are the only ones where the definition of variance explained is equivalent to variance contained.

Table 4 Cumulative explained variances of sequences of simple logratios, ILR balances, amalgamation balances and principal components

The results are also presented graphically in Fig. 14 in the style of Table 3 of Martín-Fernández et al. (2018). The results have been graphed in two separate figures for clarity. In both, the PCA sequence of cumulative explained variances is shown to give common reference points. In the left-hand figure, the first ILR, involving nine out of the ten parts, is higher by 1.5 percentage points compared to the first logratio Na₂O/MgO, involving only two parts. At steps 3, 4 and 5, the simple logratio sequence is superior to the ILR sequence, after which the two sequences converge. In the right-hand figure, the ILR balance sequence is superior to the amalgamation balance sequence for the first two steps, but afterwards, they are practically identical. Notice that the amalgamation balance sequence does not necessarily reach exactly 100% variance explained, but in this example, it reaches 99.97% variance explained using 9 balances, lacking only 0.03%.

Fig. 14

Plots of cumulative variances in Table 4. The optimal values, obtained by PCA, are shown in both plots as a reference for comparison

In conclusion, this is another example where the sequence of simple logratios seems perfectly adequate to explain the variance of the whole compositional data set. They are comparable to the ILR sequence in terms of explained variance, sometimes even outperforming it, and are much easier to compute and interpret. Using amalgamations instead of geometric means is an alternative way of defining balances, and these also have an easier interpretation in practice.

1.4 A.4 Simulation Study of the Ward Dendrogram as Parts are Sequentially Randomized

The idea of this simulation is to study how the dendrogram from the Ward clustering breaks down as parts are sequentially randomized (i.e., columns are randomly permuted) to simulate growing random noise in the data set. The values of each part (i.e., oxide in the Roman glass cups data set) are permuted in turn, the data reclosed and the Ward clustering repeated. The order of the parts randomized is from the part with the least part of variance to that of the highest part (the parts being randomized are shown in the boxes next to the dendrograms). Figure 15 is read in horizontal steps, and after three parts are randomized, the structure is still fairly stable, but starts to break down from the fourth part being randomized onwards. The element Si is kept fixed throughout, but by the last randomization, the whole data set has been effectively converted to noise.

Fig. 15

Sequence of weighted Ward clusterings of the glass cup data as elements shown in boxes are successively randomized in the compositional data matrix