Partial Correlations in Compositional Data Analysis

TL;DR

This paper proposes a natural and advantageous definition of partial correlations for compositional data, addressing issues with existing measures and enabling more accurate association analysis.

Contribution

It introduces a new partial correlation measure for compositional data based on residuals of log-ratios, improving interpretability and addressing limitations of current methods.

Findings

Partial correlations are well-defined and scale properly for compositional data.

Residuals of log-ratios are reference-independent, simplifying analysis.

The method allows for negative associations and is compatible with normalization approaches.

Abstract

Partial correlations quantify linear association between two variables adjusting for the influence of the remaining variables. They form the backbone for graphical models and are readily obtained from the inverse of the covariance matrix. For compositional data, the covariance structure is specified from log ratios of variables, so unless we try to "open" the data via a normalization, this implies changes in the definition and interpretation of partial correlations. In the present work, we elucidate how results derived by Aitchison (1986) lead to a natural definition of partial correlation that has a number of advantages over current measures of association. For this, we show that the residuals of log-ratios between a variable with a reference, when adjusting for all remaining variables including the reference, are reference-independent. Since the reference itself can be controlled for, correlations between residuals are defined for the variables directly without the necessity to recur to ratios except when specifying which variables are partialled out. Thus, perhaps surprisingly, partial correlations do not have the problems commonly found with measures of pairwise association on compositional data. They are well-defined between two variables, are properly scaled, and allow for negative association. By design, they are subcompositionally incoherent, but they share this property with conventional partial correlations (where results change when adjusting for the influence of fewer variables). We discuss the equivalence with normalization-based approaches whenever the normalizing variables are controlled for. We also discuss the partial variances and correlations we obtain from a previously studied data set of Roman glass cups.

Figures

Tables

Table 1. Table 1: Single-variable quantities for Roman glass cups. All values are given in percent.

oxide, element	av. weight	res. var/tot.	var/tot. (mean)	$R^2$ (mean)	var/tot. (SiO2)	$R^2$ (SiO2)
SiO2	72	0.86	7.7	89	-	-
Al2O3	1.9	0.37	2.5	85	2.0	90
Na2O	18	0.75	3.6	79	1.4	70
Fe2O3	0.30	1.5	5.9	75	11	93
MgO	0.46	3.4	6.4	46	8.3	78
CaO	5.6	1.2	1.9	39	3.5	81
TiO2	0.07	3.5	4.9	29	8.2	77
MnO	0.01	21	29	26	30	61
Sb	0.35	22	28	23	22	45
P2O5	0.05	4.3	5.3	19	8.0	70
K2O	0.49	3.6	4.2	15	5.9	67

Table 2. Table 2: Partial correlations and their q 𝑞 q -values of five pairs of variables

variable 1	variable 2	partial correlation	q-value
Fe2O3	Al2O3	0.73	0.03
Al2O3	SiO2	0.72	0.03
Fe2O3	SiO2	-0.66
Na2O	CaO	0.60	0.18
MgO	Fe2O3	0.43	0.46

Equations

$$r_{1,2|3,4}(X_1,\dots ,X_4)=\mathrm{corr}\left(\log \frac{X_1}{X_4},\log \frac{X_2}{X_4}|\log \frac{X_3}{X_4}\right)$$

$$\begin{gathered} r_{1,2|3,4}(X_1,\dots ,X_4)=\mathrm{corr}\left(\log \frac{X_1}{X_3},\log \frac{X_2}{X_3}|\log \frac{X_3}{X_4}\right) \\ =\mathrm{corr}\left(\log \frac{X_1}{\sqrt{X_3{X}_4}},\log \frac{X_2}{\sqrt{X_3{X}_4}}|\log \frac{X_3}{X_4}\right)=\mathrm{corr}\left(\log \frac{X_1}{X_3},\log \frac{X_2}{X_4}|\log \frac{X_3}{X_4}\right). \end{gathered}$$

$$Y=\mathbf{F}Z.$$

$$\mathbf{\Sigma }=\mathbf{F\Gamma }{\mathbf{F}}^T.$$

$${\mathbf{\Gamma }}^{-}={\mathbf{F}}^T{\mathbf{\Sigma }}^{-1}\mathbf{F}$$

$${\hat{Y}}_j(Y_{\mathcal{C}})=\mathrm{cov}(Y_j,Y_{\mathcal{C}})\mathrm{var}(Y_{\mathcal{C}}{)}^{-1}{Y}_{\mathcal{C}}.$$

