Reference-Invariant Inverse Covariance Estimation with Application to Microbial Network Recovery

Let $\mathbf{p}=(p_{1},\ldots,p_{K+2})^{\prime}$ denote a vector of compositional probabilities satisfying that $p_{1}+\cdots+p_{K+2}=1$. The additive log-ratio (alr) transformation picks an entry of this vector as the reference and transforms the compositional vector using log ratios of each entry to the reference. Without loss of generality, suppose we pick the last entry as the reference, then the alr transformed vector becomes

$$\mathbf{z}=\left[\log\left(\frac{p_{1}}{p_{K+2}}\right),\ldots,\log\left(\frac{p_{K}}{p_{K+2}}\right),\log\left(\frac{p_{K+1}}{p_{K+2}}\right)\right]^{\prime}.$$

The transformed vector $\mathbf{z}$ is often assumed to follow a multivariate continuous distribution with a mean vector $\boldsymbol{\mu}$ and a covariance matrix $\boldsymbol{\Sigma}$. For example, $\mathbf{z}\sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$. Denote further the inverse covariance matrix $\boldsymbol{\Omega}=\boldsymbol{\Sigma}^{-1}$.

Similarly, if we pick another entry as the reference, we can define another alr-transformed vector. For simplicity of illustration, suppose we choose the second last entry to be the reference and consider the following alr transformation

$$\mathbf{z}_{p}=\left[\log\left(\frac{p_{1}}{p_{K+1}}\right),\ldots,\log\left(\frac{p_{K}}{p_{K+1}}\right),\log\left(\frac{p_{K+2}}{p_{K+1}}\right)\right]^{\prime},$$

where the subscript $p$ denotes the “permuted” version of $\mathbf{z}$. Similarly, define the mean vector of $\mathbf{z}_{p}$ by $\boldsymbol{\mu}_{p}$, the covariance matrix by $\boldsymbol{\Sigma}_{p}$, and the inverse covariance matrix by $\boldsymbol{\Omega}_{p}=\boldsymbol{\Sigma}_{p}^{-1}$.

A simple derivation implies that $\mathbf{z}_{p}$ is a linear transformation of $\mathbf{z}$ as $\mathbf{z}_{p}=\mathbf{Q}_{p}\mathbf{z}$, where

$$\mathbf{Q}_{p}=\left(\begin{matrix}\mathbf{I}_{K}&\mathbf{-1}\\ \mathbf{0}^{\prime}&-1\end{matrix}\right)$$

with $\mathbf{I}_{K}$ denoting the identity matrix, and $\mathbf{0}$ and $\mathbf{-1}$ denoting the column vectors with all $0$’s and all $-1$’s, respectively. It follows that $\boldsymbol{\mu}_{p}=\mathbf{Q}_{p}\boldsymbol{\mu}$, $\boldsymbol{\Sigma}_{p}=\mathbf{Q}_{p}\boldsymbol{\Sigma}\mathbf{Q}_{p}^{\prime}$, and $\boldsymbol{\Omega}_{p}=(\mathbf{Q}_{p}^{\prime})^{-1}\boldsymbol{\Omega}\mathbf{Q}_{p}^{-1}$. It is also worth noting that $\mathbf{Q}_{p}$ is an involutory matrix, i.e., $\mathbf{Q}_{p}^{-1}=\mathbf{Q}_{p}$.

The following theorem states the reference-invariance property of the inverse covariance matrix $\boldsymbol{\Omega}$ under the alr transformation.

Theorem 1 regards the reference-invariance property of the true value of the inverse covariance matrix $\boldsymbol{\Omega}$. It can also be extended to a class of estimators of $\boldsymbol{\Omega}$. Suppose we have i.i.d. observations of the compositional vectors $\mathbf{p}_{1},\ldots,\mathbf{p}_{n}$, and consequently, their alr transformed counterparts $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$. Then, we can construct an estimator of $\boldsymbol{\Sigma}$, denoted by $\hat{\boldsymbol{\Sigma}}$, based on the i.i.d. observations $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$. Furthermore, we can construct an estimator of $\boldsymbol{\Omega}$, denoted by $\hat{\boldsymbol{\Omega}}$, by taking its inverse or generalized inverse. The following corollary presents the reference-invariance property for a class of such estimators.

