Discovering and deciphering relationships across disparate data modalities

Understanding the relationships between different properties of data, such as whether a genome or connectome has information about disease status, is increasingly important. While existing approaches can test whether two properties are related, they may require unfeasibly large sample sizes and often are not interpretable. Our approach, ‘Multiscale Graph Correlation’ (MGC), is a dependence test that juxtaposes disparate data science techniques, including k-nearest neighbors, kernel methods, and multiscale analysis. Other methods may require double or triple the number of samples to achieve the same statistical power as MGC in a benchmark suite including high-dimensional and nonlinear relationships, with dimensionality ranging from 1 to 1000. Moreover, MGC uniquely characterizes the latent geometry underlying the relationship, while maintaining computational efficiency. In real data, including brain imaging and cancer genetics, MGC detects the presence of a dependency and provides guidance for the next experiments to conduct.

Identifying the existence of a relationship between a pair of properties or modalities is the critical initial step in data science investigations.Only if there is a statistically significant relationship does it make sense to try to decipher the nature of the relationship.Discovering and deciphering relationships is fundamental, for example, in high-throughput screening [1], precision medicine [2], machine learning [3], and causal analyses [4].One of the first approaches for determining whether two properties are related to-or statistically dependent on-each other is Pearson's Product-Moment Correlation (published in 1895 [5]).This seminal paper prompted the development of entirely new ways of thinking about and quantifying relationships (see [6,7] for recent reviews and discussion).Modern datasets, however, present challenges for dependence-testing that were not addressed in Pearson's era.First, we now desire methods that can correctly detect any kind of dependence between all kinds of data, including high-dimensional data (such as 'omics), structured data (such as images or networks), with nonlinear relationships (such as oscillators), even with very small sample sizes as is common in modern biomedical science.Second, we desire methods that are interpretable by providing insight into how or why they discovered the presence of a statistically significant relationship.Such insight can be a crucial component of designing the next computational or physical experiment.
While many statistical and machine learning approaches have been developed over the last 120 years to combat aspects of the first issue-detecting dependencies-no approach satisfactorily addressed the challenges across all data types, relationships, and dimensionalities.Hoeffding and Renyi proposed non-parametric tests to address nonlinear but univariate relationships [8,9].In the 1970s and 1980s, nearest neighbor style approaches were popularized [10,11], but they were sensitive to algorithm parameters resulting in poor empirical performance."Energy statistics", and in particular the distance correlation test (Dcorr), was recently shown to be able to detect any dependency with sufficient observations [12], at arbitrary dimensions [13], and structured data [14].Another set of methods, referred to a "kernel mean embedding" approaches, including the Hilbert Schmidt Independence Crite-rion (Hsic) [15,16], have the same theoretical guarantees, which is not surprising given that they are known to be equivalent to energy statistics [17,18].Both energy statistics methods and kernel methods perform very well empirically with a relatively small sample size on high-dimensional linear data, whereas another test (Heller, Heller, and Gorfine's test, Hhg) [19] perform well on low-dimensional nonlinear data.But no test performs particularly well on high-dimensional nonlinear data with typical sample sizes, which characterizes a large fraction of real data challenges in the current big data era.
Moreover, to our knowledge, no method for detecting dependence has even addressed the issue of interpretability.On the other hand, much effort has been devoted to characterizing "point cloud data", that is, summarizing certain global properties in unsupervised settings (for example, having genomics data, but no disease data).Classic examples of such approaches include Fourier [20] and wavelet analysis [21].More recently, topological and geometric data analysis compute properties of graphs, or even higher order simplices [22].Such methods build multiscale characterization of the samples, much like recent developments in harmonic analysis [23,24].However, these tools typically lack statistical guarantees under noisy observations, and are often quite computationally burdensome.We surmised that both (i) empirical performance in high-dimensional, nonlinear, low-sample size settings, and (ii) providing insight into the discovery process, could be satisfactorily addressed via extending existing dependence tests to be adaptive to the data [25].Existing tests rely on a fixed a priori selection of an algorithmic parameter, such as the kernel bandwidth [26], intrinsic dimension [24], and/or local scale [10,11].Indeed, the Achilles Heel of manifold learning has been the requirement to manually choose these parameters [27].Post-hoc cross-validation is often used to make these methods effective adaptive, but doing so adds an undesirable computational burden, and may weaken or destroy any statistical guarantees.There is therefore a need for statistically valid and computationally efficient adaptive methods.
To illustrate the importance of adapting to different kinds of relationships, consider a simple illustrative example: investigate the relationship between cloud density and grass wetness.If this relationship were approximately linear, the data might look like those in Figure 1A (top).On the other hand, if the relationship were nonlinear-such as a spiral-it might look like those in Figure 1A (bottom).Although the relationship between clouds and grass is unlikely to be spiral, spiral relationships are prevalent in nature and mathematics (for example, shells, hurricanes, and galaxies), and are canonical in evaluations of manifold learning techniques [28], thereby motivating its use here.
Under the linear relationship (top panels), when a pair of observations are close to each other in cloud density, they also tend to be close to each other in grass wetness (for example, observations 1 and 2 highlighted in black in Figure 1A, and distances between them in Figure 1B).Similarly, when a pair of observations are far from each other in cloud density, they also tend to be far from each other in grass wetness (see for example, distances between observations 2 and 3).On the other hand, consider the nonlinear (spiral) relationship (bottom panels).Here, when a pair of observations are close to each other in cloud density, they also tend to be close to each other in grass wetness (see points 1 and 2 again).However, the same is not true for large distances (see points 2 and 3).Thus, in the linear relationship, the distance between every pair of points is informative with respect to the relationship, while under the nonlinear relationship, only a subset of the distances are.
For this reason, we juxtapose nearest neighbor methods with distance/kernel methods.Specifically, for each point, we find its k-nearest neighbors for one property (e.g., cloud density), and its l-nearest neighbors for the other property (e.g., grass wetness); we call the pair (k, l) the "scale".A priori, however, we do not know which scales will be most informative.Therefore, leveraging recent ideas from multiscale analysis, we efficiently compute the distances for all scales.The test statistics, described in detail below, summarize the correlations between distances at each scale (Figure 1C), illustrating which  2 to n for both x and y).For the linear relationship, the global scale is optimal, and is the scale that Mgc selects, resulting in a p-value identical to Dcorr.For the nonlinear relationship, the optimal scale is local in both x and y, so Mgc achieves a far larger test statistic, and a correspondingly smaller and significant p-value.Thus, Mgc uniquely detects dependence and characterizes the geometry in both relationships.
scales are relatively informative about the relationship.The key, therefore, to successfully discover and decipher relationships between disparate data modalities is to adaptively determine which scales are the most informative.Doing so not only provides an estimate of whether the modalities are related, but also provides insight into how the determination was made.This is especially important in highdimensional data, where simple visualizations do not reveal relationships to the unaided human eye.
Our method, "Multiscale Graph Correlation" (Mgc, pronounced "magic"), generalized and extends previously proposed pairwise comparison-based approaches by adaptively estimating the informative scales for any relationship -linear or nonlinear, low-dimensional or high-dimensional, unstructured or structuredin a computationally efficient and statistically consistent fashion.This adapative nature of Mgc effectively guarantees equally good or better statistical performance compared to existing global methods.
Moreover, the dependence strength across all scales is informative about how Mgc determined the existence of a statistical relationship, therefore providing further guidance for subsequent experimental or analytical steps.Mgc is thus a hypothesis-testing and insight-providing approach that builds on recent developments in manifold and kernel learning (operating on pairwise comparisons) by combining them with complementary developments in nearest-neighbor search, and multiscale analyses.It is this union of disparate disciplines spanning data science that enables improved theoretical and empirical performance.We provide an R package called Mgc distributed on the Comprehensive R Archive Network (CRAN) to enable others to use this method for a wide variety of applications [29].

The Multiscale Graph Correlation Procedure
Mgc is a multi-step procedure to discover and decipher dependencies across disparate data modalities or properties, as follows (see Appendix A and [30] for details): 1. Compute two distance matrices, one consisting of distances between all pairs of one property (e.g., cloud densities, entire genomes or connectomes) and the other consisting of distances between all pairs of the other property (e.g., grass wetnesses or disease status).Then center each matrix (by subtracting its overall mean, the column-wise mean from each column, and the row-wise mean from each row).Call the resulting n-by-n matrices A and B.
2. Compute the k-nearest neighbor graphs for one property, and the l-nearest neighbor graph for the other property, for all possible values of k and l.Let {G k } and {H l } be the nearest neighbor graphs for all k, l = 1, . . ., n, where G k (i, j) = 1 indicates that A(i, j) is within the k smallest values of the i th row of A. Note that this yields n 2 binary n-by-n matrices, 3. Estimate the local generalized correlation, that is, the correlation between distances restricted to only the (k, l) neighbors by summing the products of these matrices, for all values of k and l.H 0 : X and Y are independent The standard criterion for evaluating statistical tests is to compute the probability that it correctly rejects a false null hypothesis, the testing power, at a given type 1 error level.In a complementary manuscript [30], we established the theoretical properties of Mgc, including proving its validity and universal consistency for dependence testing against all distributions of finite second moments, meaning that it will reject false null hypotheses for any dependency with enough samples. .
Here, we address the empirical performance of Mgc as compared with multiple popular tests: (i) Dcorr, a popular approach from the statistics community [12,31], (ii) Mcorr, a modified version of Dcorr designed to be unbiased for sample data [13], (iii) Hhg, a distance-based test that is very powerful for detecting low-dimensional nonlinear relationships [19].(iv) Hsic, a kernel test [15] which is equivelant to Dcorr but using a different kernel [18], (v) Mantel, which is historically widely used in biology and ecology [32].(vi) RV coefficient [5,7], which is a multivariate generalization of Pearson's product moment correlation whose test statistic is the sum of the trace-norm of the cross-covariance matrix, and (vii) the Cca method, which is the largest (in magnitude) singular value of the cross-covariance matrix, and can be viewed as a different generalization of Pearson in high-dimensions that is more appropriate for sparse settings [33][34][35][36].Note that while we focus on high-dimensional settings, Appendix D shows further results in one-dimensional settings, also comparing to a number of tests that are limited to one dimension, including: (viii) Pearson's product moment correlation, (ix) Spearman's rank correlation [37], (x) Kendall's tau correlation [38], and (xi) Mic [39].Under the regularity condition that the data distribution has finite second moment, the first four tests are universally consistent, whereas the other tests are consistent only for linear or monotone relationships.
We generate an extensive benchmark suite of 20 relationships, including different polynomial (linear, quadratic, cubic), trigonometric (sinusoidal, circular, ellipsoidal, spiral), geometric (square, diamond, W-shape), and other functions.This suite includes and extends the simulated settings from previous dependence testing work [13,19,31,40,41].For many of them, we introduce high-dimensional variants, to more extensively evaluate the methods; function details are in Appendix C. The visualization of one-dimensional noise-free (black) and noisy (gray) samples is shown in Supplementary Figure E1.
For each relationship, we compute the power of each method relative to Mgc for ∼20 different dimensionalities, ranging from 1 up to 10, 20, 40, 100, or 1000.The high-dimensional relationships are more challenging because (1) they cannot be easily visualized, and (2) each dimension is designed to have less and less signal, so there are many noisy dimensions.Figure 2 shows that Mgc achieves the highest (or close to the highest) power given 100 samples for each relationship and dimensionality.Supplementary Figure E2 shows the same advantage in one-dimension with increasing sample size.

Multiplicative
Relative Testing Power for 20 Simulated High-Dimensional Settings Beyond simply discovering the existence of a relationship, the next goal is often to decipher the nature or structure of that relationship, thereby providing insight and guiding future experiments.A single scalar quantity (such as effect size) is inadequate given the vastness and complexities of possible relationships.Existing methods would require a secondary procedure to characterize the relationship, which introduces complicated "post selection" statistical quandaries that remain mostly unresolved [42].optimal scale is shown in green.In linear or close-to-linear relationships (first row), the optimal scale is global, i.e., the green dot is in the top right corner.Otherwise the optimal scale is non-global, which holds for the remaining dependencies.Moreover, similar dependencies often share similar Mgc-Maps and similar optimal scales, such as (10) logarithmic and ( 11) fourth root, the trigonometric functions in ( 12) and ( 13), ( 16) circle and ( 17) ellipse, and ( 14) square and ( 18) diamond.A visualization of each dependency is provided in Appendix Figure E1, and the Mgc-Maps for HD simulations are provided in Figure E5.
Instead, Mgc provides a simple, intuitive, and nonparametric (and therefore infinitely flexible) "map" of how it discovered the relationship.As described below, this map not only provides interpretability for how Mgc detected a dependence, it also partially characterize the geometry of the investigated relationship.
The Mgc-Map shows local correlation as a function of the scales of the two properties.More concretely, it is the matrix of t kl 's, as defined above.Thus, the Mgc-Map is an n-by-n matrix which encodes the strength of dependence for each possible scale.Figure 3 provides the Mgc-Map for all 20 different one-dimensional relationships; the optimal scales, t * , are shown with green dots.For the monotonic dependencies (1-5), the optimal scale is always the largest scale, i.e., the global one.For all nonmonotonic dependencies (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19), Mgc chooses smaller scales.Thus, a global optimal scale implies a close-to-linear dependency, otherwise the dependency is strongly nonlinear.In fact, this empirical observation led to the following theorem (which is proved in Appendix A.VI) : Theorem 1.When (X, Y ) are linearly related (meaning that Y can be constructed from X by rotation, scaling, translation, and/or reflectiMgc-Mapon), the optimal scale of Mgc equals the global scale.
Conversely, a local optimal scale implies a nonlinear relationship.
Thus, the Mgc-Map explains how Mgc discovers relationships, specifically, which pairwise comparisons are most informative, and how that relates to the geometrical form of the relationship.Note that Mgc provides the geometric characterization "for free", meaning that no separate procedure is required; therefore, Mgc provides both a valid test and information about the geometric relationship.Moreover, similar dependencies have similar Mgc-Maps and often similar optimal scales.For example, logarithmic (10) and fourth root (11), though very different functions analytically, are geometrically similar, and yield very similar Mgc-Maps.Similarly, ( 12) and ( 13) are trigonometric functions, and they share a narrow range of significant local generalized correlations.Both circle (16) and ellipse (17), as well as square (14) and diamond (18), are closely related geometrically and also have similar Mgc-Maps.This indicates that the Mgc-Map partially characterizes the geometry of these relationships, differentiating different dependence structures and assisting subsequent analysis steps.Moreover, in [30] we proved that the sample Mgc-Map (which Mgc estimates) converges to the true Mgc-Map provided by the underlying joint distribution of the data.In other words, each relationship has a specific map that characterizes it based on its joint distribution, and Mgc is able to accurately estimate it via sample observations.The existence of a population level characterization of the joint distribution strongly differentiates Mgc from previously proposed multi-scale geometric or topological characterizations of data, such as persistence diagrams [22].

Mgc Uniquely Reveals Relationships in Real Data
Geometric intuition, numerical simulations, and theory all provide evidence that Mgc will be useful for real data discoveries.Nonetheless, real data applications provide another necessary ingredient to justify its use in practice.Below, we describe several real data applications where we have used Mgc to understand relationships in data that other methods were unable to provide.

Mgc Discovers the Relationships between Brain and Mental Properties
The human psyche is of course dependent on brain activity and structure.Previous work has studied two particular aspects of our psyche: personality and creativity, developing quantitative metrics for evaluating them using structured interviews [43,44].However, the relationship between brain activity and structure, and these aspects of our psyche, remains unclear.We therefore utilized Mgc to published open access data to investigate.Mgc-Map provides a characterization of the dependence, for which the optimal scale indicates that the dependency is strongly nonlinear.Interestingly, the Mgc-Map does not look like any of the 20 images from the simulated data, suggesting that the nonlinearity characterizing this dependency is more complex or otherwise different from those we have considered so far.
Second, we investigate the relationship between diffusion MRI derived connectivity and creativity [44] (see Appendix E.II for details).The second row of Table 2 shows that Mgc is able to ascertain a dependency between the whole brain network and the subject's creativity.The Mgc-Map in Figure 4B closely resembles a linear relationship where the optimal scale is global.The performance in the real data closely matches the simulations in terms of the superiority of Mgc: the first dataset is a strongly nonlinear relationship, for which Mgc has the lowest p-value, followed by Hhg and Hsic and then all other methods; the second dataset is a close-to-linear relationship, for which global methods often perform the best while Hhg and Hsic are trailing.Moreover, Mgc detected a complex nonlinear relationship for brain activity versus personality, and a nearly linear relationship for brain network versus creativity, the only method able to make either of those claims.In a separate experiment, we assessed the frequency with which Mgc obtained false positive results using brain activity data, based on experiments from [46,47].Supplementary Figure E6 shows that Mgc achieves a false positive rate of 5% when using a significance level of 0.05, implying that it correctly controls for false positives, unlike typical parametric methods on these data.

Mgc Identifies Potential Cancer Proteomics Biomarkers
Mgc can also be useful for a completely complementary set of scientific questions: screening proteomics data for biomarkers, often involving the analysis of tens of thousands of proteins, peptides, or transcripts in multiple samples representing a variety of disease types.Determining whether there is a relationship between one or more of these markers and a particular disease state can be challenging, but is a necessary first step.We sought to discover new useful protein biomarkers from a quantitative proteomics technique that measures protein and peptide abundance called Selected Reaction Monitoring (SRM) [48] (see Appendix E.III for details).Specifically, we were interested in finding biomarkers that were unique to pancreatic cancer, because it is lethal and no clinically useful biomarkers are currently available.
The data consist of proteolytic peptides derived from the blood samples of 95 individuals harboring pancreatic (n = 10), ovarian (n = 24), colorectal cancer (n = 28), and healthy controls (n = 33).The processed data included 318 peptides derived from 121 proteins.Previously, we used these data and other techniques to find ovarian cancer biomarkers (a much easier task because the dataset has twice as many ovarian patients) and validated them with subsequent experiments [49].Therefore, our first step was to check whether Mgc could correctly identify ovarian biomarkers.Indeed, the pepetides that have been validated previously are also identified by Mgc (see Appendix E.III).Emboldened, using the same dataset, we applied Mgc to screen for biomarkers unique to pancreatic cancer.To do so, we first screened for a difference between pancreatic cancer and healthy controls, identifying several potential biomarkers.Then, we screened for a difference between pancreatic cancer and all other conditions, to find peptides that differentiate pancreatic cancer from other cancers.Figure 4C shows the p-value of each peptide assigned by Mgc, which reveals one particular protein, neurogranin, that exhibits a strong dependency specifically with pancreatic cancer.Subsequent literature searches reveal that neurogranin is a potentially valuable biomarker for pancreatic cancer because it is exclusively expressed in brain tissue among normal tissues and has not been linked with any other cancer type.In comparison, Hsic identified neurogranin as well, but it also identified another peptide; Hhg identified the same two by Hsic, and a third peptide.A literature evaluation of these additional peptides shows that they are upregulated in other cancers as well and are unlikely to be useful as a pancreatic biomarker.The rest of the global methods did not identify any markers.
We carried out a classification task using the biomarkers identified by the various algorithms, using a k-nearest-neighbor classifier to predict pancreatic cancer, and a leave-one-subject-out validation.
Figure 4D shows that the peptide selected by Mgc (neurogranin) works better than any other subset of the peptides selected by Hsic or Hhg, in terms of both fewer false positives and negatives.This analysis suggests Mgc can effectively be used for screening and subsequent classification.

Discussion
There are a number of connections between Mgc and other prominent statistical procedures that may be worth further exploration.First, Mgc can be thought of as a regularized or sparsified variant of distance or kernel methods.Regularization is central to high-dimensional and ill-posed problems, where dimensionality is larger than sample size.The connection made here between regularization and dependence testing opens the door towards considering other regularization techniques for correlationbased dependence testing, including Hhg.Second, Mgc can be thought of informally as learning a metric because it chooses amongst a set of n 2 truncated distances, motivating studying the relation- ship between Mgc and recent advances in metric learning [50].In particular, deep learning can be thought of as metric learning [51], and generative adversarial networks [52] are implicitly testing for equality, which is closely related to dependence [53].While Mgc searches over a two-dimensional parameter space to optimize the metric, deep learning searches over a much larger parameter space, sometimes including millions of dimensions.Probably neither is optimal, and somewhere between the two would be useful in many tasks.Third, energy statistics provide state of the art approaches to other problems, including goodness-of-fit [54], analysis of variance [55], conditional dependence [56,57], and feature selection [58,59], so Mgc can be adapted for them as well.In fact, Mgc can also implement a two-sample (or generally the K-sample) test [18,60,61], so further comparisons of Mgc to standard methods for two-sample testing will be interesting.Finally, although energy statistics have not yet been used for classification, regression, or dimensionality reduction, Mgc opens the door to these applications by providing guidance as to how to proceed.Specifically, it is well documented in machine learning literature that the choice of kernel, metric, or scale often has an undesirably strong effect on the performance of different machine learning algorithms [27].Mgc provides a mechanism to estimate scale that is both theoretically justified and computationally efficient, by optimizing a metric for a task wherein the previous methods lacked a notion of optimization.Nonlinear dimensionality reduction procedures, such as Isomap [62] and local linear embedding [63] for example, must also choose a scale, but have no principled criteria for doing so.Mgc could be used to provide insight into multimodal dimensionality reduction as well.
The fact that Mgc provides an estimate of the informative scales suggests several theoretical steps to extend this work.First, further theoretical guidance for choosing the optimal scale in finite samples, could possibly further improve performance.Second, because the Mgc-Maps provide insight into the geometry of dependence, theoretically determining a mapping from these maps to the set of all nonlinear functions to provide a formal characterization of the geometry of the dependency, would be of interest.
Mgc also addresses a particularly vexing statistical problem that arises from the fact that methods methods for discovering dependencies are typically dissociated from methods for deciphering them.This dissociation creates a problem because the statistical assumptions underlying the "deciphering" methods become compromised in the process of "discoverying"; this is called the "post-selection inference" problem [42].The most straightforward way to address this issue is to collect new data, which is costly and time-consuming.Therefore, researchers often ignore this fact and make statistically invalid claims.Mgc circumvents this dilemma by carefully constructing its permutation test to estimate the scale in the process of estimating a p-value, rather than after.To our knowledge, Mgc is the first dependence test to take a step towards valid post-selection inference.
As a separate next theoretical extension, we could reduce the computational space and time required by Mgc.Mgc currently requires space and time that are quadratic with respect to the number of samples, which can be costly for very large data.Recent advances in related work suggest that we could reduce computational time to close to linear [64], although with some weakening of the theoretical guarantees [65].Alternately, semi-external memory implementations would allow the running of Mgc on any data as long as the interpoint comparison matrix fits on disk rather than main memory [66][67][68][69].Another approach would be to derive an approximation to the asymptotic null distribution for Mgc, obviating the need for the permutation test, but at the cost of potential finite-sample bias.
Finally, Mgc is easy to use: it merely requires pairs of samples to run, and all the code is available in both MATLAB and as a R package on Comprehensive R Archive Network, available from https: //neurodata.io/tools/and https://github.com/neurodata/mgc(code for reproducing all the figures in this manuscript is also available from the above websites) [29].Because Mgc is open source and reproducible, and obtains near empirical dominance of other methods, Mgc is situated to be useful in a wide range of applications.We showed its value in diverse applications spanning neuroscience, which motivated this work, and an 'omics example.Applications in other domains facing similar questions of dependence, such as finance, pharmaceuticals, commerce, and security, could likewise benefit from the methodology proposed here.

Bibliography
Mensh of Optimize Science for acting as our intellectual consigliere, Julia Kuhl for help with figures, and Dr. Ruth Heller, Dr. Bert Vogelstein, Dr. Don Geman, and Dr. Yakir Reshef for insightful suggestions.

A Mathematical Details
This section contains essential mathematical details on independence testing, the notion of the generalized correlation coefficient and the distance-based correlation measure, how to compute the local correlations, and the smoothing technique.A more statistical treatment on MGC is in [30], which introduces the population version of Mgc and various theoretical properties.

A.I Testing Independence
Given pairs of observations (x i , y i ) ∈ R p × R q for i = 1, . . ., n, assume they are independently identically distributed as (X, Y ) iid ∼ F XY .If the two random variables X and Y are independent, the joint distribution equals the product of the marginals, i.e., F XY = F X F Y .The statistical hypotheses for testing independence is as follows: Given a test statistic, the testing power equals the probability of rejecting the independence hypothesis (i.e. the null hypothesis) when it is false.A test statistic is consistent if and only if the testing power increases to 1 as sample size increases to infinity.We would like a test to be universally consistent, i.e., consistent against all joint distributions.Dcorr, Mcorr, Hsic, and Hhg are all consistent against any joint distribution of finite second moments and finite dimension.
Note that p is the dimension for x's, q is the dimension for y's.For Mgc and all benchmark methods, there is no restriction on the dimensions, i.e., the dimensions can be arbitrarily large, and p is not required to equal q.The ability to handle data of arbitrary dimension is crucial for modern big data.
There also exist some special methods that only operate on one-dimensional data, such as [39,61,64], which are not generalizable to multidimensional data.

A.II Generalized Correlation
Instead of computing on the sample observations directly, most state-of-the-art dependence tests operate on pairwise comparisons, either similarities (such as kernels) or dissimilarities (such as distances).
denote the matrices of sample observations, and δ x be the distance function for x's and δ y for y's.One can then compute two n × n distance matrices Ã = {ã ij } and B = { bij }, where ãij = δ x (x i , x j ) and bij = δ y (y i , y j ).A common example of the distance function is the Euclidean metric (L 2 norm), which serves as the starting point for all methods in this manuscript.Note that we will use slightly different notations in the appendix: in the main paper a ij and b ij denote the Euclidean distance, while in the appendix they denote the centered distance with ãij and bij denoting the Euclidean distance.
Let A and B be the transformed (e.g., centered) versions of the distance matrices Ã and B, respectively.Any "generalized correlation coefficient" [37,38] can be written as: where z is proportional to the standard deviations of A and B, that is z = n 2 σ a σ b .In words, c is the global sample correlation across pairwise comparison matrices A and B, rather than the individual data samples.A generalized correlation always has the range [−1, 1], has expectation 0 under independence, and implies a stronger dependency when the correlation is further away from 0.
Traditional correlations such as the Pearson's correlation and the rank correlation can be written as generalized correlation coefficients, where A and B are derived from sample observations rather than distances.Distance-based methods like Dcorr and Mantel operate on the distance metric, which may be chosen on the basis of domain knowledge, or by default they use the Euclidean distance; then transform the resulting distance matrices Ã and B by certain centering schemes into A and B.
Hsic chooses the Gaussian kernel and computes two kernel matrices, then transform the kernel matrices Ã and B by the same centering scheme as Dcorr.For Mgc, A and B are always distance matrices (or can be transformed to distances from kernels by [17]), and we shall apply a slightly different centering scheme that turns out to equal Dcorr.
To carry out the hypothesis testing on sample data via a nonparametric test statistic, e.g., a generalized correlation, the permutation test is often an effective choice [70], because a p-value can be computed by comparing the correlation of the sample data to the correlation of the permuted sample data.The independence hypothesis is rejected if the p-value is lower than a pre-determined type 1 error level, say 0.05.Then the power of the test statistic equals the probability of a correct rejection at a specific type 1 error level.Note that Hhg is the only exception that cannot be cast as a generalized correlation coefficient, but the permutation testing is similarly effective for the Hhg test statistic; also note that the iid assumption is critical for permutation test to be valid, which may not be applicable in special cases like auto-correlated time series [71].

A.III Distance Correlation (Dcorr) and the Unbiased Version (Mcorr)
Define the row and column means of Ã by ā and similarly for b ij .For distance correlation, the numerator of Equation 1 is named the distance covariance (Dcov), while σ a and σ b in the denominator are named the distance variances.The centering scheme is important to guarantee the universal consistency of Dcorr, whereas Mantel uses a simple centering scheme and thus not universal consistent.
Let c(X, Y ) be the population distance correlation, that is, the distance correlation between the underlying random variables X and Y .Szekely et al. (2007) define the population distance correlation via the characteristic functions of F X and F Y , and show that the population distance correlation equals zero if and only if X and Y are independent, for any joint distribution F XY of finite second moments and finite dimensionality.They also show that as n → ∞, the sample distance correlation converges to the population distance correlation, that is, c(X n , Y n ) → c(X, Y ).Thus the sample distance correlation is consistent against any dependency of finite second moments and dimensionality.Of note, the distance covariance, distance variance, and distance correlation are always non-negative.Moreover, the consistency result holds for a much larger family of metrics, those of strong negative type [14].
It turns out that the sample distance correlation has a finite-sample bias, especially as the dimension p or q increases [13].For example, for independent Gaussian distributions, the sample distance correlation converges to 1 as p, q → ∞.By excluding the diagonal entries and slightly modifies the off-diagonal entries of A and B, Szekely and Rizzo (2013) [13,72,73] show that Mcorr is an unbiased estimator of the population distance correlation c(x, y) for all p, q, n, which is approximately normal even if p, q → ∞.Thus it enjoys the same theoretical consistency as Dcorr and always has zero mean under independence, which is the default choice Mgc is based on in this paper.

A.IV Local Generalized Correlations
Local generalized correlations can be thought of as further generalizations of generalized correlation coefficients.In particular, given any matrices A and B, we can define a set of local variants of them as follows.Let R(A •j , i) be the "rank" of x i relative to x j , that is, R(A •j , i) = k if x i is the k th closest point (or "neighbor") to x j , as determined by ranking the n − 1 distances to x j .Define R(B i• , j) equivalently for the Y 's, but ranking relative to the rows rather than the columns (see below for explanation).For any neighborhood size k around each x i and any neighborhood size l around each y j , we define the local pairwise comparisons: and then let a k ij = a k ij − āk , where āk is the mean of { a k ij }, and similarly for b l ij .The local variant of any global generalized correlation coefficient is defined to effectively excludes large distances: where z kl = n 2 σ k a σ l b , with σ k and σ l b is the standard deviations for the truncated pairwise comparisons.
Thus, c kl is the local sample generalized correlation at a given scale.The Mgc-Map can be constructed by computing all local generalized correlations, which allows the discovery of the optimal correlation.
For any aforementioned generalized correlation (Dcorr, Mcorr, Hsic, Mantel, Pearson), its local generalized correlations can be directly defined by Equation 3, by plugging in the respective a ij and b ij from Equation 1.

A.V Mgc as the Optimal Local Correlation
We define the multiscale graph correlation statistic as the optimal local correlation, for which the family Thus, we can write: Then the optimal scale equals all scales within R whose local correlations are as large as c * .The choice of τ is made explicit in the pseudo-code, with further discussion and justification offered in [30].

A.VI Proof for Theorem 1
Theorem 1.When (X, Y ) are linearly related (rotation, scaling, translation, reflection), the optimal scale of Mgc equals the global scale.Conversely, that.the optimal scale is local implies a nonlinear relationship.
Proof.It suffices to prove the first statement, then the second statement follows by contrapositive.
When (X, Y ) are linearly related, Y = W X + b for a unitary matrix W and a constant b up-to possible scaling, in which case the distances are preserved, i.e., y i − y It follows that Mcorr(X n , Y n ) = 1, so the global scale achieves the maximum possible correlation, and the largest connected region R is empty.Thus the optimal scale is global and Mgc

B Mgc Algorithms and Testing Procedures
Six algorithms are presented in order: 1. Algorithm C1 describes Mgc in its entirety (which calls most of the other algorithms as functions).
2. Algorithm C2 evaluates the testing power of Mgc by a given distribution.

4.
Algorithm C4 computes the p-value of Mgc by the permutation test.

5.
Algorithm C5 computes the local generalized correlation coefficient at a given scale (k, l), for a given choice of the global correlation coefficient.

6.
Algorithm C6 efficiently computes all local generalized correlations, in nearly the same running time complexity as computing one local generalized correlation.
For ease of presentation, we assume there are no repeating observations of X or Y , and note that Mcorr is the global correlation choice that Mgc builds on.
Pseudocode C1 Multiscale Graph Correlation (Mgc); requires O(n 2 max(log n, p, q, r)/T ) time, where r is the number of permutations and T is the number of cores available for parallelization.
Input: n samples of (x i , y i ) pairs, an integer r for the number of random permutations.Output: (i) MGC statistic c * , (ii) the optimal scale (k, l), (iii) the p-value p(c * ), function MGC((x i , y i ), for i ∈ [n]) (1) Calculate all pairwise distances: for i, j := 1, . . ., n do a ij = δ x (x i , x j ) δ x is the distance between pairs of x samples b ij = δ y (y i , y j ) δ y is the distance between pairs of y samples end for (

Algorithm C4 end function
Pseudocode C2 Power computation of Mgc against a given distribution.By repeatedly sampling from the joint distribution F XY , sample data of size n under the null and the alternative are generated for r Monte-Carlo replicates.The power of Mgc follows by computing the test statistic under the null and the alternative using Algorithm C3.In the simulations we use r = 10,000 MC replicates.Note that power computation for other benchmarks follows from the same algorithm by plugging in the respective test statistic.
Input: A joint distribution F XY , the sample size n, the number of MC replicates r, and the type 1 error level α.

Output:
The power β of Mgc.
for t := 1, . . ., r do 3: sample from alternative 6: end for 7: for i, j := 1, . . ., n do 8: pairwise distances under the null 9: pairwise distances under the alternative 10: end for 11: Mgc statistic under the null 12: Mgc statistic under the alternative 13: end for 14: the critical value of Mgc under the null 15: compute power by the alternative distribution 16: end function Pseudocode C3 Mgc test statistic.This algorithm computes all local correlations, take the smoothed maximum, and reports the (k, l) pair that achieves it.For the smoothing step, it: (i) finds the largest connected region in the correlation map, such that each correlation is significant, i.e., larger than a certain threshold to avoid correlation inflation by sample noise, (ii) take the largest correlation in the region, (iii) if the region area is too small, or the smoothed maximum is no larger than the global correlation, the global correlation is used instead.The running time is O(n 2 ).
Output: The Mgc statistic c * ∈ R, all local statistics C ∈ R n×n , and the corresponding local scale

All local correlations
variance of all negative local generalized correlations un-normalized local distance covariance 5: local distance variances 6: 7: center and normalize 10: end function Pseudocode C6 Compute the multiscale correlation map (i.e., all local generalized correlations) in O(n 2 log n/T ).Once the distances are sorted, the remaining algorithm runs in O(n 2 ).An important observation is that each product a ij b ij is included in c kl if and only if (k, l) satisfies k ≤ R(A •j , i) and l ≤ R(B •j , i), so it suffices to iterate through a ij b ij for i, j := 1, . . ., n, and add the product simultaneously to all c kl whose scales are no more than (R(A •j , i), R(B •j , i)).To achieve the above, we iterate through each product, add it to c kl at (kl) = (R(A •j , i), R(B •j , i)) only (so only one local scale is accessed for each operation); then add up adjacent c kl for k, l = 1, . . ., n.The same applies to all local covariances, variances, and expectations.for i, j := 1, . . ., n do iterate through all local scales to calculate each term 5: for k, l := 1, . . ., n do 23: end for 25: end function

C Simulation Dependence Functions
This section provides the 20 different dependency functions used in the simulations.We used essentially the exact same relationships as previous publications to ensure a fair comparison [31,40,41].We only made changes to add white noise and a weight vector for higher dimensions, thereby making them more difficult, to better compare all methods throughout different dimensions and sample sizes.A few additional relationships are also included.
For each sample x ∈ R p , we denote x [d] , d = 1, . . ., p as the d th dimension of the vector x.For the purpose of high-dimensional simulations, w ∈ R p is a decaying vector with w [d] = 1/d for each d, such that w T x is a weighted summation of all dimensions of x.Furthermore, U(a, b) denotes the uniform distribution on the interval (a, b), B(p) denotes the Bernoulli distribution with probability p, N (µ, Σ) denotes the normal distribution with mean µ and covariance Σ, U and V represent some auxiliary random variables, κ is a scalar constant to control the noise level (which equals 1 for one-dimensional simulations and 0 otherwise), and is a white noise from independent standard normal distribution unless mentioned otherwise.
For all of the below equations, (X, Y ) For each relationship, we provide the space of (X, Y ), and define F Y |X and F X , as well as any additional auxiliary distributions.

Step Function
where I is the indicator function, that is I(z) is unity whenever z true, and zero otherwise.6. Quadratic (X, Y ) ∈ R p × R: 13.Sine Period 16π (X, Y ) ∈ R p × R p : Same as above except θ = 16π and the noise on Y is changed to 0.5κ .

Multiplicative Noise
For each distribution, X and Y are dependent except (20); for some relationships (8,14,(16)(17)(18) they are independent upon conditioning on the respective auxiliary variables, while for others they are "directly" dependent.A visualization of each dependency with D = D y = 1 is shown in Figure E1.
For the increasing dimension simulation in the main paper, we always set κ = 0 and n = 100, with p increasing.Note that q = p for types 4, 10, 12, 13, 14, 18, 19, 20, otherwise q = 1.The decaying vector w is utilized for p > 1 to make the high-dimensional relationships more difficult (otherwise, additional dimensions only add more signal).For the one-dimensional simulations, we always set p = q = 1, κ = 1 and n = 100.

Independence
Figure E1: Visualization of the 20 dependencies at p = q = 1.For each, n = 100 points are sampled with noise (κ = 1) to show the actual sample data used for one-dimensional relationships (gray dots).For comparison purposes, n = 1000 points are sampled without noise (κ = 0) to highlight each underlying dependency (black dots).Note that only black points are plotted for type 19 and 20, as they do not have the noise parameter κ.

Multiplicative
Testing Power Difference for HD Settings Figure E4: The same power plots as in Figure E3, except the 20 dependencies are one-dimensional with noise, and the x-axis shows sample size increasing from 5 to 100.

E.I Brain Activity vs Personality
This experiment investigates whether there is any dependency between resting brain activity and personality.Human personality has been intensively studied for many decades; the most widely used and studied approach is the NEO Personality Inventory-Revised the characterized personality along five dimensions [43].This dataset consists of 42 subjects, each with 197 time-steps of resting-state functional magnetic resonance activity (rs-fMRI) activity, as well as the subject's five-dimensional "personality".Adelstein et al. [45] were able to detect dependence between the activity of certain brain regions and dimensions of personality, but lacked the tools to test for dependence of whole brain activity against all five dimensions of personality.For the five-factor personality modality, we used the Euclidean distance.
For the brain activity modality, we derived the following comparison function.For each scan, (i) run Configurable Pipeline for the Analysis of Connectomes pipeline [74] to process the raw brain images yielding a parcellation into 197 regions of interest, (ii) run a spectral analysis on each region and keep the power of band, (iii) bandpass and normalize it to sum to one, (iv) calculate the Kullback-Leibler divergence across regions to obtain a similarity matrix across comparing all regions.Then, use the normalized Hellinger distance to compute distances between each subject.

E.II Brain Connectivity vs Creativity
This experiment investigates whether there is any dependency between brain structural networks and creativity.Creativity has been extensively studied in psychology; the "creativity composite index" (CCI) is an index similar to an "intelligence quotient" but for creativity rather than intelligence [44].This dataset consists of 109 subjects, each with diffusion weighted MRI data as well as the subject's CCI.Neural correlates of CCI have previously been investigated, though largely using structural MRI and cortical thickness [44].Previously published results explored the relationship between graphs and CCI [75], but did not provide a valid test.We used Euclidean distance to compare CCI values.For the raw brain imaging data, we derived the following comparison function.For each scan we estimated brain networks from diffusion and structural MRI data via Migraine, a pipeline for estimating brain networks from diffusion data [76].We compute the distance between the graphs using the semi-parametric graph test statistic [77][78][79], embedding each graph into two dimensions and aligning the embeddings via a Procrustes analysis.

E.III Proteins vs Cancer
This experiment investigated whether there is any dependency between abundance levels of peptides in human plasma and the presence of cancers.Selected Reaction Monitoring (SRM) is a targeted quantitative proteomics technique for measuring protein and peptide abundance in complicated biological samples [48].In a previous study, we used SRM to identify 318 peptides from 33 normal, 10 pancreatic cancer, 28 colorectal cancer, and 24 ovarian cancer samples [49].Then, using other methods, we identifed three peptides that were implicated in ovarian cancer, and validated them as legitimate biomarkers with a follow-up experiment.
In this study, we performed the following five sets of tests on those data: 1. ovarian vs. normal for all proteins, 2. ovarian vs. normal for each individual protein, 3. pancreas vs. normal for all proteins, 4. pancreas vs. all others for each individual protein, 5. pancreas vs. normal for each individual protein.
is exclusively expressed only in brain tissues.Moreover, there exists strong evidence of tropomyosin 3 upregulated in other cancers [81][82][83][84].Therefore, initial literature search suggests that tropomyosin is likely falsely identified by Hhg and less useful as a pancreatic cancer marker, meaning that only Mgc identified putative pancreatic cancer biomarkers without also identifying likely false positives.
Furthermore, although neurogranin is not identified by other methods, it is always the most dependent peptide in all methods except Mic.Namely, all of Pearson, Dcorr, Mcorr, Mantel, Hhg, Hsic, and Mgc rank neurogranin as the most significant protein by p-value; the only difference is that the p-values are not significant enough for other methods after multiple testing adjustment.Also, the three peptides identified by Hhg are also the top three in Mgc; and if we further investigate the top three peptides in all methods, they always come from these three peptides, and another peptide (mitogen-activated protein kinase); the only exception is Mic, whose top three peptides do not coincide with all other correlation measures, which suggests it may detect too many false positives.Along with the classification result showing that neurogranin along has the best classification error, this experiment strongly indicate that Mgc, Hsic, Hhg are the top methods in dependency testing, able to amplify the signal, and do not detect false signals.

E.IV Mgc Does Not Inflate False Positive Rates in Screening
In this final experiment, we empirically determine that Mgc does not inflate false positive rates via a neuroimaging screening.To do so, we extend the work of Eklund et al. [46,47], where a number of parametric methods are shown to largely inflate the false positives.Specifically, we applied Mgc to test whether there is any dependency between brain voxel activities and random numbers.For each brain region, Mgc attempts to test the following hypothesis: Is activity of a brain region independent of the time-varying stimuli?Any region that is selected as significant is a false positive by construction.By testing each brain region separately, Mgc provides a distribution of false positive rates.If Mgc is valid, the resulting distribution should be centered around the significance level, which is set at 0.05 for these experiments.We considered 25 resting state fMRI experiments from the 1,000 Functional Connectomes Project consisting of a total of 1,583 subjects [85].Figure E6 shows the false positive rates of Mgc for each dataset, which are centered around the critical level 0.05, as it should be.In contrast, many standard parametric methods for fMRI analysis, such as generalized linear models, can significantly increase the false positive rates, depending on the data and pre-processing details [46,47].Moreover, even the proposed solutions to those issues make linearity assumptions, thereby limiting detection to only a small subset of possible dependence functions.

E.V Running Time Report
Here

Figure 1 :
Figure 1: Illustration of Multiscale Graph Correlation (Mgc) simulating cloud density (x i ) and grass wetness (y i ).We present two different relationships: linear (top) and nonlinear spiral (bottom; see Appendix C for simulation details).(A) Scatterplots of the raw data using 50 pairs of samples for each scenario.Samples 1, 2, and 3 (black) are highlighted; arrows show x distances between these pairs of points while their y distances are almost 0. (B) Scatterplots of all pairs of distances comparing x and y distances.Distances are linearly correlated in the linear relationship, whereas they are not in the spiral relationship.Dcorr uses all distances (gray dots) to compute its test statistic and p-value, whereas Mgc chooses the local scale and then uses only the local distances (green dots).(C) Heatmaps characterizing the strength of the generalized correlation at all possible scales (ranging from

Figure 3 :
Figure 3: The Mgc-Map characterizes the geometry of the dependence function.For each of the 20 panels, the abscissa and ordinate denote the number of neighbors for X and Y , respectively, and the color denotes the magnitude of each local correlation.For each simulation, the sample size is 60, and both X and Y are onedimensional.Each dependency has a different Mgc-Map characterizing the geometry of dependence, and the

Figure 4 :
Figure 4: Demonstration that Mgc successfully detects dependency, distinguishes linearity from nonlinearity, and identifies the most informative feature in a variety of real data experiments.(A) The Mgc-Map for brain activity versus personality.Mgc has a large test statistic and a significant p-value at the optimal scale (13, 4), while the global counterpart is non-significant.That the optimal scale is non-global implies a strongly nonlinear relationship.(B)TheMgc-Map for brain connectivity versus creativity.The image is similar to that of a linear relationship, and the optimal scale equals the global scale, thus both Mgc and Mcorr are significant in this case.(C) For each peptide, the x-axis shows the p-value for testing dependence between pancreatic and healthy subjects by Mgc, and the y-axis shows the p-value for testing dependence between pancreatic and all other subjects by Mgc.At critical level 0.05, Mgc identifies a unique protein after multiple testing adjustment.(D) The true and false positive counts using a k-nearest neighbor (choosing the best k ∈[1,10]) leave-one-out classification using only the significant features identified by each testing method on the peptide data.The peptide identified by Mgc achieves the best true and false positive rates, as compared to the peptides identified by Hsic or Hhg.
of local correlation is computed based on Euclidean distance and Mcorr transformation.Instead of taking a direct maximum, Mgc takes a smoothed maximum, i.e., the maximum local correlation of the largest connected component R such that all local correlations within R are significant.If no such region exists, Mgc defaults the test statistic to the global correlation (details in Algorithm C3).

3 :τ
= THRESHOLDING(C) find a threshold to determine large local correlations 4: for i, j := 1, . . ., n do r ij ← I(c ij > τ ) end for identify all scales with large correlation 5: R ← {r ij : i, j = 1, . . ., n} binary map encoding scales with large correlation 6: R = CONNECTED(R) largest connected component of the binary matrix 7: c * ← C(n, n) use the global correlation by default 8: k ← n, l ← n 9: if i,j r ij ≥ 2n then proceed when the significant region is sufficiently large 10: [c * , k, l] ← max(C • R) find the smoothed maximum and the respective scale 11: end if 12: end function Input: C ∈ R n×n .Output: A threshold τ to identify large correlations.

τ 2 : 3 :
← max{τ, 2/n, c nn } 17: end function Pseudocode C4 Permutation Test.This algorithm uses the random permutation test with r random permutations for the p-value, requiring O(rn 2 log n) for Mgc.In the real data experiment we always set r = 10,000.Note that the p-value computation for any other global generalized correlation coefficient follows from the same algorithm by replacing Mgc with the respective test statistic.Input: A pair of distance matrices (A, B) ∈ R n×n × R n×n , the number of permutations r, and Mgc statistic c * for the observed data.Output: The p-value pval ∈ [0, 1].1: function PERMUTATIONTEST(A, B, r, c * ) c * ≤ c * 0 [t]) compute p-value of Mgc 7: end function Pseudocode C5 Compute local test statistic at a given scale.This algorithm runs in O(n 2 ) once the rank information is provided, which is suitable for Mgc computation if an optimal scale is already estimated.But it would take O(n 4 ) if used to compute all local generalized correlations.Note that for the default Mgc implementation uses single centering, the centering function centers A by column and B by row, and the sorting function sorts A within column and B within row.By utilizing T = log(n) cores, the sorting function can be easily parallelized to take O(n 2 log(n)/T ) = O(n 2 ).Input: A pair of distance matrices (A, B) ∈ R n×n × R n×n , and a local scale (k, l) ∈ N × N. Output: The local generalized correlation coefficient c kl ∈ [−1, 1].1: function LOCALGENCORR(A, B, k, l) for Z := A, B do E Z = SORT(Z) end for parallelized sorting for Z := A, B do Z = CENTER(Z) end for center distance matrices 4:

4 .
Joint normal (X, Y ) ∈ R p × R p : Let ρ = 1/2p, I p be the identity matrix of size p × p, J p be the matrix of ones of size p × p, and Σ = I p ρJ p ρJ p (1 + 0.5κ)I p .Then (X, Y ) ∼ N (0, Σ).

Figure E6 :
Figure E6:We demonstrate that Mgc is a valid test that does not inflate the false positives in screening and variable selection.This figure shows the density estimate for the false positive rates of applying Mgc to select the "falsely significant" brain regions versus independent noise experiments; dots indicate the false positive rate of each experiment.The mean ± standard deviation is 0.0538 ± 0.0394.

Typically Requires Substantially Fewer Samples to Achieve the Same Power Across All Dependencies and Dimensions
Computing all local generalized correlations, the test statistic, and p-value requires O(n 2 log n) time, which is about the same running -time complexity as other methods.
4. Estimate the optimal local generalized correlation, t * by finding the smoothed maximum of all local generalized correlations, t kl .Smoothing avoids biases and provides Mgc with better finite-sample performance and stronger theoretical guarantees. 5. Determine whether the relationship is significantly dependent-that is, whether t * is more extreme than expected under the null-via a permutation test.The permutation procedure repeats steps 1-4 on each permutation, thereby eliminating the multiple hypothesis testing problem by only computing one overall p-value, rather than one p-value per scale, ensuring that it is a valid test (meaning that the false positive rate is properly controlled at the specified type I error rate).Running Mgc is straightforward-simply input n paired samples of two measured properties, or two dissimilarity matrices of size n × n.When, and to what extent, does Mgc outperform other approaches, and when does it not?To address this question, we formally pose the following hypothesis test (see Appendix A for details): [30]number of samples as Mgc to achieve the same power.More specifically, traditional correlation methods(Pearson,RV, Cca, Spearman, Kendall) always perform the best in monotonic simulations, distance-based methods including Mcorr, Dcorr, Mgc, Hhg and Hsic are slightly worse, while Mic and Mantel are the worst.Mgc's performance is equal to linear methods on monotonic relationships.For non-monotonic relationships, traditional correlations fail to detect the existence of dependencies, Dcorr, Mcorr, and Mic, do reasonably well, but Hhg and Mgc require the fewest samples.In the high-dimensional non-monotonic relationships that motivated this work, and is common in biomedicine, Mgc significantly outperforms other methods.In fact, the second best method is Mantel, and it requires 1.6× as many samples as Mgc to achieve the same power (in prior work we proved that Mantel is not a universally consistent test[30]).The second best test that is universally consistent (Hhg) requires nearly double as many samples as Mgc, demonstrating that Mgc could half the time and cost of experiments designed to discover relationships with a given effect size.Mgc extends previously proposed global methods, such as Mantel and Dcorr.The above experiments extended Mcorr, because Mcorr is universally consistent and an unbiased version of Dcorr [13].Supplementary Figure E3 directly compares multiscale generalizations of Mantel and Mcorr as dimension increases, demonstrating that empirically, Mgc nearly dominates its global variant for essen-Pearson / RV / Cca triple tially all dimensions and simulation settings considered here.Supplementary FigureE4shows a similar result for one-dimensional settings while varying sample size.Thus, not only does Mgc empirically nearly dominate existing tests, it is a framework that one can apply to future tests to further improve their performance.Table1:The median sample size for each method to achieve power 85% at type 1 error level 0.05, grouped into monotone (type 1-5) and non-monotone relationships (type 6-19) for both one-and ten-dimensional settings, normalized by the number of samples required by Mgc.In other words, a 2.0 indicates that the method requires double the sample size to achieve 85% power relative to Mgc.Pearson, RV, and Cca all achieve the same performance, as do Spearman and Kendall.Mgc requires the fewest number of samples in all settings, and for high-dimensional non-monotonic relationships, all other methods require about double or triple the number of samples Mgc requires.

Table 2 :
The p-values for brain imaging vs mental properties.Mgc always uncovers the existence of significant relationships and discovers the underlying optimal scales.Bold indicates significant p-value per dataset.
[45]t, we analyzed the relationship between resting-state functional magnetic resonance (rs-fMRI) activity and personality[45](see Appendix E.I for details).The first row of Table2compares the p-value of different methods, and Figure4Ashows the Mgc-Map for the sample data.Mgc is able to yield a significant p-value (< 0.05), whereas all previously proposed global dependence tests under consideration (Mantel, Dcorr, Mcorr, or Hhg) fail to detect dependence at a significance level of 0.05.Moreover, the The same power plots as in Figure2, except the 20 dependencies are one-dimensional with noise, and the x-axis shows sample size increasing from 5 to 100.Again, Mgc empirically achieves similar or better power than the previous state-of-the-art approaches on most problems.Note that Mic is included in 1D case; RV and Cca both equal Pearson in 1D; Kendall and Spearman are too similar to Pearson in power and thus omitted in plotting.
Figure E3: The same set-ups as in Figure 2, comparing different Mgc implementations versus its global counterparts.The default Mgc builds upon Mcorr throughout the paper, and we further consider Mgc on Mantel to illustrate the generalization.The magenta line shows the power difference between Mcorr and Mgc , and the cyan line shows the power difference between Mantel and the Mgc version of Mantel.Indeed, Mgc is able to improve the global counterpart in testing power under nonlinear dependencies, and maintains similar power under linear and independent dependencies.

Table 4 :
The Actual Testing Time (in seconds) on Real Data.