$$\mathrm{var}(Y_j|Y_{\mathcal{C}})=\mathrm{var}\left(Y_j-{\hat{Y}}_j(Y_{\mathcal{C}})\right)$$

$$\mathrm{cov}(Y_i,Y_j|Y_{\mathcal{C}})=\mathrm{cov}\left(Y_i-{\hat{Y}}_i(Y_{\mathcal{C}}),Y_j-{\hat{Y}}_j(Y_{\mathcal{C}})\right).$$

$$\mathrm{var}(Y_j|Y_{\{1,\dots ,D-1\}\backslash j})$$

$$\mathrm{corr}(Y_i,Y_j|Y_{\{1,\dots ,D-1\}\backslash \{i,j\}})$$

$$Y_j-{\hat{Y}}_j\left(Y_{{\mathcal{C}}_j}\right)$$

$$Y_j-Y_k-(\widehat{Y_j-Y_k})=Y_j-{\hat{Y}}_j-(Y_k-{\hat{Y}}_k),$$

$$\begin{gathered} {\hat{Y}}_k\left(Y_{{\mathcal{C}}_j}\right)=\mathrm{cov}(Y_k,Y_{{\mathcal{C}}_j})\mathrm{var}(Y_{{\mathcal{C}}_j}{)}^{-1}{Y}_{{\mathcal{C}}_j}={\left(\sum _{i\in {\mathcal{C}}_j}{\sigma }_{ki}{c}_{il}^{(-1)}\right)}_{l\in {\mathcal{C}}_j}^T{Y}_{{\mathcal{C}}_j}=({\delta }_{kl}{)}_{l\in {\mathcal{C}}_j}^T{Y}_{{\mathcal{C}}_j} \\ =Y_k. \end{gathered}$$

$$Z_j-{\hat{Z}}_j\left(Z_{{\mathcal{D}}_j}\right),$$

$${\gamma }_{jj}^{(-)}={\sigma }_{jj}^{(-1)},j=1,\dots ,D-1.$$

$$\mathrm{var}\left(Z_j|Z_{{\mathcal{D}}_j}\right)=1/{\gamma }_{jj}^{(-)}.$$

$$\mathrm{corr}\left(Z_i,Z_j|Z_{{\mathcal{D}}_{ij}}\right)=\frac{-{\gamma }_{ij}^{(-)}}{\sqrt{{\gamma }_{ii}^{(-)}{\gamma }_{jj}^{(-)}}},$$

$${\sigma }_{j|{\mathcal{D}}_j}^2(X)$$

$$r_{ij|{\mathcal{D}}_{ij}}(X)$$

$$\mathrm{var}({\hat{Y}}_j)/\mathrm{var}(Y_j)$$

$$\mathrm{var}({\hat{Z}}_j)/\mathrm{var}(Z_j)$$

$${\mathbf{Q}}_P=\mathbf{FP}{\mathbf{F}}^T{\mathbf{H}}^{-1},$$

$$Y^P$$

$${\mathbf{\Sigma }}_P$$

$$\begin{gathered} \widehat{Y_j}(Y_{{\mathcal{C}}_j}^P)=\mathrm{cov}\left(Y_j,Y_{{\mathcal{C}}_j}^P\right)\mathrm{var}(Y_{{\mathcal{C}}_j}^P{)}^{-1}{Y}_{{\mathcal{C}}_j}^P=\mathrm{cov}\left(Y_j,{\mathbf{Q}}_P{Y}_{{\mathcal{C}}_j}\right){\mathbf{C}}_P^{-1}{\mathbf{Q}}_P{Y}_{{\mathcal{C}}_j} \\ =\mathrm{cov}\left(Y_j,Y_{{\mathcal{C}}_j}\right){\mathbf{Q}}_P^T{\mathbf{C}}_P^{-1}{\mathbf{Q}}_P{Y}_{{\mathcal{C}}_j}=\mathrm{cov}\left(Y_j,Y_{{\mathcal{C}}_j}\right){\mathbf{Q}}_P^T{\left({\mathbf{Q}}_P\mathbf{C}{\mathbf{Q}}_P^T\right)}^{-1}{\mathbf{Q}}_P{Y}_{{\mathcal{C}}_j} \\ =\mathrm{cov}\left(Y_j,Y_{{\mathcal{C}}_j}\right){\mathbf{Q}}_P^T{\left({\mathbf{Q}}_P^T\right)}^{-1}{\mathbf{C}}^{-1}{\mathbf{Q}}_P^{-1}{\mathbf{Q}}_P{Y}_{{\mathcal{C}}_j}=\mathrm{cov}\left(Y_j,Y_{{\mathcal{C}}_j}\right){\mathbf{C}}^{-1}{Y}_{{\mathcal{C}}_j}=\widehat{Y_j}(Y_{{\mathcal{C}}_j}). \end{gathered}$$

