A Robust Spearman Correlation Coefficient Permutation Test

In this work, we show that Spearman's correlation coefficient test about $H_0:\rho_s=0$ found in most statistical software packages is theoretically incorrect and performs poorly when bivariate normality assumptions are not met or the sample size is small. The historical works about these tests make an unverifiable assumption that the approximate bivariate normality of original data justifies using classic approximations. In general, there is common misconception that the tests about $\rho_s=0$ are robust to deviations from bivariate normality. In fact, we found under certain scenarios violation of the bivariate normality assumption has severe effects on type I error control for the most commonly utilized tests. To address this issue, we developed a robust permutation test for testing the general hypothesis $H_0: \rho_s=0$. The proposed test is based on an appropriately studentized statistic. We will show that the test is theoretically asymptotically valid in the general setting when two paired variables are uncorrelated but dependent. This desired property was demonstrated across a range of distributional assumptions and sample sizes in simulation studies, where the proposed test exhibits robust type I error control across a variety of settings, even when the sample size is small. We demonstrated the application of this test in real world examples of transcriptomic data of the TCGA breast cancer patients and a data set of PSA levels and age.


Introduction
The concept of correlation and regression was originally conceived by Galton when studying how strongly the characteristics of one generation of living things manifested in the following generation [1].The ideas prompting the development of more mathematically rigorous treatment of correlation were developed by Karl Pearson in 1896, which yielded the well-known Pearson Product Moment Correlation Coefficient [2] given as where X and Y are two random variables from a non-degenerative joint distribution F XY , Cov(X, Y ) denotes the covariance, and µ X and µ Y , σ X and σ Y are the population means and standard deviations, respectively.If we let (X 1 , Y 1 ), (X 2 , Y 2 ), . . ., (X n , Y n ) denote n paired i.i.d.observations, then the sample Pearson correlation coefficient is given as Shortly thereafter Pearson's work was published Spearman introduced the rank correlation coefficient in 1904, with the advantages of being robust to extreme values and disparities of distributions between two variables [3].It should be noted however that K. Pearson, in his biography of Galton, says that the latter "dealt with the correlation of ranks before he even reached the correlation of variates, i.e. about 1875", but Galton apparently published nothing explicitly [4].Mathematically, Spearman's correlation coefficient is defined as the Pearson correlation coefficient on the ranks of (X 1 , Y 1 ), (X 2 , Y 2 ), . . ., (X n , Y n ) and denoted as where a i = Rank(X i ) and b i = Rank(Y i ) are the ranks of X i and Y i , respectively, i = 1, 2, • • • , n.In general, when discussing the Spearman correlation coefficient little attention is given to its population measure.However, if we consider that a i /n converges to F (X i ), then one may consider the population measure linked to Spearman's sample correlation coefficient as, ) where F X and F Y are the marginal cumulative distribution functions (CDFs) for X and Y , respectively.The sample estimator of ρ s can be obtained by replacing the original observations with their ranks in Equation 2, For samples from a bivariate normal population, there is also a known relation between the Spearman and Pearson correlation coefficients [5], which is While Pearson's ρ measures the linear relationship between two random variables it is often described that Spearman's ρ s measures a monotonic association between X and Y , thus it may be considered more general measure of association, albeit it does measure the linear association between F (X) and F (Y ).
Under bivariate normality assumptions the transformed Z statistic approximately follows normal distribution N (0, 1.06 n−3 ) under H 0 .For small sample size scenarios, naive permutation tests are also often used, where X and Y are randomly shuffled separately to simulate the sample distribution of r s under H 0 , which is an exact test for testing the independence between X and Y under the null exchangeability assumptions, but may be an invalid test for testing H 0 : These tests are so widely used that they are often the default options in common statistical software packages such as R [6] and SAS [7].However, there is little discussion that these tests relies on the untenable assumption that the underlying sample distribution of the ranks has a bivariate normal distribution, which is in fact an impossibility.Even among those who noted this assumption, there is a misconception that the above tests are robust to such deviations because Spearman's ρ s is rank based.This is exemplified in a discussion by Feller et al. in their article [8], Conversely, starting from any bivariate distribution φ(X, Y ) we can always find monotonic transformations X = f (x), Y = g(y) to standardized normal variates x and y.The resulting bivariate distribution ψ(X, Y ) will not necessarily be bivariate normal, but we think it likely that in practical stations it would not differ greatly from this form.This is a field in which further investigation would be of considerable interest.
However, as we will show in Section 3, all the commonly used tests about ρ s = 0 as discussed above, including the naive permutation test, are not even asymptotically valid when the non-exchangeabilty assumptions are violated under H 0 .In some cases, the type I error can severely drift away from the desired level as the sample size increases!A undesirable feature that is more notable in the era of "big-data".Another variation of this approach is the Fisher-Yates coefficient, which transforms the original X and Y to their corresponding normal quantiles before the testing [9].Although the marginal distributions are of the transformed variates take a pseudo normal form, the joint normality of these transformed values is not guaranteed.
In terms of our modified permutation test it is important to note the classic large sample result in Serfling where the "distribution free" large sample normal approximation for the sampling distribution for Pearson's sample correlation coefficient r is derived using the multivariate delta method [4].This method guarantees type I error converges to α when n → ∞ given finite fourth moments.A straightforward way to obtain a similar result for Spearman's correlation is given by replacing the i.i.d sample with corresponding ranks.The test is asymptotically valid because the ranks are asymptotically independent, as we will discuss in Section 2.Even though large sample approximations about these estimators are asymptotically valid they tend to suffer inflated type I errors in the small sample setting, e.g.n < 50.
To address this issue, we propose a studentized permutation test for Spearman's correlation ρ s , which extends the work or Diccicio and Romano for Pearson's correlation coefficient ρ [10].We will show that the proposed test is asymptotically valid under general assumptions and is exact under exchangeability assumptions when ), i.e. more simply when X and Y are independent.We show that our newly proposed test has robust Type I error controls type I error control even when the sample distribution is dependent (non-exchangeability) and non-normal.Importantly, the type I error is well controlled when the sample size is small.This will be illustrated by a set of simulation studies.Finally, we will demonstrate the application of this test in real world examples of transcriptomic data of TCGA breast cancer patients, as well as a data set of PSA levels and age.

