Predictive Data Calibration for Linear Correlation Significance Testing

Inferring linear relationships lies at the heart of many empirical investigations. A measure of linear dependence should correctly evaluate the strength of the relationship as well as qualify whether it is meaningful for the population. Pearson's correlation coefficient (PCC), the \textit{de-facto} measure for bivariate relationships, is known to lack in both regards. The estimated strength $r$ maybe wrong due to limited sample size, and nonnormality of data. In the context of statistical significance testing, erroneous interpretation of a $p$-value as posterior probability leads to Type I errors -- a general issue with significance testing that extends to PCC. Such errors are exacerbated when testing multiple hypotheses simultaneously. To tackle these issues, we propose a machine-learning-based predictive data calibration method which essentially conditions the data samples on the expected linear relationship. Calculating PCC using calibrated data yields a calibrated $p$-value that can be interpreted as posterior probability together with a calibrated $r$ estimate, a desired outcome not provided by other methods. Furthermore, the ensuing independent interpretation of each test might eliminate the need for multiple testing correction. We provide empirical evidence favouring the proposed method using several simulations and application to real-world data.


Introduction
Inferring linear relation between two observed variables is a ubiquitous task in empirical sciences. The century old Pearson's correlation coefficient (PCC) [Pearson, 1896] has become the de-facto measure to estimate bivariate linear relationships. A PCC test provides an estimate of the strength and direction of the relationship r together with a p-value. In the framework of Null Hypothesis Significance Testing (NHST), this p-value quantifies the uncertainty associated with observing a linear association (H1) when the null hypothesis (H0) of no association were true. A p-value lower than a predefined point-wise significance level α (usually set to 0.05) is taken as evidence for rejecting H0 and H1 is declared to be true. This significance testing procedure is widely adopted due to its ease of application and domain-independent nature. PCC, however, is not perfect and both the estimated strength r and p-value suffer from shortcomings. First, violation of assumptions behind PCC, such as presence of outliers and nornormality of data, leads to mischaracterization of the linear dependence [Pernet et al., 2013, Bishara and Hittner, 2015, Armstrong, 2019. Second, erroneous interpretation of a p-value as the probability of H0 being true leads to inflated Type I errors and has been criticised of its interpretation and reliability [Goodman, 2008, Gao, 2020, and blamed for the replication/reproducibility crisis [Nuzzo, 2014, Halsey et al., 2015, Anderson, 2020. In a wider context there are calls to abandon NHST altogether [Szucs and Ioannidis, 2017, Häggström, 2017, Cumming, 2014]. An associated issue is that of multiple testing (also called multiplicity) which leads to a surge in the Type I error rate (increased false positives) when multiple tests are performed simultaneously.
With large-scale data measuring thousands or even millions of variables becoming common owing to the advancements in data collection techniques and falling costs, multiple testing has become a common concern.
To alleviate the issue of incorrect estimation of the r value, there have been suggestions of alternative or specialized methods that exhibit robustness [Pernet et al., 2013, Bishara and Hittner, 2015, Wilcox, 2012. However, such methods are often designed for handling specific issues, for instance robustness to outliers using skipped-correlations [Wilcox, 2015]. Use of alternate methods, rank-based Spearman's correlation can be used which is less affected by outliers but measures monotonic dependence and not linear dependence as desired here. Other options include performing non-linear transformation of the data, e.g. Box-Cox or rankit [Bishara and Hittner, 2015].
Regarding the (mis)use of p-values, their reflected use has been advised with a caution that any interpretation alternative will also be susceptible to similar fallacies [Leek and Peng, 2015, Nickerson, 2000, Gibson, 2021, Lakens, 2021 (also see statement from the American Statistical Association [Wasserstein and Lazar, 2016]). It has been argued that point estimates can provide valuable insights and are an essential part of science for testing a claim [Mayo and Hand, 2022].
A possible remedy is "making p values work harder" [Matthews, 2021] (cf. March 2019 special issue of The American Statistician [Wasserstein et al., 2019]), e.g. estimating posterior probabilities by p value calibration [Sellke et al., 2001, Shi and Yin, 2020, Bickel, 2021a, Cabras and Castellanos, 2017, Gao, 2020, Rafi and Greenland, 2020. After this calibration the type I error probability can be interpreted as as the probability that the null hypothesis is true. As p value calibration is easy to compute and interpret it provides an attractive way to counter misinterpretation of p values.
To deal with the multiple testing, correction to the p-values arising from the individual tests is commonly employed.
Several correction methods have been proposed which rely on the number of tests and/or the distribution of p-values.
These methods differ in their strictness in penalizing p-values reflecting the trade-off between false positives and false negatives, e.g. family-wise error rate (FWER) correction provides the strictest control while false discovery rate (FDR) control is more lenient. However, any p-value correction method is not without issues and choosing a correction method is usually left to the users. As the correction is often proportional to the number of tests, the power to identify an effect decreases with the number of tests. Defining a domain for control is not straightforward, e.g. whether to exert control over data used within an experiment or an entire discipline [Trafimow and Earp, 2017]. Dependencies between the tests (usually not known a priori) require special consideration as ignoring them may increase false negatives Storey, 2008, Wilson, 2019]. Furthermore "big data" with number of tests ranging in hundreds of thousands or even millions, even lenient correction can lose power. General methods like local FDR [Efron et al., 2001] as well as domain-specific methods with higher power such as threshold free cluster enhancement (TFCE) [Smith and Nichols, 2009] for MRI images and independent hypothesis weighting (IHW) for genomics [Ignatiadis et al., 2016] have been developed. However, specialized methods can not be immediately applied to other fields. Taken together, the issues related to data assumptions, p-value interpretation and multiple testing make the trivial sounding task of inferring linear relationships rather challenging and new methods that can solve one or more of the issues are in demand.
If a hypothesis testing method could provide a p-value that can be interpreted as the probability of the null hypothesis being true and in isolation, i.e. independent of other tests, then the above mentioned issues would be resolved.
Specifically, such a test (1) would not need multiple testing correction, and (2) it will not be susceptible to dependencies between the hypotheses. A posterior probability estimate can indeed serve this purpose []. Here we take steps in this direction and propose a method which provides an alternative way to perform large scale statistical tests without having to do multiple testing corrections. The proposed correction-free method rests on the rationale that if a test can be declared generally significant then the associated p-value can be used directly for inference and does not require correction. By generally we mean if the test is repeated on other samples of data the it would yield a similar result. This implicit generalizability makes such the test independent of other tests and renders correction unnecessary, paving the way towards high-throughput analysis by identifying as many significant effects as possible, while producing relatively few false positives. In this article we propose a new machine-learning-based method for predictive data calibration.
We propose a new method where each test is treated independently by combining the notion and estimation of generalization as defined in the machine-learning field with inferential statistics. Here, we adopt the cross-validation procedure aimed at estimating the out-of-sample (i.e. generalization) prediction performance based on a finite sample. This is achieved by training a predictive model (e.g. linear regression) on a subset of the data and using the remaining subset to evaluate the prediction performance. The out-of-sample predictions then provide a generalized representation of the predicted variable. Thus, performing the statistical tests on the variables in their generalized forms provides an inference mechanism which can be interpreted as being independent of the other tests. In the proposed data calibration method dcal we replace each value with its out-of-sample prediction. The proposed method reduces the power of the test in a data-driven way solely based on the data of a particular test as opposed to traditional multiple testing correction where the correction is typically dependent on the number of tests. We claim that the proposed approach captures the true relationship and hence provides protection against the multiple comparison problem. Indeed we empirically show this to be the case by performing numerical experiments on simulated and real-world data.

Predictive Data Calibration: Proposal and Design
Our proposed method, which we call data calibration (dcal), relies on estimating generalizable information. dcal works by calibrating the data to the alternative hypothesis under test (following the convention that the default hypothesis is the null of no relationship) using out-of-sample prediction [Hastie et al., 2009]. In effect, it regularizes the data such that the samples are penalized by its deviation-i.e. the prediction error it incurs-w.r.t. the hypothesis. To this end, dcal uses a learning algorithm to replace each sample by their out-of-sample (OOS) prediction, in effect creating new data calibrated to the alternative hypothesis and then applies the classical test to this new calibrated data.
Here we focus on linear correlation and propose an alternative to the classical linear correlation test. Linear correlation is calculated for two random variables x (x 1 , x 2 , . . . , x n ∈ R) and y (y 1 , y 2 , . . . , y n ∈ R) which provides a correlation coefficient and a p-value r, p = cor(x, y). This p-value is then used as the basis for NHST.
For linear correlation, linear regression (LR) fitted using the ordinary least squares (OLS) is the natural choice as a learning algorithm. Thus for each of the two variables we calculate mutual out-of-sample predictions for each value, e.g. using the leave-one-out (LOO) cross-validation procedure: where −i indicates a subset without the i th sample and the corresponding OLS is a model induced using this subset. In words, we induce a linear regression function by leaving out a particular instance and estimate it's value by predicting it.
We then calculate the linear correlation between the predicted variables which provides the dcal version of the linear correlation test.
An illustrative and intuitive way to showcase the usefulness of the proposed method is using the Anscombe's quartet [Anscombe, 1973]. This dataset was created to demonstrate the pitfalls when using the linear correlation and consists of four pairs of variables with identical first order statistical properties and the same linear correlation but the actual  [Pernet et al., 2013]. Overall, although simplistic, analysis of the Anscombe's quartet revealed several illuminating findings (Fig. 1).
• The dcal results were different for the four datasets while, as expected, they were the same for the classical correlation tests.
• For datasets A-C, the dcal correlation values were lower and the p values values higher than their classical counterparts.
• The r dcal values decreased and the p dcal values increased from A to C as the violations by the data become more severe.
• Datasets B and C did not result in a significant dcal test with the traditionally used threshold of p < 0.05 (B would be significant at a trend-level p < 0.1).
• For dataset D, the sign of the dcal correlation reversed with a higher r dcal and a lower p dcal value. This shows that data calibration, i.e. replacing each sample with its out-of-sample prediction, reveals the underlying predictive information. When calibrated using the linear regression model, the predicted values are better suited for testing the linear relationship. Specifically, it identifies the true relationship (Fig. 1, A) and also refutes the false ones   ( Fig. 1, B-C). This suggests that the dcal test can reveal the true linear relationship between the variables and also qualify its statistical significance.
The Anscombe's dataset D exposes a critical aspect in the dcal-based test by providing a counter-intuitive result. This is because of a pitfall in the cross-validation procedure which artificially creates negative correlation when there lack of predictability. To the best of out knowledge, this issue has not been part of a publication but has been discussed in on-line forums: Perils of LOO corssvalidation 1,2 . The counter-intuitive result on the dataset D clarifies that this issue must be tackled while designing the dcal correlation test. Here we added a simple heuristic that enforces that the original correlation sign is retained after the dcal test.
Based on eq. (1), eq. (2) and the insights from the Anscombe's quartet, we devised an easy to implement algorithm for dcal correlation test (Algo. 1).
As there are several ways to perform OOS prediction [Hastie et al., 2009], it is necessary to choose a method suitable for data calibration. We evaluated three OOS methods; leave-one-out (LOO) validation, 10 times repeated 10-fold cross-validation (10x10CV), and 100 times repeated bootstrap 632 (100boot632) on various settings of two types of simulated data sets, one with no relationship at all and one with some relationships. We found that LOO and 100boot632 performed adequately while 10x10CV lead to a higher false positive rate (Fig. 2). We, therefore, deemed LOO as our default choice. However, LOO can be computationally prohibitive with very large sample sizes and 100boot632 can serve as an alternative.

Experimental Settings and Data
We used simulations to establish basic properties of our proposed method. For this reason we compared our method with various commonly used multiple testing correction methods which typically penalize the p values based on the number of tests performed. These tests included, family-wise error control using the Holm's correction (Holm), False Discover Rate control using the Benjamini-Hochberg correction (BH), permutation test using the maximum statistic across the tests (Perm_max) as well as permutation test for individual test (Perm). We employed two p-value calibration methods, a first proposed by Selke and colleagues [Sellke et al., 2001] (pcalSBB, pcal function from the pcal R package) and a second proposed by Bickel [Bickel, 2021b] (pcalBickel, implemented as pcalBickel = (1 − abs(2.7 p ln p)) p + 2 abs(2.7 p ln p)). We also included Bayesian correlation test (the correlationBF function from the R package BayesFactor) which provides a Bayes factor (BF) per test and might need multiple testing correction [Gelman et al., 2012]. There is no pre-defined and accepted cut-off for BF and two options have been proposed, 3 as moderate and 10 as strong evidence for accepting the H1 [van Doorn et al., 2020]. BF-based hypothesis testing-so called "null hypothesis Bayesian testing" (NHBT)-is prudent compared to NHST but it might be a calibration issue rather than a real disagreement between the two [Trafimow, 2003, Wetzels et al., 2011, Tendeiro and Kiers, 2019.
Instead of using the Bayes factors, we converted them to posterior probability estimates with prior probability of 0.5 (ppBF, bfactor_to_prob function from the pcal R package) [Jeffreys, 1961, Ly et al., 2016. For all the methods statistical significance for null hypotheses rejection was set to the commonly used threshold of α = 0.05. For simulation studies we did not apply any additional correction to p values from dcal, the two calibration methods or ppBF. This choice was made to gain a better understanding of "raw" p values provides by these methods. Note that use of α = 0.05 is intensely debated (e.g., see [Benjamin et al., 2018, Amrhein and Greenland, 2018, Wasserstein et al., 2019), this debate, however, is out of the scope of the current work.
The simulation studies were helpful in empirically establishing that the proposed dcal correlation test when applied to multiple hypotheses provides p values penalized for each individual test and hence can be used without additional multiple testing correction. However, the simulated data can not cover all the complexities of real world data. We, therefore, applied the proposed test together with the comparison methods to a high-dimensional RNA-seq dataset ARXHS4 [Lachmann et al., 2018] to identify co-expressed genes.

Comparison With p-value Calibration and Bayesian Methods
Using simulated data with different sample sizes and different effect sizes in the range 0.2 − 0.8, we first compared our dcal test with two p-value calibration methods (pcalSBB, pcalBickel), and the posterior probability based on the Bayes factor (ppBF). We observed that the dcal test behaved rather conservatively at lower effect sizes-i.e., high p dcal values-but then approaches the original p values at higher effect sizes (Fig. 3A). An important advantage provided by our dcal test is that it also provides conservative correlation value estimates at lower effect sizes where the Bayesian method remains close to the original values (Fig. 3B). Similarly to the p dcal values the r dcal estimates approach the true values at higher effect sizes and with increasing sample size, a behaviour desirable in real-world applications.

Simulations With and Without Correlations
We sought to establish how dcal tests perform in a controlled scenario. For this we first performed a baseline simulation where there were no correlation between the variables, i.e. they were generated completely randomly using a Gaussian distribution with µ = 0 and σ 2 = 1. We generated 100 independent x variables to test against a single dependent variable y, a usual setting in real-world scenarios. Here a well-calibrated method should not provide any significant results. Indeed, dcal performed well and did not produce many false positives (Fig. 4 A). The same was true for other methods employing explicit multiple testing correction.
We then generated data with 100 variables correlated with the outcome variable with effect size of 0.5. A good method should provide 100 rejections if it can correctly handle the multiple tests. We indeed observed that dcal, pcalSBB, and FDR generated close to 100 rejections when performing 1000 or 10000 tests (Fig. 4 B). Two methods, Holm and Perm_max, were overly conservative with the smaller sample size (N=50).
Overall, Perm and Uncorr produced many more false positives, and more tests (10000)

Gene Co-expression Analysis
We analysed RNA-seq data to to identify which of the available 32831 genes are significantly correlated with the human specific ARHGAP11B that is involved neocortex expansion [Florio et al., 2015].
The three methods were consistent such that the most conservative correction Holm (11439 significant) was subsumed by dcal (17882 significant) which in turn was a subset of FDR (23411 significant) (5 A).

Effect of outliers
We then investigated how the dcal method deals with outliers in general. To this end, we performed simulations following previous work [Pernet et al., 2013, Bishara andHittner, 2015] and compared with a robust correlation method, skipped correlation [Wilcox, 2015]. In all cases 10% of the data was sampled from the corresponding outlier distribution. For the next two kinds of outliers (univariate and bivariate), we varied the correlation between the two variables. For the univariate outliers ( Fig. 6 B), the skipped correlation estimated the r values close to the expected values at all sample sizes while both cor and dcal underestimated them with the latter being more pessimistic. Similar observations could be made while looking only at the significant correlations. While the sensitivity of all the methods increases with increasing value of true r, skipped was most sensitive followed by cor and dcal.
With bivariate outliers (Fig. 6 C), some of the aforementioned trends reversed. skipped estimates of r were close to the true value while both cor and dcal overestimated them. The sensitivity of skipped was lower while the other two methods showed near-perfect sensitivity.
In line with the Anscombe's quartet analysis (Fig. 1), dcal showed reduced sensitivity in presence of outliers with reduced effect sizes. However, with bivariate outliers dcal behaved similar to cor suggesting that the type of outliers can impact its behaviour, which can be expected of any method. The reduced sensitivity can be a desirable property to avoid false positives in presence of outliers.

Discussion and Conclusions
We proposed a novel ML-based predictive method to perform linear correlation NHST. Our method, dcal, calibrates the data using the alternative hypothesis model (here linear regression) such that each calibrated data point encodes information regarding its conformation to that hypothesis. Subsequent application of the conventional test to this calibrated data yields calibrated test results. Importantly, both the correlation estimate r dcal and the p-value p dcal are calibrated. This provides an advantage compared to the probability calibration methods which only calibrate the p values.
Analysis of the Anscombe's quartet dataset show that dcal-based correlation test can handle scenarios where both probability calibration and Bayesian methods fail (Fig. 1, Table 1). These results are further validated by several simulations as well as real-world gene co-expression data analysis. Our results regarding multiple testing suggest that the output of dcal does not require additional correction for multiple tests.
Our method has important implications for scientific and industrial applications. A major limitation of traditional multiple testing correction methods is that they become conservative with increasing number of tests. As the number of variables keep growing, e.g. number of probes in a micro-array experiment, resolution of neuorimgaing and satellite data, traditional frameworks become too conservative. dcal can potentially alleviate this limitation as p dcal can be used directly for qualifying a test as significant or not. Taken together, our proposed method can help scientific investigation by providing an easy to use alternative to traditional linear correlation analysis and multiple testing correction.
Another important consideration in the big data era is the dependence of p values on the sample size [Gómez-de Mariscal et al., 2021]. As increasing sample size leads to decreasing p-value it also increases the chances of false positives. How this dependence plays out with dcal needs further investigation.
Statistics is a founding discipline of machine learning and many algorithms and theoretical results have benefited from a close interaction with statistics. In this work, we use ML to make a contribution to statistics.

Limitations and Future Work
Here we developed a data-calibration-based test for linear correlation. However, the proposed dcal paradigm provides a framework which can be extended to other statistical tests using appropriate algorithms, for instance t-test using logistic regression. We envision a general data-calibration-based NHST framework applicable to various statistical tests which we fondly call Significance testing INcorporating Generalization (SING).
Our current algorithm (Algo. 1) uses LOO validation for obtaining out-of-sample predictions. It is possible to use other schemes like bootstrap to obtain out-of-sample predictions but here we used LOO as it provides the most unbiased estimates which is beneficial for our purpose. However, LOO can be computationally expensive for large data sets and further tests should be performed to identify optimal choices.
The p dcal values were observed to be close to that of p value calibration methods. Given the low computational requirement of p value calibration, they offer a tempting preference. However, they do not perform well on cornercases-as we showed using the Anscombe's quarter-and they do not offer calibrated estimates of the statistic estimate.
Investigating similarities and differences between these two approaches may provide novel insights and possibly faster algorithms.
We acknowledge that our use of Bayesian methods is not Bayesian in spirit as we converted the Bayes factors into posterior probability estimates and that too with a fixed prior probability of 0.5. Although resulting point estimates are non-Bayesian in spirit they serve well the comparative purpose as we intended in this study.
Taken together, our proposed method provides an alternative way to look at statistical significance testing. We hope that the conceptual and theoretical novelty helps in developing new methods making the application of the widely used NHST framework more robust and less prone to misinterpretation.