$${\hat{Y}}_j(Y)=\mathrm{cov}(Y_j,Y){\mathbf{\Sigma }}^{-1}Y=\mathrm{cov}(Y_j,\mathbf{F}Z){\mathbf{\Sigma }}^{-1}\mathbf{F}Z=\mathrm{cov}(Y_j,Z){\mathbf{F}}^T{\mathbf{\Sigma }}^{-1}\mathbf{F}Z=\mathrm{cov}(Y_j,Z){\mathbf{\Gamma }}^{-}Z={\hat{Y}}_j(Z).$$

$$Y_j-{\hat{Y}}_j=Z_j-Z_D-\widehat{(Z_j-Z_D)}=Z_j-{\hat{Z}}_j-Z_D+{\hat{Z}}_D,$$

$$\begin{gathered} {\hat{Z}}_D(Y)=\mathrm{cov}\left(\log \frac{X_D}{g(X)},\log \frac{X}{X_D}\right){\mathbf{\Sigma }}^{-1}Y=\mathrm{cov}\left(-\log \frac{g(X)}{X_D},\log \frac{X}{X_D}\right){\mathbf{\Sigma }}^{-1}Y \\ =\mathrm{cov}\left(-\frac{1}{D}\sum _{i=1}^{D-1}{Y}_i,Y\right){\mathbf{\Sigma }}^{-1}Y=-\frac{1}{D}{\left(\sum _{i=1}^{D-1}{\sigma }_{il}\right)}_{l=1,\dots ,D-1}^T{\mathbf{\Sigma }}^{-1}Y=-\frac{1}{D}{\left(\sum _{l=1}^{D-1}\sum _{i=1}^{D-1}{\sigma }_{il}{\sigma }_{lk}^{(-1)}\right)}_{k=1,\dots ,D-1}^{T}Y \\ =-\frac{1}{D}{\left(\sum _{i=1}^{D-1}{\delta }_{ik}\right)}_{k=1,\dots ,D-1}^{T}Y=-\frac{1}{D}\sum _{k=1}^{D-1}{Y}_k=-\frac{1}{D}\sum _{k=1}^{D-1}\log \frac{X_k}{X_D}=\log \frac{X_D}{g(X)}=Z_D. \end{gathered}$$

$$\begin{gathered} \mathbf{P}{\mathbf{\Gamma }}^{-}{\mathbf{P}}^T=\mathbf{P}{\mathbf{F}}^T{\mathbf{\Sigma }}^{-1}\mathbf{F}{\mathbf{P}}^T=\mathbf{P}{\mathbf{F}}^T\left({\mathbf{Q}}_P^T({\mathbf{Q}}_P^T{)}^{-1}\right){\mathbf{\Sigma }}^{-1}\left({\mathbf{Q}}_P^{-1}{\mathbf{Q}}_P\right)\mathbf{F}{\mathbf{P}}^T \\ =\mathbf{P}{\mathbf{F}}^T{\mathbf{Q}}_P^T{\left({\mathbf{Q}}_P\mathbf{\Sigma }{\mathbf{Q}}_P^T\right)}^{-1}{\mathbf{Q}}_P\mathbf{F}{\mathbf{P}}^T=\mathbf{P}{\mathbf{F}}^T{\mathbf{Q}}_P^T{\mathbf{\Sigma }}_P^{-1}{\mathbf{Q}}_P\mathbf{F}{\mathbf{P}}^T, \end{gathered}$$

$${\mathbf{F}}^T{\mathbf{Q}}_P^T={\mathbf{F}}^T{\left(\mathbf{FP}{\mathbf{F}}^T{\mathbf{H}}^{-1}\right)}^T={\mathbf{F}}^T({\mathbf{H}}^{-1}{)}^T{\left(\mathbf{FP}{\mathbf{F}}^T\right)}^T={\mathbf{F}}^T{\mathbf{H}}^{-1}\mathbf{F}{\mathbf{P}}^T{\mathbf{F}}^T,$$

$$\mathbf{G}$$

$$\mathbf{PG}{\mathbf{P}}^T$$