The above results imply an important property for the additive log-ratio transformation in the compositional data analysis. It can be extended to a more general situation as follows. In general, suppose we have selected a set of entries $\mathbf{p}_{\mathcal{R}}$ as “candidate references” in a compositional vector $\mathbf{p}$ and write $\mathbf{p}=(\mathbf{p}_{\mathcal{R}^{c}}^{\prime},\mathbf{p}_{\mathcal{R}}^{\prime})^{\prime}$. Then, for any alr transformed vector $\mathbf{z}$ based on a reference in the set of candidate references $\mathbf{p}_{\mathcal{R}}$, the $|\mathcal{R}^{c}|\times|\mathcal{R}^{c}|$ upper-left sub-matrix of the (estimated) inverse covariance matrix of $\mathbf{z}$ is invariant with respect to the choice of the reference.

In the following subsections, we will incorporate the reference-invariance property into the estimation of a sparse inverse covariance matrix for compositional count data, such as the OTU abundance data in microbiome research.

Consider an OTU abundance data set with $n$ independent samples, each of which composes observed counts of $K+2$ taxa, denoted by $\mathbf{x}_{i}=(x_{i,1},\ldots,x_{i,K+2})^{\prime}$ for the $i$-th sample, $i=1,\ldots,n$. Due to the compositional property of the data, the sum of all counts for each sample $i$ is a fixed number, denoted by $M_{i}$. Naturally, a multinomial distribution is imposed on the observed counts as

$$\mathbf{x}_{i}|\mathbf{p}_{i}\sim\text{Multinomial}(M_{i},\mathbf{p}_{i}),$$

(1)

where $\mathbf{p}_{i}=(p_{i,1},\ldots,p_{i,K+2})^{\prime}$ are the multinomial probabilities with $\sum_{k=1}^{K+2}p_{i,k}=1$.

In addition, we choose one taxon, without loss of generality, the $(K+2)$-th taxon as the reference and then apply the alr transformation (Aitchison, 1986) on the multinomial probabilities as follows

$$\mathbf{z}_{i}=\left[\log\left(\frac{p_{i,1}}{p_{i,K+2}}\right),\ldots,\log\left(\frac{p_{i,K}}{p_{i,K+2}}\right),\log\left(\frac{p_{i,K+1}}{p_{i,K+2}}\right)\right]^{\prime},\ i=1,\ldots,n.$$

(2)

Further assume that $\mathbf{z}_{i}$’s follow an i.i.d. multivariate normal distribution

$$\mathbf{z}_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(\boldsymbol{\mu},\boldsymbol{\Sigma}),\ i=1,\ldots,n,$$

(3)

where $\boldsymbol{\mu}$ is the mean, $\boldsymbol{\Sigma}$ is the covariance matrix, and $\boldsymbol{\Omega}=\boldsymbol{\Sigma}^{-1}$ is the inverse covariance matrix. The above model in (1)–(3) is called a logistic normal multinomial model and has been applied to analyze the microbiome abundance data (Xia et al., 2013).

Tian et al. (2022) proposed a method called compositional graphical lasso that aims to find a sparse estimator of the inverse covariance matrix $\boldsymbol{\Omega}$, in which the following objective function is minimized

	$\displaystyle\ell(\mathbf{z}_{1},\ldots,\mathbf{z}_{n},\boldsymbol{\mu},\boldsymbol{\Omega})={}$	$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}\left[\mathbf{x}_{i,-(K+2)}^{\prime}\mathbf{z}_{i}-M_{i}\log\{\mathbf{1}^{\prime}\exp(\mathbf{z}_{i})+1\}\right]$
		$\displaystyle-\frac{1}{2}\log[\det(\boldsymbol{\Omega})]+\frac{1}{2n}\sum_{i=1}^{n}(\mathbf{z}_{i}-\boldsymbol{\mu})^{\prime}\boldsymbol{\Omega}(\mathbf{z}_{i}-\boldsymbol{\mu})+\lambda\|\boldsymbol{\Omega}\|_{1},$		(4)