Background information
In this section, we start by reviewing the robust permutation test for Pearson's correlation coefficient as proposed by Diciccio and Romano [10].Towards this end, we define G n to be the set of all permutations π of {1, . . ., n}.For testing independence between two random variables X and Y , the permutation distribution of any given test statistic where Y n π represents {Y π(1) , . . ., Y π(n) }.In this setting, the permutation G n is all possible pairwise combinations between X n and Y n .A level α one-sided permutation test rejects if T n (X n , Y n π ) is larger than the 1 − α quantile of the permutation distribution.The permutation test is exact when exchangeability assumptions hold, that is, the distribution of (X n , Y n ) is invariant under the group of transformations G n .The test using the Pearson correlation coefficient ρ is exact when used a metric of dependence for testing the null hypothesis of independence given as where P X and P Y are marginal distributions of X and Y , respectively.The null hypothesis of independence is not equivalent to the test about zero correlation given as H 0 : ρ = 0 with the exception of limiting assumptions such as the data are distributed as bivariate normal random variables.In other words, in the general setting two random variables can be dependent but uncorrelated.In such cases, DiCiccio and Romano [10] have shown that, with finite fourth moments, the permutation distribution of ρ converges to N (0, 1), but its sampling distribution converges to N (0, τ 2 ), where τ 2 = µ 22 µ 20 µ 02 , and . Thus the test will not be level α unless τ = 1.In light of this result, DiCiccio and Romano proposed a studentized correlation test statistic, which has been shown to control Type I error asymptotically at α when two random variables are dependent but uncorrelated [10].Specifically, the studentized statistic is defined as S n = √ nρ n /τ n , where The permutation distribution and sampling distribution of S n both converge to the standard normal distribution asymptotically.It should be noted that even though the results presented in DiCiccio and Romano [10] are based on large sample approximations, the behavior of this test for small to moderate sample sizes is quite good as born out in their simulation results.