$$\mathbf{FG}$$

Full text

Partial Correlations in Compositional Data Analysis

**Ionas Erb1

1**Centre for Genomic Regulation (CRG),

The Barcelona Institute of Science and Technology,

Barcelona, Spain; ionas.erb@crg.eu

Summary

Partial correlations quantify linear association between two variables adjusting for the influence of the remaining variables. They form the backbone for graphical models and are readily obtained from the inverse of the covariance matrix. For compositional data, the covariance structure is specified from log ratios of variables, so unless we try to “open” the data via a normalization, this implies changes in the definition and interpretation of partial correlations. In the present work, we elucidate how results derived by Aitchison (1986) lead to a natural definition of partial correlation that has a number of advantages over current measures of association. For this, we show that the residuals of log-ratios between a variable with a reference, when adjusting for all remaining variables including the reference, are reference-independent. Since the reference itself can be controlled for, correlations between residuals are defined for the variables directly without the necessity to recur to ratios except when specifying which variables are partialled out. Thus, perhaps surprisingly, partial correlations do not have the problems commonly found with measures of pairwise association on compositional data. They are well-defined between two variables, are properly scaled, and allow for negative association. By design, they are subcompositionally incoherent, but they share this property with conventional partial correlations (where results change when adjusting for the influence of fewer variables). We discuss the equivalence with normalization-based approaches whenever the normalizing variables are controlled for. We also discuss the partial variances and correlations we obtain from a previously studied data set of Roman glass cups.

**Key words: ** Compositional covariance structure, inverse log-ratio covariance, residual log-ratio variance, compositional pairwise association, partial proportionality.

1 Introduction

1.1 Background and outline

Since the publication of Aitchison’s book on compositional analysis (Aitchison, 1986), awareness has increased about the fact that correlations between variables are problematic when these variables are parts of compositional data. To remedy the problems, a number of alternatives to correlation have been suggested. The most prominent is log-ratio variance (Aitchison, 1986), but there are also bounded (Aitchison, 2003) or scaled (Lovell et al., 2015) versions of it. All of these measures have their own problems: one would like to judge the value of log-ratio variance according to the intrinsic variability of its constituents, not only bound it between zero and one. Scalings however can lead to spurious results that become apparent when the underlying absolute data are known (Erb and Notredame, 2016). Apart from this, negative associations remain beyond the scope of all these measures of proportionality. Interestingly, the correlations between log-ratios that use the geometric mean as a reference, i.e., in clr-transformed data, have found little acceptance although they are perhaps the solution coming closest to a genuine correlation for compositional data. Its main drawback is subcompositional incoherence. However, the insistence on subcompositional coherence, and perhaps a certain reluctance to the use of correlation in the community, seem to have blocked the way towards embracing partial correlations as a valid way of analyzing pairwise association between variables summing to a constant. Clearly, here the coherence of results obtained on subsets of variables is neither required nor possible: all partial correlations change when removing variables we control for. To the best of our knowledge, this fact has not been exploited in compositional analysis, and work on partial correlation is lacking. This is especially surprising as the necessary mathematical results were derived by Aitchison in the 1980s already.

Indeed, when considering pairwise association between variables while adjusting for dependence on all other variables including the reference, the correlation that remains seems the best we can expect from a linear pairwise association measure of compositional parts. In this paper, we use some results from Aitchison (1986) to derive that correlating log-ratios when their reference variables are controlled for makes the correlations reference-independent in the sense that any reference from the variables we control for can be used interchangeably without altering the result. For a simple example, take the case of four compositional random variables $X_{1},\dots,X_{4}$. Let us denote by

[TABLE]

the partial correlation between the logarithms of the first two variables referenced by $X_{4}$ adjusting for the influence of the log-ratio $X_{3}$ with the reference. The notation on the left-hand side seems ambiguous but it turns out that

[TABLE]

Although partial correlations of the form (1) are standard (in the sense that the theory of partial correlations applies directly to log ratios), it takes some results from the statistics of compositional data analysis to show that there is no dependence on the reference variable as long as it is included in the variables we control for. Thus partial correlations of log ratios are simpler than they appear at first sight, with some interesting implications. First, the number of correlations to interrogate is restricted to pairs and is thus drastically reduced with respect to the case where references have to be taken into account. Second, attempts to “open” the data, i.e. to analyze them in absolute terms after normalization with an unchanged reference appear of little use if for most cases the same results can be derived without normalization. Third, the covariance matrix of the geometric-mean-referenced parts (i.e., of the clr-transformed data) can be used via its pseudoinverse to obtain all necessary results.