where $\mathbf{x}_{i,-(K+2)}=(x_{i,1},\ldots,x_{i,K+1})^{\prime}$ and $\mathbf{1}=(1,\ldots,1)^{\prime}$. The above objective function has two parts: The first term in (4) is the negative log-likelihood of the multinomial distribution in (1) and the remaining terms are the regular objective function of graphical lasso for the multivariate normal distribution in (3) regarding $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$ as known quantities.

Similar to 2.1, if we choose another taxon, for simplicity of illustration, the $(K+1)$-th taxon as the reference, then the alr transformation in (2) becomes

$$\mathbf{z}_{i,p}=\left[\log\left(\frac{p_{i,1}}{p_{i,K+1}}\right),\ldots,\log\left(\frac{p_{i,K}}{p_{i,K+1}}\right),\log\left(\frac{p_{i,K+2}}{p_{i,K+1}}\right)\right]^{\prime}.$$

As in Sections 2.1, $\mathbf{z}_{i,p}=\mathbf{Q}_{p}\mathbf{z}_{i}$. Therefore, $\mathbf{z}_{i,p}\stackrel{{\scriptstyle iid}}{{\sim}}N(\boldsymbol{\mu}_{p},\boldsymbol{\Sigma}_{p})$, $i=1,\ldots,n$, where $\boldsymbol{\mu}_{p}=\mathbf{Q}_{p}\boldsymbol{\mu}$, $\boldsymbol{\Sigma}_{p}=\mathbf{Q}_{p}\boldsymbol{\Sigma}\mathbf{Q}_{p}^{\prime}$, and $\boldsymbol{\Omega}_{p}=(\mathbf{Q}_{p}^{\prime})^{-1}\boldsymbol{\Omega}\mathbf{Q}_{p}^{-1}$. The reference-invariance property in 2.1 implies that $\boldsymbol{\Omega}_{1:K,1:K}=\boldsymbol{\Omega}_{p,1:K,1:K}$.

The different choice of the reference also leads to a different objective function for the compositional graphical lasso method (Comp-gLASSO) as follows

	$\displaystyle\ell_{p}(\mathbf{z}_{1,p},\ldots,\mathbf{z}_{n,p},\boldsymbol{\mu}_{p},\boldsymbol{\Omega}_{p})={}$	$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}\left[\mathbf{x}_{i,-(K+1)}^{\prime}\mathbf{z}_{i,p}-M_{i}\log\{\mathbf{1}^{\prime}\exp(\mathbf{z}_{i,p})+1\}\right]$
		$\displaystyle-\frac{1}{2}\log[\det(\boldsymbol{\Omega}_{p})]+\frac{1}{2n}\sum_{i=1}^{n}(\mathbf{z}_{i,p}-\boldsymbol{\mu}_{p})^{\prime}\boldsymbol{\Omega}_{p}(\mathbf{z}_{i,p}-\boldsymbol{\mu}_{p})+\lambda\|\boldsymbol{\Omega}_{p}\|_{1},$		(5)

where $\mathbf{x}_{i,-(K+1)}=(x_{i,1},\ldots,x_{i,K},x_{i,K+2})^{\prime}$ and $\mathbf{1}=(1,\ldots,1)^{\prime}$. Comparing (4) and (5), their first terms are the same as they are both equal to the negative log-likelihood of the multinomial distribution: $-\frac{1}{n}\sum_{i=1}^{n}\sum_{k=1}^{K+2}x_{i,k}\log p_{i,k}$. In addition, from Aitchison (1986), $\det(\boldsymbol{\Omega})=\det(\boldsymbol{\Omega}_{p})$ and $\sum_{i=1}^{n}(\mathbf{z}_{i}-\boldsymbol{\mu})^{\prime}\boldsymbol{\Omega}(\mathbf{z}_{i}-\boldsymbol{\mu})=\sum_{i=1}^{n}(\mathbf{z}_{i,p}-\boldsymbol{\mu}_{p})^{\prime}\boldsymbol{\Omega}_{p}(\mathbf{z}_{i,p}-\boldsymbol{\mu}_{p})$ as known properties of the alr transformation. However, the $L_{1}$ penalties in (4) and (5) are different because $\boldsymbol{\Omega}$ is not necessarily equal to $\boldsymbol{\Omega}_{p}$. The reference-invariance property only implies that $\boldsymbol{\Omega}_{1:K,1:K}=\boldsymbol{\Omega}_{p,1:K,1:K}$.