Spearman's permutation correlation test
Spearman's coefficient ρ s (X, Y ) is the Pearson correlation coefficient of the ranks of X and Y , that is ρ(a, b).When there are ties in the data, their ranks are typically taken as an average.Unlike Pearson's correlation coefficient ρ, which measures the linear relationship between two random variables, Spearman's correlation coefficient ρ s measures a monotonic association, thus is far less restrictive.Note that Spearman's correlation coefficient is also the linear measure between F (X) and F (Y ).It is also less sensitive to non-normality or extreme values.
Despite the above advantages, it is a misconception that the tests of ρ s = 0 based on the bivariate normality assumptions underlying the original data will be robust to the deviation from this assumption.In fact, when the original data are independent, their ranks will will be dependent.Thus, tests of ρ s typically suffer similar issue as for ρ.We also emphasize that "normality" refers to the joint normality as opposed to marginal normality, because two random variables that are marginally normal can have a joint non-normal distribution.Therefore, the Fisher-Yates coefficient, which back transforms a variables rank through the normal quantile function does not provide what heuristically one may consider as a simple correction.In Section 3, we will empirically show that violation of the joint normality assumption will have severe effect on type I error control.In addition, it is in fact impossible for the joint distribution of the ranks to be bivariate normal.
Our approach is to replace (X i , Y i ) with their ranks (a i , b i ) in order to develop a Spearman's correlation permutation test analog to the Pearson's correlation permutation test, with some subtle differences.The studentized permutation test of the Pearson's ρ only requires finite fourth moments and that observations are i.i.d.Although the pairs of ranks (a i , b i ) are no longer i.i.d observations, we can show that they asymptotically satisfy this condition. When Therefore, the paired observations (a i , b i ) are asymptotically i.i.d.Consequently, the exchangeability condition will hold at least asymptotically and the test will be asymptotically exact.Intuitively, when n is sufficiently large, knowing (a i , b i ) will lend little knowledge on the ranks of another pair (a j , b j ).It can also be shown that the correlation between a i and a j (i = j) is approximately − 1 n−1 .Therefore, for a sample sequence X n , the correlation matrix of the ranks converges to I n when n → ∞.
Specifically, the one sided studentized permutation test for testing H 0 : ρ s = 0 versus H 1 : ρ s > 0 is performed by the following steps.The test is implemented in the R perk (permutation tests of correlation c(k)oefficients) package, which will be available on CRAN (The Comprehensive R Archive Network, https://cran.r-project.org/) and GitHub (https://github.com/hyu-ub/perk).
• Estimate the variance of sample estimates r s by • Calculate the studentized statistic R s = r s /τ n .

Simulations
We examined the Type I error control across all of the tests introduced above using distributions commonly found in the literature for these examinations across a wide range of settings [10,11].For our simulation study, we focused on testing H 0 : ρ c = 0 versus H 1 : ρ c > 0, with sample sizes n = 10, 25, 50, 100, 200.Each simulation utilized 10, 000 Monte Carlo replications and the number of permutations used is 1, 000.We compared the t test, Fisher's Z-transformation (Fisher's Z), Fisher-Yates method, Serfling's large sample normal approximation (Asymp Norm), naive permutation test (Permute), and studentized permutation test (Stu Permute).The Type I error control for α = 0.05 was examined.The simulation scenarios 1 through 5 from DiCiccio and Romano.Two additional distributions were studied as well: 1. Multivariate normal (MVN) with mean zero and identity covariance.The results in Table 1 and Figure 2 show that the large sample asymptotic normal approximation has inflated type I error rates for all distributions when n ≤ 50.The t test, Fisher's Z test, Fisher-Yates, and naive permutation tests tend to be over-conservative for the exponential and circular distributions.While for t 4.1 , the type I error is consistently inflated.Note that, for these tests, such deviation cannot be corrected as sample size increases.Instead, they may converge to an arbitrary level, either lower or higher than α.
For MVN 1-9, we simulated a range of dependency among uncorrelated X and Y , where MVN 1 has the weakest and MVN 9 has the strongest dependency (Figure 1).The above four tests showed the type I error rate inflation becomes increasingly severe as the dependency increases.This demonstrates the failure in controlling type I error results from the data being dependent, which can occur when the underlying distribution is non-normal.
The MVN 45 is a case where the dependency of original data is remedied by using the ranks.In this case, the ranks will distribute as if it comes from a bivariate normal distribution, regardless of the distance between the centers of individual Gaussian sub-populations.Therefore, all four tests show well control of the type I error rate.
On the other hand, the studentized permutation test robustly control type I error for all distributions examined, even when the n is as small as 10.This demonstrates a clear advantage of the proposed test over all other commonly used tests for Spearman's correlation coefficient.

