Truncated rank correlation (TRC) as a robust measure of test-retest reliability in mass spectrometry data

Comparison of Two MS Data sets and IID Bivariate samples

In Introduction, we mention that the repeatedly measured two MS samples are not from the independent and identically distributed bivariate random samples as assumed in the use of Kendall’s τ for truncated data. Here, we provide more details about the differences between the underlying models of these two data sets. For the sake of simplicity, we suppose U=(U1,U2,…,Up) and V=(V1,V2,…,Vp) are repeatedly measured log-transformed mass spectral data of p compounds. We further assume that M is a set of sites where the intensity (abundance) of analytes is measured. In this paper, we assume that each element of U_ks and V_ks are distributed as

and

where Z_ks denote shared intensities between U and V; ϵ_ik are independent random noise from (0,σϵ2); and Z_ks and ϵ_iks are also independent of one another.

If either all shared signals are zero or M is an empty set (i.e., Z_k = 0 for k ∈ M or M = ∅), the underlying distribution of two MS data is equivalent to IID bivariate distribution of (ϵ1k,ϵ2k) for k = 1, 2, … , p, where (ϵ1k,ϵ2k) are from the IID distribution of (0,σϵ2I2) and I₂ is the two-dimensional identity matrix. However, the above conditions are unusual in real MS experiments since the analytes contain at least one chemical compound and then its intensity is measured by the instrument in the well-designed MS experiment.

The other possible condition that leads to the equivalence of underlying models of two MS data and IID bivariate samples is that M={1,2,…,p} and Z_k = μ for 1 ≤ k ≤ p. Then, two MS data can be represented as Uk=μ+ϵ1k and Vk=μ+ϵ2k and (U_k, V_k) follows IID bivariate distribution (μ12,σϵ2I2), where 1₂ is a two dimensional vector of ones and I₂ is the two dimensional identity matrix. This condition is also unusual in the real MS experiments since the intensities of different chemical compounds are usually different.

In summary, the major difference between the underlying models of two MS data sets and IID bivariate samples is in the sampling unit; the sampling unit of the MS data is the p-dimensional vector and that of IID bivariate sample is a two-dimensional vector (X_k, Y_k). Note that IID bivariate samples should have identical mean but the elements in the MS data are allowed to have different means.

From this observation, the de-centered two MS data (Uk−Zk,Vk−Zk) for k = 1, 2, … , p follow the IID bivariate distribution (0,σϵ2I2). Thus, as we stated in Section 2.1, there is no reason to believe that log-transformed intensities (U_k, V_k) are correlated over k (m/z) in the sense that (Uk−Zk,Vk−Zk) are independent for k = 1, 2, … , p. To investigate this, we show the sample autocorrelation function of (Uk−μ^k) and plots of two original MS data sets X=exp⁡(U) and Y=exp⁡(V) in Figure 6, where U_k and V_k are obtained from the first and second samples (X₁ and X₂) in the real data example of Section 5.1 and μ^k=(Uk+Vk)/2. Here, to make one-dimensional data U and V from the two-dimensional observations X₁ and X₂ in Section 5.1, we summed intensities of X₁ and X₂ along with the retention time and log-transformed the data. From Figure 6, the autocorrelations are not significantly away from 0 as we explained above.

Figure 6:

Plots of the ACF of de-centered X, the two MS data (exp(X)), (exp(Y)). Blue dotted line in (a) denotes the 95% confidence intervals of the ACF. (A) ACF of de-centered X, (B) Original MS (exp(X)) and (C) Original MS (exp(Y)).

Real data example of density estimation for generating synthetic data

In this section, we describe the distribution of the MS data through the gas chromatograph-mass spectrometry (GC-MS) data from Yu et al. (2015). The samples were collected from 30 species of Korean and 27 species of Chinese Schizandra chinensis to classify the country of origin. In this experiment, the observed intensities are pre-processed to generate the chromatogram by summing intensities in the pre-determined intervals from the instrument. After pre-processing, each MS data has 3209 intensities observed at retention time points from 6:20 to 60:00. Figure 7A and B show the plots of the intensities and log-transformed intensities along the retention times, respectively. More details on the data can be found from Yu et al. (2015). Note that we use the additional MS data to describe the distribution of the MS data instead of using the data sets in Section 5 since the data sets in Section 5 are the pre-processed data sets where the intensities from noise are removed by applying the peak detection methods.

Figure 7:

Plots of the intensities (A) and log-transformed intensities (B) along the retention times of the Korean Schizandra chinensis GC/MS data. (A) Original MS and (B) Log-transformed MS.

To describe the distribution of the MS data, we show the histograms of the intensities and log-transformed intensities in Figure 8A and B, respectively. As stated in Hastings, Norton, and Roy (2002), the log-transformed intensities in Figure 8B seems to follow the mixture of the normal distributions. Thus, in order to generate synthetic data sets that resemble the real data sets, we have applied the density estimation of the finite normal mixture model for each Korean Schizandra chinensis GC-MS data with the R package mixtools (Benaglia et al. 2009). To be more specific, let X₁, … , X_p be random samples from a finite mixture F of k normal distributions with mean μ_k and variance σk2. The density of the mixture of k normal distributions for X_i can be written as

Figure 8:

Histograms of the intensities (A) and log-transformed intensities (B) along the retention times of the Korean Schizandra chinensis GC/MS data. (A) Original MS and (B) Log-transformed MS.