1.2 Prerequisites

1.2.1 Data matrices as instances of compositional random vectors

We start with a hypothetical, real-valued but positive $N\times D$ data matrix with elements $a_{ij}$, where samples are indexed in the rows and variables in the columns. The $N$ rows of this matrix can be considered instances of a $D$-dimensional random vector $A=(A_{1},\dots,A_{D})$. In order to avoid problems with zeros when dealing with log ratios, we also assume $A_{j}>0$ throughout. Let us now consider the closure of these instances, i.e. the data matrix $(x_{ij})$ resulting when we divide each element $a_{ij}$ by its row sum. The corresponding random vector $X=A/\sum(A_{j})$ we call a compositional random vector. Compositional analysis concerns data where the total sum over a sample has no relevance to the analyst. This sum can vary between samples. The closure operation is just a convenient way of incorporating the fact that the total sum over the constituent random variables is arbitrary. The “absolute” data $a_{ij}$ are considered unavailable here and were only mentioned to clarify that they can underlie $x_{ij}$. Data other than $a_{ij}$ can of course lead to the same $x_{ij}$. Put another way, we use the random vector $X$ as our representative for the equivalence class of random vectors with the same relative relationships between their random variables.

Let us now introduce yet another kind of random vector, $Z$, which we call the log-ratio transformed random vector. It is defined by $Z=\log(X/g(X))$, where $g$ denotes the geometric mean over the variables in the random vector. We can restrict the variables to an index set $\mathcal{A}$ and write this projection as $X_{\mathcal{A}}$. The resulting denominator $g(X_{\mathcal{A}})$ is called the reference. In the case of $\mathcal{A}$ containing all indices of the composition under consideration, the transformation is called the centered log-ratio transformation (clr), and the resulting random vectors we denote by $Z$, or $Z_{\mathcal{A}}$, whatever subset is considered. Throughout the manuscript, we let the set indexing $Z$ also indicate the variables over which the geometric mean is taken (unless a single variable is indicated, like for $Z_{j}$). In the case of the reference containing a single variable, usually $X_{D}$, the log-ratio transformation is called additive (alr). Note that the use of the index $D$ is a matter of convenience, and a variable of interest can just be moved to this last position. In general, results will depend on this choice. We will denote these alr-transformed random vectors by $Y$, and their restrictions to a subset $\mathcal{A}$ as $Y_{\mathcal{A}}$.

The two types of log-ratio transformed random vectors can be transformed into each other using a matrix designed for this purpose (Aitchison, 1986). Let us denote by $\mathbf{I}_{d}$ the $d$-dimensional identity matrix and by $\mathbf{j}_{d}$ the $d$-dimensional vector with entries 1. Let us define the matrix $\mathbf{F}=[\mathbf{I}_{D-1},-\mathbf{j}_{D-1}]$, i.e. $\mathbf{F}$ is the $(D-1)\times D$ matrix resulting from writing $\mathbf{I}_{D-1}$ and $-\mathbf{j}_{D-1}$ side by side. It is then easy to verify that we have

[TABLE]

1.2.2 Compositional covariance specifications and inverse variance

Aitchison (1986) introduced three different matrices, each (equivalently) specifying the full covariance structure of a compositional data set. The first matrix has the log-ratio variances $\mathrm{var}(\log(X_{i}/X_{j}))$ as elements and will not concern us further. The other two specifications are covariance matrices of the log-ratio transformed random vectors. With the notation from the previous section, the $(D-1)\times(D-1)$ matrix of alr-transformed random vectors is $\mathbf{\Sigma}=(\sigma_{ij})=\left(\mathrm{cov}(Y_{i},Y_{j})\right)=\mathrm{var}(Y)$. Finally, the third matrix is the covariance of clr-transformed random vectors $\mathbf{\Gamma}=(\gamma_{ij})=\left(\mathrm{cov}(Z_{i},Z_{j})\right)=\mathrm{var}(Z)$. The change in log-ratio transformation results in a singular $D\times D$ matrix, with the singularity resulting from the constraint $\sum_{j}Z_{j}=0$. Similarly to (3), we can also transform the covariance matrices into each other:

[TABLE]

With partial correlations in mind, an important result derived by Aitchison (1986) concerns the relationship between the pseudoinverse $\mathbf{\Gamma^{-}}$ of $\mathbf{\Gamma}$ and the inverse variance $\mathbf{\Sigma}^{-1}$. We have