TCGA breast cancer data
As an illustration of our approach, we tested H 0 : ρ s = 0 versus H 1 : ρ s > 0 using The Cancer Genome Atlas (TCGA) breast cancer RNA sequencing (RNA-seq) data.The gene abundance was RSEM normalized [12].Fibroblast growth factor (FGF)2, FGF4, FGF7 and FGF20 are representative paracrine FGFs binding to heparan-sulfate proteoglycan and fibroblast growth factor receptors (FGFRs), whereas FGF19, FGF21 and FGF23 are endocrine FGFs binding to Klotho and FGFRs.FGFR1 is relatively frequently amplified and overexpressed in breast and lung cancer, and FGFR2 in gastric cancer.Moreover, FGF2 activates human dermal fibroblasts through transcriptional downregulation of the TP53 gene [13].In this application, we examine whether the transcriptomic abundance of FGFR1 is correlated with that of TP53.To investigate the performance in small sample settings, we selected 18 samples from 17 mucinous carcinoma patients.The scatter plot of log-transformed TP53 and FGFR1 abundances is shown in Figure 3 (left).The marginal normality of data was examined by Shapiro-Wilk test and the bivariate normality was examined by Henze-Zikler test.The p values of Shapiro-Wilk tests for log-transformed TP53 and FGFR1 abundances are 0.6077 and 0.0644, respectively.The p value of Henze-Zikler test is 0.3478.Although there is no statistical significant, the marginal of FGFR1 abundance likely deviates from normal distribution.
The estimated Spearman's correlation is ρs = 0.4056.Table 2 shows the results of hypothesis testing.Only the result of studentized permutation test is non-significant at α = 0.05 and suggests there is no evidence of positive correlation between TP53 and FGFR1.In fact the biology does not support a positive correlation either, since FGFR1 mediates negative regulation of TP53 by FGF2 at transcriptional level [13].Indeed, if we include all samples (n ≈ 1200) from TCGA breast cancer cohort, then all tests will fail to reject H 1 with p-values over 0.5 except for Fisher-Yates test.
Together with the results from the simulations, the result by studentized permutation test is clearly more reliable.shows the scatter plot of age versus − log(PSA).Table 2 shows that all tests rejects the H 1 and conclude there is a non-zero correlation between age and PSA.The is an example where all tests have consistent results.Although the normality tests are significant, such deviation may have been remedied by using the ranks in this specific example.

Discussion
Conventional tests of the Spearman's correlation rely on normality assumption, including t-test, Fisher's Z transformation, and naive permutation test, which fails to control Type I error rates when the assumption is violated.This was illustrated in our simulations studies (Section 3).Such defect cannot be remedied by transforming the marginal distributions such as by Fisher-Yates coefficient.Notably, the deviation from bivariate normality can result in a convergence of type I error rate to an arbitrary level when n → ∞.This indicates that, under scenarios when two random variables are uncorrelated but dependent, the type I error will not be controlled at desired level no matter how large the sample size is.On the other hand, the Serfling's test based on delta method guarantees that the type I error rate converges to α as long as the fourth order moment is finite.However, it typically suffers an inflated type I error when sample size is under 50.
In this work, we present a robust Spearman's correlation permutation test based on studentized statistic for testing H 0 : ρ s = 0 versus H 1 : ρ s > 0. The proposed approach is inspired by the work by DiCiccio and Romano [10], which was developed for Pearson's correlation.Through extensive simulation studies and real world application, we show the proposed test controls type I error even when sample size is as small as 10 and normality assumption is violated.Therefore, the test is valid in general cases.In addition, the studentized statistic can also be used for bootstrapping tests, so as to test for more general point null hypotheses [11].In conclusion, the proposed studentized permutation test should be used as a routine for testing non-zero Spearman's correlation coefficient.

Figure 1 :
Figure 1: Density plots of the a synthetic (n = 1, 000) for each simulation distributions.

Figure 3 :
Figure 3: Scatter plot of log-transformed TP53 versus FGFR1 abundance for TCGA data (left), and age versus − log(PSA) for the PSA data (right).

Table 2 :
[14]lts of testing HThe testing methods were also applied to a data set of age and baseline prostate-specific antigen (PSA) levels[14].The data consists of age and PSA levels of 480 subjects, of which 473 have complete paired observations.The sample Spearman's correlation coefficient between age and PSA is ρs = −0.1622.Since the alternative hypothesis of proposed test is H 1 : ρ s > 0, we applied a negative log transformation on PSA levels.Similar as the TCGA example, the marginal normality of data was examined by Shapiro-Wilk test, and the bivariate normality was examined by Henze-Zikler test.The p values of Shapiro-Wilk tests for log-transformed age and PSA levels are 0.0208 and 0.0301, respectively.The p value of Henze-Zikler test is < 0.0001.The results indicates the distribution is not bivariate normal.Figure 3 (right) 0 : ρ s = 0 versus H 1 : ρ s > 0 for TCGA breast cancer data and PSA data.