where θ=(π1,…,πk,μ1,…,μk,σ12,…,σk2) denotes parameters of the mixture distribution F, π_js are mixture weights that add up to unity. For a given number of components k, the parameters of the mixture distribution F are estimated by the EM algorithm (Dempster, Laird, and Rubin 1977), which is implemented in mixtools as an R function normalmixEM(). In practice, the number of components k is usually unknown. Thus, we need to select the number of components k in the mixture distribution based on the observed data. In this paper, the number of components k is chosen by the likelihood ratio test (LRT) based on the parametric boostrap (McLachlan 1987), which is implemented in mixtools as an R function boot.comp().

From the results of LRT using boot.comp() with 200 bootstrap realizations, three or four components are chosen in 76.7% of samples (23 out of 30, n₃ = 11 and n₄ = 12), where n_k denotes the number of data sets having k mixture components. Figure 9A and B show the estimated mixture densities of the three and four components, respectively. Table 6 reports the averages of the estimated mixture weights, means and standard deviations of components for k = 3,4. As reported in Table 6, the observed data sets do not have the same mixture components and can be grouped by the chosen number of components. The difference of the mixture components might be caused by the differences of the cultivation areas of the samples that were collected from several local regions in Korea. Note that this section aims at the estimation of the mixture density in order to generate a realistic synthetic data for the numerical study. We only focus on the estimated parameters of the mixture of three normal distributions (i.e., k = 3).

Figure 9:

Histograms and plots of estimated mixture density functions for the 19th (k = 3) and 24th (k = 4) samples, where k denote the number of components chosen by the LRT. (A) k = 3 and (B) k = 4.

To specify the parameters of the synthetic data, recall the underlying model of the TRC:

where Zk∼(μ(k),σz2); ϵ_k are independent random noise from (μϵ,σϵ2); and Z_ks and ϵ1ks are also independent of one another. Motivated by the mixture density estimation of the real data, we assume that the random variable U_k is from the finite mixture of three normal mixture distributions with mixture weights (π1,π2,π3), mean (μϵ,μϵ+μ1,μϵ+μ2) and variance (σϵ2,σϵ2+σz2,σϵ2+σz2). Based on the estimates in Table 6, we set

Then, we generate U₁, … , U_k from the underlying model of the TRC with the above parameters, m0≡m01+m02=⌊p∗π2⌋+⌊p∗π3⌋, μ(k)=μ1 for k=1,2,…,⌊p∗π2⌋ and μ(k)=μ2 for k=⌊p∗π2⌋+1,…,m0. In addition, we randomly permute the generated X₁, … , X_p to avoid systematic ordering of samples due to the generation scheme. Figure 10A and B show a histogram of the generated synthetic samples and a plot of the exponential value of the generated synthetic data that resembles the real data set depicted in Figure 7A and Figure 9A, respectively. We re-order the synthetic samples based on the order of the real data set shown in Figure 7 in order to reflect the site information of the retention time.

Figure 10:

Histogram of {X1,…,Xp} and plot of (exp⁡(X1),…,exp⁡(Xp)), where {X1,…,Xp} are synthetic samples from the finite mixture of three normal distributions. (A) Histogram of X and (B) Plot of exp(X).

Approximated null distributions of correlation measures by permutation

In the numerical study, we calculate p-values of the independence test based on the correlation measures using three methods including the asymptotic null distribution, the approximated null distributions generated from the permutation method, and the empirical null distribution of the simulated samples under null. As shown in the numerical study, the sizes and powers of the rank-based measures (Kendall’s tau and TRC) are well-approximated by the permutation method but the sizes and powers of the Pearson’s correlation are poorly estimated by the permutation for relatively large noise level (σ_ϵ = 1, 1.5).

In this section, we visually compare the approximated null distributions of each correlation measure for noise levels σ_ϵ = 1, 1.5. For Pearson’s correlation and Kendall’s tau are transformed the following equations to compare their asymptotic distributions and the approximated null distribution with generated samples:

In the simulation, we also consider the true null case (γ = 0) that two MS data are independent. Thus, the approximated null distribution of the null case using the permutation can be considered as a good approximation of the true null distribution. To investigate the quality of the approximation visually, we apply the kernel density estimation to the correlation measures calculated by the permutation. The estimated densities of the null distribution of Pearson’s correlation, Kendall’s tau, TRC with m₀ and TRC with m^ are depicted in Figure 11–Figure 14, respectively. From these visualizations, the null distribution of the Pearson’s correlation is poorly estimated by the permutation when the observed samples are from the non-null cases (γ > 0). On the other hand, the null distributions of the rank-based measures (Kendall’s tau and TRC) are well-estimated by the permutation. In addition, from Figure 12 and Figure 14, we can see that the mean of the null distribution of TRC tau with m^ is not centered at 0 and is close to 0.05, while the null distribution of TRC tau with m₀ is centered at 0. This could be caused by the proposed selection rule that tends to choose relatively large m value to contain more true shared signals as described in the numerical study.

Figure 11:

Estimated density of the transformed Pearson’s correlation Z_ρ from the approximated null distribution by the permutation. (A) Z_ρ (σ_ϵ = 1) and (B) Z_ρ (σ_ϵ = 1.5).

Figure 12:

Estimated density of the transformed Kendall’s tau Z_τ from the approximated null distribution by the permutation. (A) Z_τ (σ_ϵ = 1) and (B) Z_τ (σ_ϵ = 1.5).

Figure 13:

Estimated density of the TRC tau with true m₀τ^(m0) from the approximated null distribution by the permutation. (A) τ^(m0) (σ_ϵ = 1) and (B) τ^(m0) (σ_ϵ = 1.5).

Figure 14:

Estimated density of the TRC tau with m^ by the proposed selection rule τ^(m^) from the approximated null distribution by the permutation. (A) τ^(m^) (σ_ϵ = 1) and (B) τ^(m^) (σ_ϵ = 1.5).