[TABLE]

(see Property 5.6 (a) in Aitchison, 1986). In this article, we will essentially elucidate the implications of this relationship.

1.2.3 Partial correlation: The standard setting

Due to its unconstrained nature, the covariance matrix $\mathbf{\Sigma}$ can be used to obtain partial correlations (between log ratios having $X_{D}$ as a reference) in the standard way. Let us quickly review this procedure (see, e.g., Whittaker, 1990). In the following, to achieve more economical expressions, we assume that all log-ratio transformed random variables are centered, i.e. their averages are zero (in case they are not, we just have to subtract their average from them). Let the random vector $Y_{\mathcal{C}}$ be composed of the set of random variables having indices in $\mathcal{C}$, i.e. $(Y_{i})_{i\in\mathcal{C}}$. The linear least squares predictor (LLSP) of $Y_{j}$ given $Y_{\mathcal{C}}$ is then defined by

[TABLE]

Here, $\mathrm{cov}(Y_{j},Y_{\mathcal{C}})$ is the row vector containing the covariances of $Y_{j}$ with the scalar random variables in $Y_{\mathcal{C}}$, and $\mathrm{var}(Y_{\mathcal{C}})^{-1}$ is the inverse covariance matrix of these random variables. Note that, as a consequence of this definition, $\hat{Y}_{j}$ lies in the space spanned by $Y_{\mathcal{C}}$ and the residual $Y_{j}-\hat{Y}_{j}(Y_{\mathcal{C}})$ is orthogonal to that space.111More precisely, the predictor and the residual generate $N$-sample vectors that lie in/are orthogonal to the space spanned by the $N$-sample vectors generated from the random variables in $Y_{\mathcal{C}}$. The residual variance

[TABLE]

is a useful summary statistic that tells us how well $Y_{j}$ can be predicted from the variables $Y_{\mathcal{C}}$. It is also known as partial variance. More generally, the partial covariance between two scalar random variables is defined as the covariance between their residuals:

[TABLE]

Partial correlations are obtained scaling by the corresponding residual variances in the standard way. There is, however, a shortcut to this. An important result states that partial variances and correlations adjusting for all the remaining variables can be obtained from the inverse covariance matrix. We have

[TABLE]

where $\sigma^{(-1)}_{ij}$ denote the elements of $\mathbf{\Sigma}^{-1}$.

2 Partial correlation on compositional data

2.1 The residual is independent of the choice of alr transformation

While $\mathbf{\Sigma}$ is a standard covariance matrix, in a compositional context its usefulness may not be apparent as it depends on the reference chosen. Rather surprisingly, it turns out that the inverse $\mathbf{\Sigma}^{-1}$ is essentially independent of the choice of reference. This is implied by (5), but although the equation may seem unsurprising given (4), the fact that a matrix containing all the parts can be obtained from one where a part is sacrificed as reference appears somewhat mysterious. This can be understood better when looking from another angle. Remember that inverse covariance matrices have diagonal elements related to the variance of the residuals, see (9). The fact that for $\mathbf{\Sigma}^{-1}$ the choice of reference has little importance can be made sense of if the residuals themselves are independent of this choice.

Remember that $Y_{j}$ is an alr-transformed random variable with reference $X_{D}$ and that we have $\mathbf{\Sigma}=\mathrm{var}(Y)$ with $Y$ having $D-1$ components. For removing $Y_{j}$ form $Y$, let us define the index set $\mathcal{C}_{j}=\{1,\dots,D-1\}\backslash j$. We now show that the residual

[TABLE]

remains the same irrespective of which reference is chosen from the set $\{X_{k}|k=1,\dots,D,k\neq j\}$ for constructing its constituent random variables $\log X_{i}/X_{k}$, $i\neq j,k$. For this, we first have to clarify that the space generated from the explanatory log ratios (the random vector $Y_{\mathcal{C}_{j}}$ we are adjusting for) remains unchanged when choosing another variable for reference. Intuitively, a change of reference in the variables we adjust for corresponds to an overall subtraction of one of these variables (the one that has the new reference in the numerator) and neither reduces nor moves us outside the space generated by $Y_{\mathcal{C}_{j}}$. As prediction is equivalent with linear projection into this space, the LLSP will remain the same. A formal proof that permuting the variables in $Y_{\mathcal{C}_{j}}$ does not change the LLSP is presented in the Appendix.222Note that this does not imply that the LLSP is independent of the choice of reference, i.e., in general we have $\widehat{\log(X_{j}/X_{D})}\neq\widehat{\log(X_{j}/X_{k})}$. Now let us come back to the residual. A change from reference $X_{D}$ to reference $X_{k}$ is achieved by subtracting $Y_{k}$ from $Y_{j}$ ($k\neq j,D$):