Motivated by the reference-invariance property, we can impose the $L_{1}$ penalties only on the invariant entries of $\boldsymbol{\Omega}$ instead of all entries of $\boldsymbol{\Omega}$ as in (4), which leads to

	$\displaystyle\ell_{inv}(\mathbf{z}_{1},\ldots,\mathbf{z}_{n},\boldsymbol{\mu},\boldsymbol{\Omega})={}$	$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}\left[\mathbf{x}_{i,-(K+2)}^{\prime}\mathbf{z}_{i}-M_{i}\log\{\mathbf{1}^{\prime}\exp(\mathbf{z}_{i})+1\}\right]$
		$\displaystyle-\frac{1}{2}\log[\det(\boldsymbol{\Omega})]+\frac{1}{2n}\sum_{i=1}^{n}(\mathbf{z}_{i}-\boldsymbol{\mu})^{\prime}\boldsymbol{\Omega}(\mathbf{z}_{i}-\boldsymbol{\mu})+\lambda\|\boldsymbol{\Omega}_{1:K,1:K}\|_{1}.$		(6)

With the previous arguments, we showed that $\ell_{\text{inv}}(\mathbf{z}_{1},\ldots,\mathbf{z}_{n},\boldsymbol{\mu},\boldsymbol{\Omega})$ is reference-invariant; in other words, the objective function $\ell_{\text{inv}}$ stays the same regardless of whether the $(K+1)$-th or the $(K+2)$-th taxa is selected as the reference. This is summarized in the following theorem.

We call $\ell_{inv}(\mathbf{z}_{1},\ldots,\mathbf{z}_{n},\boldsymbol{\mu},\boldsymbol{\Omega})$ in (6) the reference-invariant compositional graphical lasso objective function. To obtain a sparse estimator of $\boldsymbol{\Omega}$, we minimize the objective function $\ell_{inv}(\mathbf{z}_{1},\ldots,\mathbf{z}_{n},\boldsymbol{\mu},\boldsymbol{\Omega})$ with respect to $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$, $\boldsymbol{\mu}$, and $\boldsymbol{\Omega}$. We call this estimation approach the reference-invariant compositional graphical lasso (Inv-Comp-gLASSO) method.

In general, suppose we have selected a set of taxa $\mathbf{x}_{\mathcal{R}}$ as “candidate references” and write $\mathbf{x}=(\mathbf{x}_{\mathcal{R}^{c}}^{\prime},\mathbf{x}_{\mathcal{R}}^{\prime})^{\prime}$. Then, the reference-invariant objective function becomes

	$\displaystyle\ell_{\text{inv}}(\mathbf{z}_{1},\ldots,\mathbf{z}_{n},\boldsymbol{\mu},\boldsymbol{\Omega})={}$	$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}\left[\mathbf{x}_{i,\mathcal{R}^{c}}^{\prime}\mathbf{z}_{i}-M_{i}\log\{\mathbf{1}^{\prime}\exp(\mathbf{z}_{i})+1\}\right]$
		$\displaystyle-\frac{1}{2}\log[\det(\boldsymbol{\Omega})]+\frac{1}{2n}\sum_{i=1}^{n}(\mathbf{z}_{i}-\boldsymbol{\mu})^{\prime}\boldsymbol{\Omega}(\mathbf{z}_{i}-\boldsymbol{\mu})+\lambda\|\boldsymbol{\Omega}_{\mathcal{R}^{c},\mathcal{R}^{c}}\|_{1}.$		(7)

In other words, (7) is invariant regardless of which reference is selected in the set of candidate references, and so is the invariant part of its minimizer.

It is noteworthy that the trick we played in defining the reference-invariant version of Comp-gLASSO is to revise the penalty term from the regular lasso penalty on the whole inverse covariance matrix to that only on the invariant part of the inverse covariance matrix. Using the same trick, we can define the reference-invariant version of other methods such as reference-invariant graphical lasso (Inv-gLASSO). The objective function of Inv-gLASSO is defined as follows when $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$ are observed instead of $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$:

$$\ell_{\text{inv}}(\boldsymbol{\mu},\boldsymbol{\Omega})=-\frac{1}{2}\log[\det(\boldsymbol{\Omega})]+\frac{1}{2n}\sum_{i=1}^{n}(\mathbf{z}_{i}-\boldsymbol{\mu})^{\prime}\boldsymbol{\Omega}(\mathbf{z}_{i}-\boldsymbol{\mu})+\lambda\|\boldsymbol{\Omega}_{\mathcal{R}^{c},\mathcal{R}^{c}}\|_{1}.$$

(8)

The objective function (7) includes naturally three sets of parameters $(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})$, $\boldsymbol{\mu}$, and $\boldsymbol{\Omega}$, which motivates us to apply a block coordinate descent algorithm. A block coordinate descent algorithm minimizes the objective function iteratively for each set of parameters given the other sets. Given the initial values $(\mathbf{z}_{1}^{(0)},\ldots,\mathbf{z}_{n}^{(0)})$, $\boldsymbol{\mu}^{(0)}$, and $\boldsymbol{\Omega}^{(0)}$, a block coordinate algorithm repeats the following steps cyclically for iteration $t=0,1,2,\ldots$ until the algorithm converges.

1. 1.

   Given $\boldsymbol{\mu}^{(t)}$ and $\boldsymbol{\Omega}^{(t)}$, find $(\mathbf{z}_{1}^{(t+1)},\ldots,\mathbf{z}_{n}^{(t+1)})$ that maximizes (7).

2. 2.

   Given $(\mathbf{z}_{1}^{(t+1)},\ldots,\mathbf{z}_{n}^{(t+1)})$ and $\boldsymbol{\Omega}^{(t)}$, find $\boldsymbol{\mu}^{(t+1)}$ that maximizes (7).

3. 3.

   Given $(\mathbf{z}_{1}^{(t+1)},\ldots,\mathbf{z}_{n}^{(t+1)})$ and $\boldsymbol{\mu}^{(t+1)}$, find $\boldsymbol{\Omega}^{(t+1)}$ that maximizes (7).

Except for the mild modification on the objective function, the details of the algorithm is essentially the same as the one described in Tian et al. (2022).

To illustrate the reference-invariance property under the aforementioned framework, we conduct a simulation study and evaluated the performance of Inv-Comp-gLASSO as well as Inv-gLASSO.

Following Tian et al. (2022), we generate three types of inverse covariance matrices $\boldsymbol{\Omega}=(\omega_{kl})_{1\leq k,l\leq K+1}$ as follows:

1. 1.

   Chain: $\omega_{kk}=1.5$, $\omega_{kl}=0.5$ if $|k-l|=1$, and $\Omega_{kl}=0$ if $|k-l|>1$. Every node is connected to the adjacent node(s), and therefore the degree is $2$ for all but two nodes.

2. 2.

   Random: $\omega_{kl}=1$ with probability $3/(K+1)$ for $k\neq l$. Every two nodes are connected with a fixed probability, and the expected degree is the same for all nodes.

3. 3.

   Hub: Nodes are randomly partitioned into $\lceil(K+1)/20\rceil$ groups, and there’s one “hub node” in each group. For the other nodes in the group, they are only connected to the hub node but not each other. There’s no connection among groups. The degree of connectedness is much higher for the hub nodes , and is $1$ for the rest of the nodes.

In the simulations, we also vary two other factors that play a crucial role in the performances of the methods:

1. 1.

   Sequencing depth. $M_{i}$’s are simulated from $\text{Uniform}(20K,40K)$ or $\text{Uniform}(100K,200K)$, denoted by “low” and “high” sequencing depth.

2. 2.

   Compositional variation. For each aforementioned inverse covariance matrix $\boldsymbol{\Omega}$ (“low” compositional variation), we also divide each of them by a factor of $5$ to obtain another set of inverse covariance matrices, i.e., $\boldsymbol{\Omega}$/5 (“high” compositional variation).

The data are simulated from the logistic normal multinomial distribution in (1)–(3). In detail, $\mathbf{z}_{i}\sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are first generated independently for $i=1,\ldots,n$; then, the softmax transformation (the inverse of the alr transformation) was applied to get the multinomial probabilities $\mathbf{p}_{i}$ with the $(K+2)$-th entry serving as the true reference; last, the multinomial random variables $\mathbf{x}_{i}$ were simulated from $Multinomial(M_{i};\mathbf{p}_{i})$, for $i=1,\ldots,n$. We set $n=100$ and $K=49$ throughout the simulations.

The simulation results are based on 100 replicates of the simulated data. Both Inv-Comp-gLASSO and Inv-gLASSO are applied with two choices of reference, the $(K+1)$-th entry serving as the false reference and the $(K+2)$-th entry serving as the true reference. and only the reference-invariant sub-network $\boldsymbol{\Omega}_{1:K,1:K}$ is used in the evaluations. For Inv-gLASSO, we estimate $\mathbf{p}_{1},\ldots,\mathbf{p}_{n}$ with $\mathbf{x}_{1}/M_{1},\ldots,\mathbf{x}_{n}/M_{n}$, and performe the alr transformation to get the estimates of $\mathbf{z}_{1},\ldots,\mathbf{z}_{n}$, which are denoted by $\tilde{\mathbf{z}}_{1},\ldots,\tilde{\mathbf{z}}_{n}$. We then apply Inv-gLASSO to $\tilde{\mathbf{z}}_{1},\ldots,\tilde{\mathbf{z}}_{n}$ directly to find the inverse covariance matrix $\boldsymbol{\Omega}$, which also serves as the starting value for Inv-Comp-gLASSO. For both methods, we implement them with a common sequence of $70$ tuning parameters of $\lambda$.

We empirically validate the invariance property of Inv-Comp-gLASSO and Inv-gLASSO by comparing the estimators of the two sub-networks $\boldsymbol{\Omega}_{1:K,1:K}$ resulted from choosing the true and false reference separately in each method, which have been shown to be theoretically invariant in Section 2. The comparison is assessed under four criteria as follows.

1. 1.

   Normalized Manhattan Similarity

   For two matrices $\mathbf{A}$ and $\mathbf{B}$, we define the normalized Manhattan similarity (NMS) as

   $$\text{NMS}(\mathbf{A},\mathbf{B})=1-\frac{\norm{\mathbf{A}-\mathbf{B}}_{1}}{\norm{\mathbf{A}}_{1}+\norm{\mathbf{B}}_{1}},$$

   where $\norm{\cdot}_{1}$ represents the entrywise $L_{1}$ norm of a matrix. Note that $0\leq\text{NMS}\leq 1$ due to the non-negativity of norms and the triangle inequality.

2. 2.

   Jaccard Index
   For two networks with the same nodes, denote their set of edges by $\mathcal{A}$ and $\mathcal{B}$. Then the Jaccard Index (Jaccard, 1901) is defined as follows:

   $$J(\mathcal{A},\mathcal{B})=\frac{|\mathcal{A}\cap\mathcal{B}|}{|\mathcal{A}\cup\mathcal{B}|}.$$

   Obviously, it also holds that $0\leq J(\mathcal{A},\mathcal{B})\leq 1$.

3. 3.

   Normalized Hamming Similarity
   In the context of network comparison, the normalized Hamming similarity for two adjacency matrices $\mathbf{A}$ and $\mathbf{B}$ with the same $N$ nodes are defined as follows (Hamming, 1950)

   $$H(\mathbf{A},\mathbf{B})=1-\frac{\norm{\mathbf{A}-\mathbf{B}}_{1}}{N(N-1)},$$

   where $\norm{\cdot}_{1}$ denotes the entrywise $L_{1}$ norm of a matrix. Since there are at most $N(N-1)$ $1$’s in an adjacency matrix with $N$ nodes, this metric is also between $0$ and $1$.

4. 4.

   ROC Curve
   The ROC curves on which true positive rate (TPR) and false positive rate (FPR) are plotted and also compared between the choices of true and false reference.