[TABLE]

where the equality comes form the linearity of the predictor. With the definition of the LLSP (6) and denoting the elements of $\mathrm{var}(Y_{\mathcal{C}_{j}})^{-1}$ by $c^{(-1)}_{il}$ we have

[TABLE]

Here, $\delta$ denotes the Kronecker delta. It appears because the elements of $\mathbf{\Sigma}$ and of $\mathrm{var}(Y_{\mathcal{C}_{j}})$ coincide for indices in $\mathcal{C}_{j}$. Comparing with the right-hand side of (12), we see that the term in brackets vanishes, thus establishing the independence of the residual from a reference $X_{k}$, $k\neq j$. Figure 1 shows an intuitive rendering of the argument we provided. For the argument it is crucial that the reference is coming from the constituent variables of the ratios we adjust for.

2.2 Connection with clr transformation

We have seen that our residuals do not depend on the particular alr transformation chosen (as long as the reference occurs as a variable in the ratios we are adjusting for). What about the reference coming from a clr transformation? Note that for an equivalent formulation, we need to consider the clr on the subspace of variables we are adjusting for. A similar argument as in the previous section (see Appendix) shows that an equivalent expression of our residual (11) that is symmetric with respect to the reference variables is

[TABLE]

where the geometric mean reference of $Z_{j}$ goes over the indices in $\mathcal{D}_{j}=\{1,\dots,D\}\backslash j$. Now from (5) we know that the diagonal elements of the pseudoinverse of $\mathbf{\Gamma}$ are

[TABLE]

Due to the symmetry of $\mathbf{\Gamma}$ with respect to reference, this also has to hold after permuting labels when constructing $\mathbf{\Sigma}$ (see Appendix for a proof). Thus, with (9), we conclude that for all indices $j=1,\dots,D$, the residual variance can be calculated from

[TABLE]

Similarly, the off-diagonal elements of $\mathbf{\Gamma}^{-}$ coincide with the ones of $\mathbf{\Sigma}^{-1}$ when they have indices $i,j<D$. Again invoking invariance with respect to label permutations, we conclude that for all $i,j=1,\dots,D$, partial correlations can be calculated from

[TABLE]

where we used the notation $\mathcal{D}_{ij}=\{1,\dots,D\}\backslash\{i,j\}$. Note that here the $Z_{i}$, $Z_{j}$ are different from the ones in (16) because the geometric mean reference is taken over $X_{\mathcal{D}_{ij}}$, not $X_{\mathcal{D}_{j}}$.

2.3 A general expression, correlation matrices, $R^{2}$

To summarize, we have shown that we can use as a reference either variable occurring in the explanatory ratios or the geometric mean taken over them. For compositional data it is thus unambiguous to define partial variances and correlations directly for the parts $X_{i}$ of the composition via

[TABLE]

where we mean the variances and correlations evaluated on residuals of the form (11) or (14) and with any reference coming from the log ratios we are adjusting for.

Correlation matrices are easier to interpret than covariance matrices, and the inverses of correlation matrices are often used in practice. Although $\mathbf{\Sigma}^{-1}$ and $\mathbf{\Gamma}^{-}$ coincide where they can, their correlation counterparts do not. Indeed, when scaling $\mathbf{\Sigma}$ and $\mathbf{\Gamma}$ to have unit entries in the diagonal,333Transforming a covariance into a correlation matrix is achieved by multiplying from the left and the right a suitably scaled diagonal matrix. after inverting the scaled version of $\mathbf{\Sigma}$ and pseudo-inverting the scaled version of $\mathbf{\Gamma}$, we have to re-scale both to unit diagonal to make them coincide again. Their multiple correlation coefficients $R^{2}$, however, remain different. As they quantify the amount of variance of $Y_{j}$ or $Z_{j}$ explained by their respective predictors (which depend on the reference), they have to be reference-dependent:

[TABLE]

Although $Y_{j}$ and $Z_{j}$ have the same partial variance $\sigma^{2}_{j|\mathcal{D}_{j}}(X)=1/\gamma_{jj}^{(-)}=1/\sigma_{jj}^{(-1)}$, their zero-order variances are different. When an $R^{2}$ needs to be reported, (21) would be a likely candidate or otherwise one would have to have a good reason to choose a more specific reference $X_{D}$.

2.4 Comparison with normalization-based partial correlations

An approach that is sometimes pursued (e.g., in genomics) consists in trying to recover the underlying absolute data $a_{ij}$ mentioned in section 1.2.1. This is achieved by means of a specific normalization, i.e., by multiplication of the rows of $(x_{ij})$ with factors $s_{i}$ that are proportional to $\sum_{j}a_{ij}$. These factors are unknown in principle (they got lost by the closure operation). They can sometimes be inferred from assumptions concerning the knowledge of (on average) unchanged $a_{ij}$ across rows (see the supplement to Quinn et al., 2018). If the assumptions are fulfilled, the resulting data matrix is no longer relative and there is no need for analyzing ratios. Note however that we can write the normalized data as a special case of log-ratio transformed data $Z=\log(X/g(X_{\mathcal{U}}))$, where the set $\mathcal{U}$ contains only indices of variables for which we have $g(A_{\mathcal{U}})=$const. (We have just transformed to logarithms of the normalized data, a common procedure in genomics.) When comparing with the partial correlations for log-ratio transformed data with more general references (19), we see that they coincide (and can be obtained assumption-free, without normalization) as long as the log ratios we are adjusting for also contain the variables indexed by $\mathcal{U}$.

3 Application to Roman glass-cup data

We evaluate (18) and (19) on a data set of 11 oxides and elements composing the glass of 47 Roman cups excavated in Colchester (Cool and Price, 1994). We use the version of the data presented in (Baxter et al., 1990). The few previous analyses have focused on ordination. One of the conclusions of Greenacre’s recent analysis (Greenacre, 2018) is that the ratio of SiO2/CaO is the significant dimension of variation, followed by the ratios of SiO2/Sb and Na2O/Sb, which were found to be of lesser influence for the multivariate structure. In Table 1, beside the average weight percentages of the various oxides and elements, we show the partial variances and log-ratio variances (with respect to the geometric mean). For better readability, both of them are divided by total variance as obtained from the trace of $\mathbf{\Gamma}$. Then, $R^{2}$ with the mean reference is shown. In the last two columns we also present log-ratio variance and $R^{2}$ with respect to SiO2. The oxides occurring in high weight percentages, with the exception of SiO2, tend to have low variances. The relatively high variance of SiO2 is however well explained by the other variables (highest $R^{2}$ with geometric mean reference). This may explain the good properties as a reference reported by Greenacre (2018). The two variables with the highest residual variance (MnO and Sb) also contribute most variance to the data. The values of $R^{2}$ with respect to SiO2 show that Fe2O3 and Al2O3 are best predicted using this reference.

In table 2 we show the five strongest partial correlations as well as a plug-in estimate of their local false discovery rates obtained under permutations of the samples for each variable in the data matrix (see Appendix). Of note, permuting the residuals directly leads to a less severe test and much lower $q$-values, presumably because a natural dependence among the residuals gets destroyed.

While Greenacre’s analysis stresses the importance of SiO2 and the oxides with high weight percentages for the multivariate structure, partial correlations and $R^{2}$ seem to point us to the importance of Fe2O3 and Al2O3 for the understanding of mutual dependencies. Of note is the conspicuous negative correlation between SiO2 and Fe2O3, the only negative association of importance that we find in this data set. It would remain undetected with current association measures.

4 Discussion

A proposal for partial proportionality by Erb and Notredame (2016) only controls for one ratio and thus introduces an additional parameter. Already the reference can be considered a nuisance parameter and needs justification when discussing results. The natural definition of partial correlation elucidated in this paper treats the reference like just another variable that is controlled for. While we cannot quite get around discussing its choice for $R^{2}$, partial variances and partial correlations are defined reference-free and can be evaluated as efficiently as in absolute analysis. Application to graphical models is a natural next step. For applications in genomics, e.g., the equivalence of the compositional approach with the one using a normalization of the data is interesting. The equivalence holds as long as the normalizing variables are controlled for in the sense discussed in section 2.4. In genomics it will also be necessary to modify the analysis for the many-variables setting. Here, the empirical covariance matrix is a poor estimator that can be improved by regularization, also allowing for the necessary covariance inversion.

Acknowledgements

I thank Cedric Notredame and Thom Quinn for their support and encouragement.