Discovering and deciphering relationships across disparate data modalities
Abstract
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 knearest 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 highdimensional 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.
https://doi.org/10.7554/eLife.41690.001eLife digest
If you want to estimate whether height is related to weight in humans, what would you do? You could measure the height and weight of a large number of people, and then run a statistical test. Such ‘independence tests’ can be thought of as a screening procedure: if the two properties (height and weight) are not related, then there is no point in proceeding with further analyses.
In the last 100 years different independence tests have been developed. However, classical approaches often fail to accurately discern relationships in the large, complex datasets typical of modern biomedical research. For example, connectomics datasets include tens or hundreds of thousands of connections between neurons that collectively underlie how the brain performs certain tasks. Discovering and deciphering relationships from these data is currently the largest barrier to progress in these fields. Another drawback to currently used methods of independence testing is that they act as a ‘black box’, giving an answer without making it clear how it was calculated. This can make it difficult for researchers to reproduce their findings – a key part of confirming a scientific discovery. Vogelstein et al. therefore sought to develop a method of performing independence tests on large datasets that can easily be both applied and interpreted by practicing scientists.
The method developed by Vogelstein et al., called Multiscale Graph Correlation (MGC, pronounced ‘magic’), combines recent developments in hypothesis testing, machine learning, and data science. The result is that MGC typically requires between one half to one third as big a sample size as previously proposed methods for analyzing large, complex datasets. Moreover, MGC also indicates the nature of the relationship between different properties; for example, whether it is a linear relationship or not.
Testing MGC on real biological data, including a cancer dataset and a human brain imaging dataset, revealed that it is more effective at finding possible relationships than other commonly used independence methods. MGC was also the only method that explained how it found those relationships.
MGC will enable relationships to be found in data across many fields of inquiry – and not only in biology. Scientists, policy analysts, data journalists, and corporate data scientists could all use MGC to learn about the relationships present in their data. To that extent, Vogelstein et al. have made the code open source in MATLAB, R, and Python.
https://doi.org/10.7554/eLife.41690.002Introduction
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 highthroughput screening (Zhang et al., 1999), precision medicine (Prescott, 2013), machine learning (Hastie et al., 2001), and causal analyses (Pearl, 2000). One of the first approaches for determining whether two properties are related to—or statistically dependent on—each other is Pearson’s ProductMoment Correlation (published in 1895; Pearson, 1895). This seminal paper prompted the development of entirely new ways of thinking about and quantifying relationships (see Reimherr and Nicolae, 2013 and Josse and Holmes, 2013 for recent reviews and discussion). Modern datasets, however, present challenges for dependencetesting 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 highdimensional 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 nonparametric tests to address nonlinear but univariate relationships (Hoeffding, 1948; Rényi, 1959). In the 1970s and 1980s, nearest neighbor style approaches were popularized (Friedman and Rafsky, 1983; Schilling, 1986), 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, at arbitrary dimensions, and structured data under a proper distance metric (Székely et al., 2007; Székely and Rizzo, 2009; Szekely and Rizzo, 2013; Lyons, 2013). Another set of methods, referred to a ‘kernel mean embedding’ approaches, including the Hilbert Schmidt Independence Criterion (Hsic) (Gretton and Gyorfi, 2010; Muandet et al., 2017), have the same theoretical guarantees, which is shown to be a kernel version of the energy statistics (Sejdinovic et al., 2013; Shen and Vogelstein, 2018). The energy statistics can perform very well with a relatively small sample size on highdimensional linear data, whereas the kernel methods and another test (Heller, Heller, and Gorfine’s test, Hhg) (Heller et al., 2013) perform well on lowdimensional nonlinear data. But no test performs particularly well on highdimensional 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, existing dependency tests do not attempt to further characterize the dependency structure. 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 (Bracewell and Bracewell, 1986) and wavelet analysis (Daubechies, 1992). More recently, topological and geometric data analysis compute properties of graphs, or even higher order simplices (Edelsbrunner and Harer, 2009). Such methods build multiscale characterization of the samples, much like recent developments in harmonic analysis (Coifman and Maggioni, 2006; Allard et al., 2012). However, these tools typically lack statistical guarantees under noisy observations and are often computationally burdensome.
We surmised that both (i) empirical performance in all dependency structures, in particular highdimensional, nonlinear, lowsample size settings, and (ii) providing insight into the discovery process, can be addressed via extending existing dependence tests to be adaptive to the data (Zhang et al., 2012). Existing tests rely on a fixed a priori selection of an algorithmic parameter, such as the kernel bandwidth (Gretton et al., 2006), intrinsic dimension (Allard et al., 2012), and/or local scale (Friedman and Rafsky, 1983; Schilling, 1986). Indeed, the Achilles Heel of manifold learning has been the requirement to manually choose these parameters (Levina and Bickel, 2004). Posthoc crossvalidation is often used to make these methods effectively 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 (Lee and Verleysen, 2007), 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 mechanism with distance 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. We compute all distance pairs, then efficiently compute the distance correlations for all scales. The local correlations (Figure 1C, described in detail below) illustrate which 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, and the geometric implication for the most informative scales. 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 comparisonbased approaches by adaptively estimating the informative scales for any relationship — linear or nonlinear, lowdimensional or highdimensional, unstructured or structured—in a computationally efficient and statistically valid and consistent fashion. This adaptive nature of Mgc effectively guarantees an improved statistical performance. Moreover, the dependency strength across all scales is informative about the structure of a statistical relationship, therefore providing further guidance for subsequent experimental or analytical steps. Mgc is thus a hypothesistesting and insightproviding approach that builds on recent developments in manifold and kernel learning, with complementary developments in nearestneighbor search, and multiscale analyses.
The multiscale graph correlation procedure
Mgc is a multistep procedure to discover and decipher dependencies across disparate data modalities or properties. Given $n$ samples of two different properties, proceed as follows (see Materials and methods and (Shen et al., 2018) for details):
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 distance matrix (by subtracting its overall mean, the columnwise mean from each column, and the rowwise mean from each row), and denote the resulting nbyn matrices $A$ and $B$.
For all possible values of $k$ and $l$ from $1$ to $n$:
Compute the $k$nearest neighbor graphs for one property, and the $l$nearest neighbor graphs for the other property. Let ${G}_{k}$ and ${H}_{l}$ be the adjacency matrices for the nearest neighbor graphs, so that ${G}_{k}(i,j)=1$ indicates that $A(i,j)$ is within the $k$ smallest values of the $i}^{th$ row of $A$, and similarly for ${H}_{l}$.
Estimate the local correlations—the correlation between distances restricted to only the $(k,l)$ neighbors—by summing the products of the above matrices, ${c}^{kl}=\sum _{ij}A(i,j){G}_{k}(i,j)B(i,j){H}_{l}(i,j)$.
Normalize $c}^{kl$ such that the result is always between $1$ and $+1$ by dividing by$\sqrt{\sum _{ij}{A}^{2}(i,j){G}_{k}(i,j)\times \sum _{ij}{B}^{2}(i,j){H}_{l}(i,j)}.$
Estimate the optimal local correlation ${c}^{*}$ by finding the smoothed maximum of all local correlations $\{{c}^{kl}\}$. Smoothing mitigates biases and provides Mgc with theoretical guarantee and better finitesample performance.
Determine whether the relationship is significantly dependent—that is, whether ${c}^{*}$ 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 pvalue, rather than one pvalue 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).
Computing all local correlations, the test statistic, and the pvalue requires $O({n}^{2}\mathrm{log}n)$ time, which is about the same running time complexity as other methods (Shen et al., 2018).
Results
Mgc typically requires substantially fewer samples to achieve the same power across all dependencies and dimensions
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 Materials and methods for details):
The standard criterion for evaluating statistical tests is the testing power, which equals the probability that a test correctly rejects the null hypothesis at a given type one error level, that is power = Prob(${H}_{0}$ is rejected $\phantom{\rule{thinmathspace}{0ex}}{H}_{0}$ is false). The higher the testing power, the better the test procedure. A consistent test has power converging to $1$ under dependence, and a valid test controls the type one error level under independence. In a complementary manuscript (Shen et al., 2018), we established the theoretical properties of Mgc, proving its validity and universal consistency for dependence testing against all distributions of finite second moments.
Here, we address the empirical performance of Mgc as compared with multiple popular tests: (i) Dcorr, a popular approach from the statistics community (Székely et al., 2007; Székely and Rizzo, 2009), (ii) Mcorr, a modified version of Dcorr designed to be unbiased for sample data (Szekely and Rizzo, 2013), (iii) Hhg, a distancebased test that is very powerful for detecting lowdimensional nonlinear relationships (Heller et al., 2013). (iv) Hsic, a kernel dependency measure (Gretton and Gyorfi, 2010) formulated in the same way as Dcorr except operating on kernels, (v) Mantel, which is historically widely used in biology and ecology (Mantel, 1967). (vi) RV coefficient (Pearson, 1895; Josse and Holmes, 2013), which is a multivariate generalization of Pearson’s product moment correlation whose test statistic is the sum of the tracenorm of the crosscovariance matrix, and (vii) the Cca method, which is the largest (in magnitude) singular value of the crosscovariance matrix, and can be viewed as a different generalization of Pearson in highdimensions that is more appropriate for sparse settings (Hotelling, 1936; Witten et al., 2009; Witten and Tibshirani, 2011). Note that while we focus on highdimensional settings, Appendix 1 shows further results in onedimensional 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 (Spearman, 1904), (x) Kendall’s tau correlation (Kendall, 1970), and (xi) Mic (Reshef et al., 2011). Under the regularity condition that the data distribution has finite second moment, the first four tests are universally consistent, whereas the other tests are not.
We generate an extensive benchmark suite of 20 relationships, including different polynomial (linear, quadratic, cubic), trigonometric (sinusoidal, circular, ellipsoidal, spiral), geometric (square, diamond, Wshape), and other functions. This suite includes and extends the simulated settings from previous dependence testing work (Székely et al., 2007; Simon and Tibshirani, 2012; Gorfine et al., 2012; Heller et al., 2013; Szekely and Rizzo, 2013). For many of them, we introduce highdimensional variants, to more extensively evaluate the methods; function details are in Materials and methods. The visualization of onedimensional noisefree (black) and noisy (gray) samples is shown in Figure 2—figure supplement 1. 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 highdimensional 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. Figure 2—figure supplement 2 shows the same advantage in onedimension with increasing sample size.
Moreover, for each relationship and each method we compute the required sample size to achieve power 85% at error level 0.05, and summarize the median size for monotone relationships (type 1–5) and nonmonotone relationships (type 6–19) in Table 1. Other methods typically require double or triple the 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, distancebased 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 nonmonotonic 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 highdimensional nonmonotonic relationships that motivated this work, and are common in biomedicine, Mgc significantly outperforms other methods. 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 at 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 (Szekely and Rizzo, 2013). Figure 2—figure supplement 3 directly compares multiscale generalizations of Mantel and Mcorr as dimension increases, demonstrating that empirically, Mgc nearly dominates its global variant for essen tially all dimensions and simulation settings considered here. Figure 2—figure supplement 4 shows a similar result for onedimensional 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.

Table 1—source data 1
 https://doi.org/10.7554/eLife.41690.010

Table 1—source data 2
 https://doi.org/10.7554/eLife.41690.011
Mgc deciphers latent dependence structure
Beyond simply testing the existence of a relationship, the next goal is often to decipher the nature or structure of the 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 (Berk et al., 2013). 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 MgcMap shows local correlation as a function of the scales of the two properties. More concretely, it is the matrix of $c}^{kl$’s, as defined above. Thus, the MgcMap is an nbyn matrix which encodes the strength of dependence for each possible scale. Figure 3 provides the MgcMap for all 20 different onedimensional relationships; the optimal scale to achieve ${\widehat{t}}_{*}$ is marked with a green dot. For the monotonic dependencies (15), the optimal scale is always the largest scale, that is the global one. For all nonmonotonic dependencies (619), Mgc chooses smaller scales. Thus, a global optimal scale implies a closetolinear dependency, otherwise the dependency is strongly nonlinear. In fact, this empirical observation led to the following theorem (which is proved in Materials and methods):
Theorem 1. When $(X,Y)$ are linearly related (meaning that $Y$ can be constructed from $X$ by rotation, scaling, translation, and/or reflection), the optimal scale of Mgc equals the global scale. Conversely, a local optimal scale implies a nonlinear relationship.
Thus, the MgcMap explains how Mgc discovers relationships, specifically, which scale has the most informative pairwise comparisons, and how that relates to the geometry 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 MgcMaps and often similar optimal scales. For example, logarithmic (10) and fourth root (11), although very different functions analytically, are geometrically similar, and yield very similar MgcMaps. Similarly, (12) and (13) are trigonometric functions, and they share a narrow range of significant local scales. Both circle (16) and ellipse (17), as well as square (14) and diamond (18), are closely related geometrically and also have similar MgcMaps. This indicates that the MgcMap partially characterizes the geometry of these relationships, differentiating different dependence structures and assisting subsequent analysis steps. Moreover, in Shen and Vogelstein, 2018, we proved that the sample MgcMap (which Mgc estimates) converges to the true MgcMap 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 multiscale geometric or topological characterizations of data, such as persistence diagrams (Edelsbrunner and Harer, 2009).
Mgc is computationally efficient
Mgc does not incur large computational costs and has a similar complexity as existing methods. Though a naïve implementation of Mgc requires $\mathcal{\mathcal{O}}({n}^{4})$ operations, we devised a nested implementation that requires only $\mathcal{\mathcal{O}}({n}^{2}\mathrm{log}n)$ operations. Moreover, obtaining the MgcMap costs no additional computation, whereas other methods would require running a secondary computational step to decipher geometric properties of the relationship. Mgc can also trivially be parallelized, reducing computation to $\mathcal{\mathcal{O}}({n}^{2}\mathrm{log}n/T)$, where $T$ is the number of cores (see Algorithm C1 for details). Since $T$ is often larger than $\mathrm{log}n$, in practice, Mgc can be $\mathcal{\mathcal{O}}({n}^{2})$, meaning only a constant factor slower than Dcorr and Hsic, which is illustrated in Figure 6 of Shen and Vogelstein, 2018. For example, at sample size $n=5000$ and dimension $p=1$, on a typical laptop computer, Dcorr requires around 0.5 s to compute the test statistic, whereas Mgc requires no more than $5$ s. But the cost and time to obtain 2.5× more data (so Dcorr has same average power as Mgc) typically far exceeds a few seconds. In comparison, the cost to compute a persistence diagram is typically $\mathcal{\mathcal{O}}({n}^{3})$, which is orders of magnitude slower when $n\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}10$. The running time of each method on the real data experiments are reported in Materials and methods.
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 (Costa and McCrae, 1992; Jung et al., 2009). However, the relationship between brain activity and structure, and these aspects of our psyche, remains unclear (DeYoung et al., 2010; Xu and Potenza, 2012; Bjørnebekk et al., 2013; Sampaio et al., 2014). For example, prior work did not evaluate the relationship between entire brain connectivity and all five factors of the standard personality model (Costa and McCrae, 1992). We therefore utilized Mgc to investigate published open access data (see Materials and methods for details).
First, we analyzed the relationship between restingstate functional magnetic resonance (rsfMRI) activity and personality (Adelstein et al., 2011). The first row of Table 2 compares the pvalue of different methods, and Figure 4A shows the MgcMap for the sample data. Mgc is able to yield a significant pvalue (< 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 MgcMap provides a characterization of the dependence, for which the optimal scale indicates that the dependency is strongly nonlinear. Interestingly, the MgcMap 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.

Table 2—source data 1
 https://doi.org/10.7554/eLife.41690.016

Table 2—source data 2
 https://doi.org/10.7554/eLife.41690.017
Second, we investigated the relationship between diffusion MRI derived connectivity and creativity (Jung et al., 2009). 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 MgcMap in Figure 4B closely resembles a linear relationship where the optimal scale is global. The closetolinear relationship is also supported from the pvalue table as all methods except Hsic are able to detect significant dependency, which suggests that there is relatively little to gain by pursuing nonlinear regression techniques, potentially saving valuable research time by avoiding tackling an unnecessary problem. The test statistic for both Mgc and Mcorr equal 0.04, which is quite close to zero despite a significant pvalue, implying a relatively weak and noisy relationship. A prediction of creativity via linear regression turns out to be nonsignificant, which implies that the sample size is too low to obtain useful predictive accuracy (not shown), indicating that more data are required for single subject predictions. If one had first directly estimated the regression function, obtaining a null result, it would remain unclear whether a relationship existed. This experiment demonstrates that for highdimensional and potentially structured data, Mgc is able to reveal dependency with relatively small sample size while parametric techniques and directly estimating regression functions can often be ineffective.
The performance in the real data closely matches the simulations results: the first dataset exhibits a strongly nonlinear relationship, for which Mgc has the lowest pvalue, followed by Hhg and Hsic and then all other methods; the second dataset exhibits a closetolinear relationship, for which global methods 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 but noisy 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 Eklund et al. (2012); Eklund et al. (2016). Appendix 1—figure 1 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 (Frantzi et al., 2014). We sought to discover new useful protein biomarkers from a quantitative proteomics technique that measures protein and peptide abundance called Selected Reaction Monitoring (SRM) (Wang et al., 2011). 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 (Bhat et al., 2012).
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 (Wang et al., 2017). Therefore, our first step was to check whether Mgc could correctly identify ovarian biomarkers. Indeed, the peptides that have been validated previously are also identified by Mgc. 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 pvalue 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 (Yang et al., 2015; Willemse et al., 2018). 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 (Helfman et al., 2018; Lam et al., 2012). The rest of the global methods did not identify any markers at significance level $0.05$, see Materials and methods for more details and Appendix 1—table 2 for identified peptide information using each method.
Since there is no ground truth yet in this experiment, we further carried out a classification task using the biomarkers identified by the various algorithms, using a knearestneighbor classifier to predict pancreatic cancer, and a leaveonesubjectout 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 highdimensional and illposed problems, where dimensionality is larger than sample size. Second, Mgc can also be thought of as learning a metric because it chooses the optimal scale amongst a set of ${n}^{2}$ truncated distances, motivating studying the relationship between Mgc and recent advances in metric learning (Xing et al., 2003). In particular, deep learning can be thought of as metric learning (Giryes et al., 2015), and generative adversarial networks (Goodfellow et al., 2014) are implicitly testing for equality, which is closely related to dependence (Sutherland et al., 2016). While Mgc searches over a twodimensional 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 goodnessoffit (Székely and Rizzo, 2005), analysis of variance (Rizzo and Székely, 2010), conditional dependence (Székely and Rizzo, 2014; Wang et al., 2015), and feature selection (Li et al., 2012; Zhong and Zhu, 2015), so Mgc can be adapted for them as well. Indeed, Mgc can also implement a twosample (or generally the Ksample) test (Szekely and Rizzo, 2004; Heller et al., 2016; Shen and Vogelstein, 2018). Specifically, for more than two modalities, one may use summation of pairwise Mgc test statistics, similar to how energy statistic is generalized to Ksample testing from twosample testing (Rizzo and Székely, 2010; Rizzo and Székely, 2016; Shen and Vogelstein, 2018), or how canonical correlation analysis is generalized into more than two modalities (Kettenring, 1971; Tenenhaus and Tenenhaus, 2011; Shen et al., 2014). Finally, although energy statistics have not yet been explicitly 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 a strong effect on the performance of different machine learning algorithms (Levina and Bickel, 2004). 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 (Tenenbaum et al., 2000) and local linear embedding (Roweis and Saul, 2000) 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 default metric choice of Mgc in this paper is always the Euclidean distance, but other metric choices may be more appropriate in different fields, and using the strong negative type metric as specified in Lyons (2013) guarantees consistency. However, if multiple metric choices are experimented to yield multiple Mgc pvalues, then the optimal pvalue should be properly corrected for multiple testing. Alternatively, one may use the maximum Mgc statistic among multiple metric choices, apply the same procedure in each permutation (i.e. in each permutation, use the same number of metric choices and take the maximum Mgc as the permuted statistic), then derive a single pvalue. Such a testing procedure properly controls the type one error level without the need for additional correction.
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 ‘postselection inference’ problem (Berk et al., 2013). The most straightforward way to address this issue is to collect new data, which is costly and timeconsuming. 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 pvalue, rather than after. To our knowledge, Mgc is the first dependence test to take a step towards valid postselection 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 demonstrated that one could reduce computational time of distancebased tests to close to linear via faster implementation, subsampling, random projection, and null distribution approximation (Huo and Székely, 2016; Huang and Huo, 2017; Zhang et al., 2018; Chaudhuri and Hu, 2018), making it feasible for large amount of data. Alternately, semiexternal memory implementations would allow running Mgc even as the interpoint comparison matrix exceeds the size of main memory (Da Zheng et al., 2015; Da Zheng et al., 2016a; Da Zheng et al., 2016b; Da Zheng et al., 2016c).
Finally, Mgc is easy to use. Source code is available in MATLAB, R, and Python from https://mgc.neurodata.io/ (Bridgeford et al., 2018; experiments archived at https://github.com/elifesciencespublications/MGCpaper). Code for reproducing all the figures in this manuscript is also available from the above websites. We showed Mgc’s 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 Mgc.
Materials and methods
Mathematical details
This section contains essential mathematical details on independence testing, the notion of the generalized correlation coefficient and the distancebased correlation measure, how to compute the local correlations, and the smoothing technique. A statistical treatment on MGC is in Shen and Vogelstein, 2018, which introduces the population version of Mgc and various theoretical properties.
Testing independence
Request a detailed protocolGiven pairs of observations $({\mathit{x}}_{i},{\mathit{y}}_{i})\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{q}$ for $i=1,\mathrm{\dots},n$, assume they are independently identically distributed as $(X,Y)\stackrel{iid}{\sim}{F}_{XY}$. If the two random variables $X$ and $Y$ are independent, the joint distribution equals the product of the marginals, that is $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, that is 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 $\mathit{x}$’s, $q$ is the dimension for $\mathit{y}$’s. For Mgc and all benchmark methods, there is no restriction on the dimensions, that is 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 onedimensional data, such as (Reshef et al., 2011; Heller et al., 2016; Huo and Székely, 2016), which are not directly applicable to multidimensional data.
Correlation measures
Request a detailed protocolTo achieve consistent testing, most stateoftheart dependence measures operate on pairwise comparisons, either similarities (such as kernels) or dissimilarities (such as distances).
Let $\mathcal{\mathcal{X}}}_{n}=\{{\mathit{x}}_{1},\cdots ,{\mathit{x}}_{n}\}\in {\mathbb{R}}^{p\times n$ and $\mathcal{\mathcal{Y}}}_{n}=\{{\mathit{y}}_{1},\cdots ,{\mathit{y}}_{n}\}\in {\mathbb{R}}^{q\times n$ denote the matrices of sample observations, and ${d}_{x}$ be the distance function for $\mathit{x}$’s and ${d}_{y}$ for $\mathit{y}$’s. One can then compute two $n\times n$ distance matrices $\stackrel{~}{A}=\{{\stackrel{~}{a}}_{ij}\}$ and $\stackrel{~}{B}=\{{\stackrel{~}{b}}_{ij}\}$, where ${\stackrel{~}{a}}_{ij}={\delta}_{x}({\mathit{x}}_{i},{\mathit{x}}_{j})$ and ${\stackrel{~}{b}}_{ij}={\delta}_{y}({\mathit{y}}_{i},{\mathit{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.
Let $A$ and $B$ be the transformed (e.g., centered) versions of the distance matrices $\stackrel{~}{A}$ and $\stackrel{~}{B}$, respectively. Any ‘generalized correlation coefficient’ (Spearman, 1904; Kendall, 1970) can be written as:
where $z$ is proportional to the standard deviations of $A$ and $B$, that is $z={n}^{2}{\sigma}_{a}{\sigma}_{b}$. In words, $c$ is the global sample correlation across pairwise comparison matrices $A$ and $B$, and is normalized into the range $[1,1]$, which usually 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 via the above correlation formulation, by using $A$ and $B$ directly from sample observations rather than distances. Distancebased methods like Dcorr and Mantel operate on the Euclidean distance by default, or other metric choices on the basis of domain knowledge; then transform the resulting distance matrices $\stackrel{~}{A}$ and $\stackrel{~}{B}$ by certain centering schemes into $A$ and $B$. Hsic chooses the Gaussian kernel and computes two kernel matrices, then transform the kernel matrices $\stackrel{~}{A}$ and $\stackrel{~}{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 Sejdinovic et al. (2013)), 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, for example a generalized correlation, the permutation test is often an effective choice (Good, 2005), because a pvalue 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 pvalue is lower than a predetermined 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 autocorrelated time series (Guillot and Rousset, 2013).
Distance correlation (Dcorr) and the Unbiased Version (Mcorr)
Request a detailed protocolDefine the row and column means of $\stackrel{~}{A}$ by $\overline{a}}_{\cdot j}=\frac{1}{n}\sum _{i=1}^{n}{\stackrel{~}{a}}_{ij$ and $\overline{a}}_{i\cdot}=\frac{1}{n}\sum _{j=1}^{n}{\stackrel{~}{a}}_{ij$. Dcorr defines
and similarly for $b}_{ij$. For distance correlation, the numerator of Equation 1 is named the distance covariance (Dcov), while ${s}_{a}$ and ${s}_{b}$ in the denominator are the square root of each distance variance. The centering scheme is important to guarantee the universal consistency of Dcorr, whereas Mantel uses a simple centering scheme and thus not universally consistent.
Let $c(X,Y)$ be the population distance correlation, that is, the distance correlation between the underlying random variables $X$ and $Y$. Székely 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\to \mathrm{\infty}$, the sample distance correlation converges to the population distance correlation, that is, $c({\mathcal{\mathcal{X}}}_{n},{\mathcal{\mathcal{Y}}}_{n})\to 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 nonnegative. Moreover, the consistency result holds for a much larger family of metrics, those of strong negative type (Lyons, 2013).
It turns out that the sample distance correlation has a finitesample bias, especially as the dimension $p$ or $q$ increases (Szekely and Rizzo, 2013). For example, for independent Gaussian distributions, the sample distance correlation converges to $1$ as $p,q\to \mathrm{\infty}$. By excluding the diagonal entries and slightly modifies the offdiagonal entries of $\mathcal{\mathcal{A}}$ and $B$, Szekely and Rizzo (Szekely and Rizzo, 2013; Székely and Rizzo, 2014) show that Mcorr is an unbiased estimator of the population distance correlation $c(\mathit{x},\mathit{y})$ for all $p,q,n$, which is approximately normal even if $p,q\to \mathrm{\infty}$. Thus it enjoys the same theoretical consistency as Dcorr and always has zero mean under independence.
Local correlations
Request a detailed protocolGiven any matrices $A$ and $B$, we can define a set of local correlations as follows. Let $R({A}_{\cdot j},i)$ be the ‘rank’ of $\mathit{x}}_{i$ relative to $\mathit{x}}_{j$, that is, $R({A}_{\cdot j},i)=k$ if $\mathit{x}}_{i$ is the $k}^{th$ closest point (or ‘neighbor’) to $\mathit{x}}_{j$, as determined by ranking the $n1$ distances to ${x}_{j}$. Define $R({B}_{i\cdot},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 $\mathit{x}}_{i$ and any neighborhood size $l$ around each $\mathit{y}}_{j$, we define the local pairwise comparisons:
and then let $a}_{ij}^{k}={\stackrel{~}{a}}_{ij}^{k}{\overline{a}}^{k$, where ${\overline{a}}^{k}$ is the mean of $\{{\stackrel{~}{a}}_{ij}^{k}\}$, and similarly for $b}_{ij}^{l$.
The local correlation coefficient at a given scale is defined to effectively exclude large distances:
where $z}_{kl}={n}^{2}{\sigma}_{a}^{k}{\sigma}_{b}^{l$, with ${s}_{a}^{k}$ and ${s}_{b}^{l}$ is the standard deviations for the truncated pairwise comparisons. The MgcMap can be constructed by computing all local correlations, which allows the discovery of the optimal correlation. For any aforementioned correlation (Dcorr, Mcorr, Hsic, Mantel, Pearson), one can define its local correlations by using Equation 3 and plugging in the respective $a}_{ij$ and $b}_{ij$ from Equation 1.
As most nonlinear relationships intrinsically exhibit a local linear structure, considering the nearestneighbors is able to amplify the dependency signal over the global correlation. There could be two other scenarios: when the small distances in one modality mostly correspond to large distances in another modality, or when the large distances in one modality correspond to large distance in another modality. For the first scenario, the small distances become negative terms after centering while the large distances become positive terms after centering, so adding their product to $c}^{kl$ will cause the test statistic to be smaller — in fact, as distance correlation is shown to be > 0 under dependence (Székely et al., 2007), the first scenario cannot happen for all distances pairs. For the second scenario, one can experiment using the large distances (or the furthest neighbors) only by reversing the ranking scheme in local correlation to descending order. However, whenever the large distances are highly correlated, the small distances must also be highly correlated after centering by the mean distances, so global correlation coefficient like Dcorr already handles this scenario. Therefore considering the nearestneighbor may significantly improve the performance over global correlation, while considering the other scenarios does not.
Mgc as the optimal local correlation
Request a detailed protocolWe define the multiscale graph correlation statistic as the optimal local correlation, for which the family of local correlation is computed based on Euclidean distance and Mcorr transformation.
Instead of taking a direct maximum, Mgc takes a smoothed maximum, that is 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 C2). Thus, we can write:
Then the optimal scale equals all scales within $R$ whose local correlations are as large as $c}^{\ast$. The choice of $\tau$ is made explicit in the pseudocode, with further discussion and justification offered in Shen and Vogelstein, 2018.
Proof for theorem 1
Request a detailed protocolTheorem 1. When $\mathrm{(}X\mathrm{,}Y\mathrm{)}$ 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=WX+b$ for a unitary matrix $W$ and a constant $b$ upto possible scaling, in which case the distances are preserved, that is $\Vert {y}_{i}{y}_{j}\Vert =\Vert W{x}_{i}W{x}_{j}\Vert =\Vert {x}_{i}{x}_{j}\Vert$. It follows that $\mathtt{\text{M}}\mathtt{\text{corr}}({\mathcal{\mathcal{X}}}_{n},{\mathcal{\mathcal{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 $\mathtt{\text{M}}\mathtt{\text{gc}}({\mathcal{\mathcal{X}}}_{n},{\mathcal{\mathcal{Y}}}_{n})=\mathtt{\text{M}}\mathtt{\text{corr}}({\mathcal{\mathcal{X}}}_{n},{\mathcal{\mathcal{Y}}}_{n})=1$.
Computational complexity of each step
Request a detailed protocolThe distance computation takes $\mathcal{\mathcal{O}}({n}^{2}max\{p,q\})$, and the ranking process takes $\mathcal{\mathcal{O}}({n}^{2}\mathrm{log}n)$. Once the distance and ranking are completed, computing one local generalized correlation requires $\mathcal{\mathcal{O}}({n}^{2})$ (see Algorithm C4). Thus, a naive approach to compute all local generalized correlations requires at least $\mathcal{\mathcal{O}}({n}^{2}max\{{n}^{2},p,q\})$ by going through all possible scales, meaning possibly $\mathcal{\mathcal{O}}({n}^{4})$ which would be computationally prohibitive. However, given the distance and ranking information, we devised an algorithm that iteratively computes all local correlations in $\mathcal{\mathcal{O}}({n}^{2})$ by reusing adjacent smaller local generalized correlations (see Algorithm C5). Therefore, when including the distance computation and ranking overheads, the MGC statistic is computed in $\mathcal{\mathcal{O}}({n}^{2}max\{\mathrm{log}n,p,q\})$), which has the same running time as the Hhg statistic, and the same running time up to a factor of $\mathrm{log}n$ as global correlations like Dcorr and Mcorr, which require $\mathcal{\mathcal{O}}({n}^{2}max\{p,q\})$ time. By utilizing a multicore architecture, Mgc can be computed in $\mathcal{\mathcal{O}}({n}^{2}max\{\mathrm{log}n,p,q\}/T)$ instead. As $T=\mathrm{log}(n)$ is often a small number, for example $T$ is no more than $30$ at $1$ billion samples, thus Mgc can be effectively computed in the same complexity as Dcorr. Note that the permutation test adds another $r$ random permutations to the ${n}^{2}$ term, so computing the pvalue requires $\mathcal{\mathcal{O}}({n}^{2}max\{\mathrm{log}n,p,q,r\}/T)$.
Mgc algorithms and testing procedures
Request a detailed protocolSix algorithms are presented in order:
Algorithm C1 describes Mgc in its entirety (which calls most of the other algorithms as functions).
Algorithm C2 computes the Mgc test statistic.
Algorithm C3 computes the pvalue of Mgc by the permutation test.
Algorithm C4 computes the local generalized correlation coefficient at a given scale $(k,l)$, for a given choice of the global correlation coefficient.
Algorithm C5 efficiently computes all local generalized correlations, in nearly the same running time complexity as computing one local generalized correlation.
Algorithm C6 evaluates the testing power of Mgc by a given distribution.
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 $\mathcal{\mathcal{O}}({n}^{2}max(\mathrm{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 pvalue $p({c}^{\ast})$,  
function MG$(({x}_{i},{y}_{i})$, for $i\in [n])$  
(1) Calculate all pairwise distances:  
for $i,j:=1,\mathrm{\dots},n$ do  
${a}_{ij}={\delta}_{x}({x}_{i},{x}_{j})$  ${d}_{x}$ is the distance between pairs of $x$ samples  
${b}_{ij}={\delta}_{y}({y}_{i},{y}_{j})$  ${d}_{y}$ is the distance between pairs of $y$ samples  
end for  
Let $A=\{{a}_{ij}\}$ and $B=\{{b}_{ij}\}$.  
(2) Calculate Multiscale Correlation Map $\mathcal{\mathcal{C}}$ & Mgc Test Statistic:  
$[{c}^{\ast},\mathcal{\mathcal{C}},k,l]=\mathrm{M}\mathrm{G}\mathrm{C}\mathrm{S}\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{L}\mathrm{E}\mathrm{S}\mathrm{T}\mathrm{A}\mathrm{T}(A,B)$  Algorithm C2  
(3) Calculate the pvalue  
$pval({c}^{\ast})=\mathrm{P}\mathrm{E}\mathrm{R}\mathrm{M}\mathrm{U}\mathrm{T}\mathrm{A}\mathrm{T}\mathrm{I}\mathrm{O}\mathrm{N}\mathrm{T}\mathrm{E}\mathrm{S}\mathrm{T}(A,B,r,{c}^{\ast})$  Algorithm C3  
end Function 
Pseudocode C2 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, that is 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 $\mathcal{\mathcal{O}}({n}^{2})$.  
Input: A pair of distance matrices $(A,B)\in {\mathbb{R}}^{n\times n}\times {\mathbb{R}}^{n\times n}$.  
Output: The Mgc statistic $c}^{\ast}\in \mathbb{R$, all local statistics $\mathcal{\mathcal{C}}\in {\mathbb{R}}^{n\times n}$, and the corresponding local scale $(k,l)\in \mathbb{N}\times \mathbb{N}$.  
1:  function MGCSampleStat$(A,B)$  
2:  $\mathcal{\mathcal{C}}=\mathrm{M}\mathrm{G}\mathrm{C}\mathrm{A}\mathrm{L}\mathrm{L}\mathrm{L}\mathrm{O}\mathrm{C}\mathrm{A}\mathrm{L}(A,B)$  All local correlations 
3:  $\tau =\mathrm{T}\mathrm{H}\mathrm{R}\mathrm{E}\mathrm{S}\mathrm{H}\mathrm{O}\mathrm{L}\mathrm{D}\mathrm{I}\mathrm{N}\mathrm{G}(\mathcal{\mathcal{C}})$  find a threshold to determine large local correlations 
4:  for $i,j:=1,\dots ,n$do ${r}_{ij}\leftarrow \mathbb{I}({c}^{ij}>\tau )$end for  identify all scales with large correlation 
5:  $\mathcal{\mathcal{R}}\leftarrow \{{r}_{ij}:i,j=1,\dots ,n\}$  binary map encoding scales with large correlation 
6:  $\mathcal{\mathcal{R}}=\mathrm{C}\mathrm{O}\mathrm{N}\mathrm{N}\mathrm{E}\mathrm{C}\mathrm{T}\mathrm{E}\mathrm{D}(\mathcal{\mathcal{R}})$  largest connected component of the binary matrix 
7:  ${c}^{\ast}\leftarrow \mathcal{\mathcal{C}}(n,n)$  use the global correlation by default 
8:  $k\leftarrow n,l\leftarrow n$  
9:  if $\left(\sum _{i,j}{r}_{ij}\right)\ge 2n$ then  proceed when the significant region is sufficiently large 
10:  $[{c}^{\ast},k,l]\leftarrow max(\mathcal{\mathcal{C}}\circ \mathcal{\mathcal{R}})$  find the smoothed maximum and the respective scale 
11:  end if  
12:  end Function  
Input: $\mathcal{\mathcal{C}}\in {\mathbb{R}}^{n\times n}$.  
Output: A threshold $t$ to identify large correlations.  
13:  function Thresholding $\mathcal{\mathcal{C}}$  
14:  $\tau \leftarrow {\sum}_{{c}^{ij}<\phantom{\rule{thinmathspace}{0ex}}0}({c}^{ij}{)}^{2}/{\sum}_{{c}^{ij}<\phantom{\rule{thinmathspace}{0ex}}0}1$  variance of all negative local generalized correlations 
15:  $\tau \leftarrow max\{0.01,\sqrt{\tau}\}\times 3.5$  threshold based on negative correlations 
16:  $\tau \leftarrow max\{\tau ,2/n,{c}^{nn}\}$  
17:  end Function 
Pseudocode C3 Permutation Test. This algorithm uses the random permutation test with $r$ random permutations for the pvalue, requiring $\mathcal{\mathcal{O}}(r{n}^{2}\mathrm{log}n)$ for Mgc. In the realdata experiment, we always set $r=10$,$000$. Note that the pvalue 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)\in {\mathbb{R}}^{n\times n}\times {\mathbb{R}}^{n\times n}$, the number of permutations $r$, and Mgc statistic ${c}^{*}$ for the observed data.  
Output: The pvalue $pval\in [0,1]$.  
1:  function PermutationTest($A$, $B$, $r$, ${c}^{*}$)  
2:  for $t:=1,\mathrm{\dots},r$ do  
3:  $\pi =\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{D}\mathrm{P}\mathrm{E}\mathrm{R}\mathrm{M}(n)$  generate a random permutation of size $n$ 
4:  $c}_{0}^{\ast}[t]=\mathrm{M}\mathrm{G}\mathrm{C}\mathrm{S}\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{L}\mathrm{E}\mathrm{S}\mathrm{T}\mathrm{A}\mathrm{T}(A,B(\pi ,\pi ))\phantom{\rule{0ex}{0ex}$  calculate the permuted Mgc statistic 
5:  end for  
6:  $pval({c}^{\ast})\leftarrow \frac{1}{t}\sum _{t=1}^{r}\mathit{I}({c}^{\ast}\le {c}_{0}^{\ast}[t])$  compute pvalue of Mgc 
7:  end function 
Pseudocode C4 Compute local test statistic at a given scale. This algorithm runs in $\mathcal{\mathcal{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 $\mathcal{\mathcal{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=\mathrm{log}(n)$ cores, the sorting function can be easily parallelized to take $\mathcal{\mathcal{O}}({n}^{2}\mathrm{log}(n)/T)=\mathcal{\mathcal{O}}({n}^{2})$.  
Input: A pair of distance matrices $(A,B)\in {\mathbb{R}}^{n\times n}\times {\mathbb{R}}^{n\times n}$, and a local scale $(k,l)\in \mathbb{N}\times \mathbb{N}$.  
Output: The local generalized correlation coefficient ${c}^{kl}\in [1,1]$.  
1:  function LocalGenCorr($A$, $B$, $k$, $l$)  
2:  for $Z:=A,B$ do ${\mathcal{\mathcal{E}}}^{Z}=\mathrm{S}\mathrm{O}\mathrm{R}\mathrm{T}(Z)$ end for  parallelized sorting 
3:  for $Z:=A,B$ do $Z=\mathrm{C}\mathrm{E}\mathrm{N}\mathrm{T}\mathrm{E}\mathrm{R}(Z)$ end for  center distance matrices 
4:  ${\stackrel{~}{c}}^{kl}\leftarrow tr((A\circ {\mathcal{\mathcal{E}}}^{A}{)}^{\mathrm{T}}\times (B\circ ({\mathcal{\mathcal{E}}}^{B}{)}^{\mathrm{T}}))$  unnormalized local distance covariance 
5:  ${v}^{A}\leftarrow tr((A\circ {\mathcal{\mathcal{E}}}^{A}{)}^{\mathrm{T}}\times (A\circ ({\mathcal{\mathcal{E}}}^{A}{)}^{\mathrm{T}}))$  local distance variances 
6:  ${v}^{B}\leftarrow tr((B\circ {\mathcal{\mathcal{E}}}^{B}{)}^{\mathrm{T}}\times (B\circ ({\mathcal{\mathcal{E}}}^{B}{)}^{\mathrm{T}}))$  
7:  $e}^{A}\leftarrow \sum _{i,j=1}^{n}(A\circ {\mathcal{\mathcal{E}}}^{A}{)}_{ij$  sample means 
8:  $e}^{B}\leftarrow \sum _{i,j=1}^{n}(B\circ {\mathcal{\mathcal{E}}}^{B}{)}_{ij$  
9:  $c}^{kl}\leftarrow \left({\stackrel{~}{c}}^{kl}{e}^{A}{e}^{B}/{n}^{2}\right)/\sqrt{\left({v}^{A}({e}^{A}/n{)}^{2}\right)\left({v}^{B}({e}^{B}/n{)}^{2}\right)$  center and normalize 
10:  end function 
Pseudocode C5 Compute the multiscale correlation map (i.e., all local generalized correlations) in $\mathcal{\mathcal{O}}({n}^{2}\mathrm{log}n/T)$. Once the distances are sorted, the remaining algorithm runs in $\mathcal{\mathcal{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\le R({A}_{\cdot j},i)$ and $l\le R({B}_{\cdot j},i)$, so it suffices to iterate through $a}_{ij}{b}_{ij$ for $i,j:=1,\mathrm{\dots},n$, and add the product simultaneously to all $c}^{kl$ whose scales are no more than $(R({A}_{\cdot j},i),R({B}_{\cdot j},i))$. To achieve the above, we iterate through each product, add it to $c}^{kl$ at $(kl)=(R({A}_{\cdot j},i),R({B}_{\cdot j},i))$ only (so only one local scale is accessed for each operation); then add up adjacent $c}^{kl$ for $k,l=1,\mathrm{\dots},n$. The same applies to all local covariances, variances, and expectations.  
Input: A pair of distance matrices $(A,B)\in {\mathbb{R}}^{n\times n}\times {\mathbb{R}}^{n\times n}$.  
Output: The multiscale correlation map $\mathcal{\mathcal{C}}\in [1,1{]}^{n\times n}$for $k,l=1,\mathrm{\dots},n$.  
1:  function MGCAllLocal($A$, $B$)  
2:  for $Z:=A,B$ do ${\mathcal{\mathcal{E}}}^{Z}=\mathrm{S}\mathrm{O}\mathrm{R}\mathrm{T}(Z)$ end for  
3:  for $Z:=A,B$ do $Z=\mathrm{C}\mathrm{E}\mathrm{N}\mathrm{T}\mathrm{E}\mathrm{R}(Z)$end for  
4:  for $i,j:=1,\dots ,n$ do  iterate through all local scales to calculate each term 
5:  $k\leftarrow {\mathcal{\mathcal{E}}}_{ij}^{Z}$  
6:  $l\leftarrow {\mathcal{\mathcal{E}}}_{ij}^{Z}$  
7:  $\stackrel{~}{c}}^{kl}\leftarrow {\stackrel{~}{c}}^{kl}+{a}_{ij}{b}_{ij$  
8:  $v}_{k}^{A}\leftarrow {v}_{k}^{A}+{a}_{ij}^{2$  
9:  $v}_{l}^{B}\leftarrow {v}_{l}^{B}+{b}_{ij}^{2$  
10:  $e}_{k}^{A}\leftarrow {e}_{k}^{A}+{a}_{ij$  
11:  $e}_{l}^{B}\leftarrow {e}_{l}^{B}+{b}_{ij$  
12:  end for  
13:  for $k:=1,\mathrm{\dots},n1$ do  iterate through each scale again and add up adjacent terms 
14:  $\stackrel{~}{c}}^{1,k+1}\leftarrow {\stackrel{~}{c}}^{1,k}+{\stackrel{~}{c}}^{1,k+1$  
15:  $\stackrel{~}{c}}^{k+1,1}\leftarrow {\stackrel{~}{c}}^{k+1,1}+{\stackrel{~}{c}}^{k+1,1$  
16:  for $Z:=A,B$ do $v}_{k+1}^{Z}\leftarrow {v}_{k}^{Z}+{v}_{k+1}^{Z$ end for  
17:  for $Z:=A,B$ do $e}_{k+1}^{Z}\leftarrow {e}_{k}^{Z}+{e}_{k+1}^{Z$ end for  
18:  end for  
19:  for $k,l:=1,\mathrm{\dots},n1$ do  
20:  $\stackrel{~}{c}}^{k+1,l+1}\leftarrow {\stackrel{~}{c}}^{k+1,l}+{\stackrel{~}{c}}^{k,l+1}+{\stackrel{~}{c}}^{k+1,l+1}{\stackrel{~}{c}}^{k,l$  
21:  end for  
22:  for $k,l:=1,\mathrm{\dots},n$ do  
23:  $c}^{kl}\leftarrow \left({\stackrel{~}{c}}^{kl}{e}_{k}^{A}{e}_{l}^{B}/{n}^{2}\right)/\sqrt{\left({v}_{k}^{A}{{e}_{k}^{A}}^{2}/{n}^{2}\right)\left({v}_{l}^{B}{{e}_{l}^{B}}^{2}/{n}^{2}\right)$  
24:  end for  
25:  end function 
Pseudocode C6 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$ MonteCarlo replicates. The power of Mgc follows by computing the test statistic under the null and the alternative using Algorithm C2. 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 $a$.  
Output: The power $\xdf$ of Mgc.  
1:  function MGCPower($F}_{XY$, $n$, $r$, $a$)  
2:  for $t:=1,\mathrm{\dots},r$ do  
3:  for $i:=[n]$ do  
4:  $x}_{i}^{0}\stackrel{iid}{\sim}{F}_{X},\phantom{\rule{thinmathspace}{0ex}}{y}_{i}^{0}\stackrel{iid}{\sim}{F}_{Y$  sample from null 
5:  $({x}_{i}^{1},{y}_{i}^{1})\stackrel{iid}{\sim}{F}_{XY},$  sample from alternative 
6:  end for  
7:  for $i,j:=1,\mathrm{\dots},n$ do  
8:  ${a}_{ij}^{0}={\delta}_{x}({x}_{i}^{0},{x}_{j}^{0})$, ${b}_{ij}^{0}={\delta}_{y}({y}_{i}^{0},{y}_{j}^{0})$  pairwise distances under the null 
9:  ${a}_{ij}^{1}={\delta}_{x}({x}_{i}^{1},{x}_{j}^{1})$, ${b}_{ij}^{1}={\delta}_{y}({y}_{i}^{1},{y}_{j}^{1})$  pairwise distances under the alternative 
10:  end for  
11:  ${c}_{0}^{\ast}[t]=\mathrm{M}\mathrm{G}\mathrm{C}\mathrm{S}\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{L}\mathrm{E}\mathrm{S}\mathrm{T}\mathrm{A}\mathrm{T}({A}^{0},{B}^{0})$  Mgc statistic under the null 
12:  ${c}_{1}^{\ast}[t]=\mathrm{M}\mathrm{G}\mathrm{C}\mathrm{S}\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{L}\mathrm{E}\mathrm{S}\mathrm{T}\mathrm{A}\mathrm{T}({A}^{1},{B}^{1})$  Mgc statistic under the alternative 
13:  end for  
14:  ${\omega}_{\alpha}\leftarrow {\mathrm{C}\mathrm{D}\mathrm{F}}_{1\alpha}({c}_{0}^{\ast}[t],t\in [r])$  the critical value of Mgc under the null 
15:  $\beta \leftarrow \sum _{t=1}^{r}({c}_{1}^{\ast}[t]>{\omega}_{\alpha})/r$  compute power by the alternative distribution 
16:  end function 
Simulation dependence functions
Request a detailed protocolThis 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 (Székely et al., 2007; Simon and Tibshirani, 2012; Gorfine et al., 2012). 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 $\mathit{x}\in {\mathbb{R}}^{p}$, we denote ${\mathit{x}}_{[d]},d=1,\dots ,p$ as the $d}^{th$ dimension of the vector $\mathit{x}$. For the purpose of highdimensional simulations, $w\in {\mathbb{R}}^{p}$ is a decaying vector with ${w}_{[d]}=1/d$ for each $d$, such that $w}^{\mathsf{T}}\mathit{x$ is a weighted summation of all dimensions of $\mathit{x}$. Furthermore, $\mathcal{\mathcal{U}}(a,b)$ denotes the uniform distribution on the interval $(a,b)$, $\mathcal{\mathcal{B}}(p)$ denotes the Bernoulli distribution with probability $p$, $\mathcal{\mathcal{N}}(\mu ,\mathrm{\Sigma})$ denotes the normal distribution with mean $\mu $ and covariance $\mathrm{S}$, $U$ and $V$ represent some auxiliary random variables, $\kappa$ is a scalar constant to control the noise level (which equals $1$ for onedimensional simulations and 0 otherwise), and $\u03f5$ is a white noise from independent standard normal distribution unless mentioned otherwise.
For all the below equations, $(X,Y)\stackrel{iid}{\sim}{F}_{XY}={F}_{YX}{F}_{X}$. For each relationship, we provide the space of $(X,Y)$, and define $F}_{Y\mathit{X}$ and ${F}_{X}$, as well as any additional auxiliary distributions.
1. Linear $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}$,
2. Exponential $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}$:
3. Cubic $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}$:
4. Joint normal $(X,Y)\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}$: Let $\rho =1/2p$, ${I}_{p}$ be the identity matrix of size $p\times p$, ${J}_{p}$ be the matrix of ones of size $p\times p$, and $\mathrm{\Sigma}=\left[\begin{array}{cc}{I}_{p}& \rho {J}_{p}\\ \rho {J}_{p}& (1+0.5\kappa ){I}_{p}\end{array}\right]$. Then
5. Step Function $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}$
where $\mathit{I}$ is the indicator function, that is $\mathit{I}(z)$ is unity whenever $z$ true, and zero otherwise.
6. Quadratic $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}$:
7. W Shape $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}:\phantom{\rule{thinmathspace}{0ex}}U\sim \mathcal{\mathcal{U}}(1,1{)}^{p}$,
8. Spiral $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}:\phantom{\rule{thinmathspace}{0ex}}U\sim \mathcal{\mathcal{U}}(0,5)$, $\u03f5\sim \mathcal{\mathcal{N}}(0,1)$
9. Uncorrelated Bernoulli $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}:\phantom{\rule{thinmathspace}{0ex}}U\sim \mathcal{\mathcal{B}}(0.5)\phantom{\rule{thinmathspace}{0ex}}{\u03f5}_{1}\sim \mathcal{\mathcal{N}}(0,{I}_{p}),\phantom{\rule{thinmathspace}{0ex}}{\u03f5}_{2}\sim \mathcal{\mathcal{N}}(0,1),$
10. Logarithmic $(X,Y)\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}:\phantom{\rule{thinmathspace}{0ex}}\u03f5\sim \mathcal{\mathcal{N}}(0,{I}_{p})$
11. Fourth Root $(X,Y)\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}:$
12. Sine Period $4\pi (X,Y)\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}:U\sim \mathcal{\mathcal{U}}(1,1),V\sim \mathcal{\mathcal{N}}(0,1{)}^{p},\theta =4\pi$,
13. Sine Period $16\pi \phantom{\rule{thinmathspace}{0ex}}(X,Y)\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}$: Same as above except $\theta =16\pi$ and the noise on $Y$ is changed to $0.5\kappa \u03f5$.
14. Square $(X,Y)\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}$: Let $U\sim \mathcal{\mathcal{U}}(1,1),\phantom{\rule{thinmathspace}{0ex}}V\sim \mathcal{\mathcal{U}}(1,1),\phantom{\rule{thinmathspace}{0ex}}\u03f5\sim \mathcal{\mathcal{N}}(0,1{)}^{p},\phantom{\rule{thinmathspace}{0ex}}\theta =\frac{\pi}{8}$. Then
for $d=1,\dots ,p.$
15. Two Parabolas $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}$: $\u03f5\sim \mathcal{\mathcal{U}}(0,1),\phantom{\rule{thinmathspace}{0ex}}U\sim \mathcal{\mathcal{B}}(0.5)$,
16. Circle $\begin{array}{l}(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}:U\sim \mathcal{\mathcal{U}}(1,1{)}^{p},\u03f5\sim \mathcal{\mathcal{N}}(0,{I}_{p}),r=1,\end{array}$
17. Ellipse $(X,Y)\in {\mathbb{R}}^{p}\times \mathbb{R}$: Same as above except $r=5$.
18. Diamond $(X,Y)\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}$: Same as 'Square' except $\theta =\frac{\pi}{4}$.
19. Multiplicative Noise $(\mathit{x},\mathit{y})\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}:u\sim \mathcal{\mathcal{N}}(0,{I}_{p}),$
20. Multimodal Independence $(X,Y)\in {\mathbb{R}}^{p}\times {\mathbb{R}}^{p}:\mathrm{L}\mathrm{e}\mathrm{t}\phantom{\rule{thinmathspace}{0ex}}U\sim \mathcal{\mathcal{N}}(0,{I}_{p}),V\sim \mathcal{\mathcal{N}}(0,{I}_{p}),$$U}^{\prime}\sim \mathcal{\mathcal{B}}(0.5{)}^{p},{V}^{\prime}\sim \mathcal{\mathcal{B}}(0.5{)}^{p$. Then
For each distribution, $X$ and $Y$ are dependent except (20); for some relationships (8,14,1618) 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 2—figure supplement 1.
For the increasing dimension simulation in the main paper, we always set $\kappa =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\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$ to make the highdimensional relationships more difficult (otherwise, additional dimensions only add more signal). For the onedimensional simulations, we always set $p=q=1$, $\kappa =1$ and $n=100$.
Appendix 1
Real data processing
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 InventoryRevised the characterized personality along five dimensions (Costa and McCrae, 1992).
This dataset consists of 42 subjects, each with 197 timesteps of restingstate functional magnetic resonance activity (rsfMRI) activity, as well as the subject’s fivedimensional 'personality'. Adelstein et al. (Adelstein et al., 2011) 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 fivefactor 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 (Craddock et al., 2013) 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 KullbackLeibler divergence across regions to obtain a similarity matrix across comparing all regions. Then, use the normalized Hellinger distance to compute distances between each subject.
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 (Jung et al., 2009).
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 (Jung et al., 2009). Previously published results explored the relationship between graphs and CCI (Koutra et al., 2015), 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 (Roncal et al., 2013). We compute the distance between the graphs using the semiparametric graph test statistic (Sussman et al., 2012; Shen et al., 2017; Tang et al., 2017), embedding each graph into two dimensions and aligning the embeddings via a Procrustes analysis.
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 (Wang et al., 2011). 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 (Wang et al., 2017). Then, using other methods, we identifed three peptides that were implicated in ovarian cancer, and validated them as legitimate biomarkers with a followup experiment.
In this study, we performed the following five sets of tests on those data:
Ovarian vs. normal for all proteins,
Ovarian vs. normal for each individual protein,
Pancreas vs. normal for all proteins,
Pancreas vs. all others for each individual protein,
Pancreas vs. normal for each individual protein.
These tests are designed to first validate the Mgc method from ovarian cancer, then identify biomarkers unique to pancreatic cancer, that is, find a protein that is able to tell the difference between pancreas and normals, as well as pancreas vs all other cancers. For each of the five tests, we create a binary label vector, with $1$ indicating the cancer type of interest for the corresponding subject, and 0 otherwise. Then each algorithm is applied to each task. For all tests we used Euclidean distances and the type 1 error level is set to $a=0.05$. The three test sets assessing individual proteins provide 318 pvalues; we used the BenjaminiHochberg procedure (Benjamini and Hochberg, 1995) to control the false discovery rate. A summary of the results are reported in Appendix 1—table 1.

Appendix 1—table 1—source data 1
 https://doi.org/10.7554/eLife.41690.021

Appendix 1—table 1—source data 2
 https://doi.org/10.7554/eLife.41690.022

Appendix 1—table 1—source data 3
 https://doi.org/10.7554/eLife.41690.023
All methods are able to successfully detect a dependence between peptide abundances in ovarian cancer samples versus normal samples (Appendix 1—table 1, line 1). This is likely because there are so many individual peptides that have different abundance distributions between ovarian and normal samples (Appendix 1—table 1, line 2). Nonetheless, Mgc identified more putative biomarkers than any of the other methods. While we have not checked all of them with subsequent experiments to identify potential false positives, we do know from previous experiments that three peptides in particular are effective biomarkers.
All three peptides have pvalue ≈ 0 for all methods including Mgc, that is, they are all correctly identified as significant. However, by ranking the peptides based on the actual test statistic of each peptide, Mgc is the method that ranks the three known biomarkers the lowest, suggesting that it is the least likely to falsely identify peptides.
We then investigated the pancreatic samples in an effort to identify biomarkers that are unique to pancreas. We first checked whether the methods could identify a difference using all the peptides. Indeed, three methods found a dependence at the 0.05 level, with Mgc obtaining the lowest pvalue (Appendix 1—table 1, line 3). We then investigated how many individual peptides the methods identified; all of them found 6 to 11 peptides with a significant difference between pancreatic and normal samples (Appendix 1—table 1, line 4). Because we were interested in identifying peptides that were uniquely useful for pancreatic cancer, we then compared pancreatic samples to all others. At significance level $0.05$, only Mgc, Hsic, and Hhg identified peptides that expressed different abundances in this more challenging case, and we list the top four peptides in Appendix 1—table 2 along with the corrected pvalue for each peptide.
All three methods reveal the same unique protein for pancreas: neurogranin. Hsic identifies another peptide (tropomyosin alpha3 chain isoform 4), and Hhg identifies a third peptide (fibrinogenlike protein 1 precursor). However, fibrinogenlike protein 1 precursor is not significant for pvalue testing between pancreatic and normal subjects. On the other hand, tropomyosin is a ubiquitously expressed protein, since normal tissues and other cancers will also express tropomyosin and leak it into blood, whereas neurogranin is exclusively expressed only in brain tissues. Moreover, there exists strong evidence of tropomyosin 3 upregulated in other cancers (Karsani et al., 2014; Sun et al., 2016; Lee et al., 2012; Lam et al., 2012). Therefore, it suggests that the other two peptides identified by Hhg and Hsic are likely false positives.
In fact, neurogranin is always one of the top 4 candidates in all methods except Mic; the only difference is that the corrected pvalues are not significant enough for other methods. Along with the classification result in Figure 4D showing that neurogranin alone has the best classification error, Mgc discovers an ideal candidate for potential biomarker. Moreover, the fact that Mgc, Hhg and Hsic discover the dependency while others cannot implies a nonlinear relationship.
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. (Eklund et al., 2012; Eklund et al., 2016), 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 timevarying 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 1000 Functional Connectomes Project consisting of a total of 1583 subjects (Biswal et al., 2010). Appendix 1—figure 1 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 preprocessing details (Eklund et al., 2012; Eklund et al., 2016). Moreover, even the proposed solutions to those issues make linearity assumptions, thereby limiting detection to only a small subset of possible dependence functions.
Running time report in experiments
Appendix 1—table 3 lists the actual running time of Mgc versus other methods for testing on the real data, based on a modern desktop with a six core I76850K CPU and 32 GB memory on MATLAB 2017a on Windows 10. The first two experiments are timed based on 1000 permutations, while the screening experiment is timed without permutation, that is compute the test statistic only. Pearson runs the fastest, trailed by Mic and then Dcorr. Pearson and Mic are only possible to run in the screening experiment, as the other two experiments are multivariate. The running time of Mgc is a constant times (about 10) higher than that of Dcorr, and Hhg is implemented in a running time of $O({n}^{3})$ and thus significantly slower.
Data availability
To facilitate reproducibility, we make all datasets available from: https://github.com/neurodata/MGCpaper/tree/master/Data/Preprocessed (copy archived at https://github.com/elifesciencespublications/MGCpaper).
References

Multiscale geometric methods for data sets II: geometric MultiResolution analysisApplied and Computational Harmonic Analysis 32:435–462.https://doi.org/10.1016/j.acha.2011.08.001

Controlling the false discovery rate: a practical and powerful approach to multiple testingJournal of the Royal Statistical Society: Series B 57:289–300.https://doi.org/10.1111/j.25176161.1995.tb02031.x

Valid postselection inferenceThe Annals of Statistics 41:802–837.https://doi.org/10.1214/12AOS1077

Advances in biomarker research for pancreatic cancerCurrent Pharmaceutical Design 18:2439–2451.https://doi.org/10.2174/13816128112092439

Diffusion waveletsApplied and Computational Harmonic Analysis 21:53–94.https://doi.org/10.1016/j.acha.2006.04.004

ConferenceFlashGraph: processing BillionNode graphs on an array of commodity SSDsUSENIX Conference on File and Storage Technologies.

Testing predictions from personality neuroscience. brain structure and the big fivePsychological Science 21:820–828.https://doi.org/10.1177/0956797610370159

Clinical proteomic biomarkers: relevant issues on study design & technical considerations in biomarker developmentClinical and Translational Medicine 3:7.https://doi.org/10.1186/2001132637

GraphTheoretic measures of multivariate association and predictionThe Annals of Statistics 11:377–391.https://doi.org/10.1214/aos/1176346148

Advances in Neural Information Processing System2672–2680, Generative adversarial nets, Advances in Neural Information Processing System, MIT Press.

ReportComment on Detecting Novel Associations in Large Data SetsIsrael Institute of Technology.

Advances in Neural Information Processing Systems513–520, A kernel method for the twosampleproblem, Advances in Neural Information Processing Systems, MIT Press.

Consistent nonparametric tests of independenceJournal of Machine Learning Research 11:1391–1423.

Dismantling the mantel testsMethods in Ecology and Evolution 4:336–344.https://doi.org/10.1111/2041210x.12018

Tropomyosin as a regulator of cancer cell transformationAdvances in Experimental Medicine and Biology 644:124–131.https://doi.org/10.1007/9780387857664_10

Consistent distributionfree sample and independence tests for univariate random variablesJournal of Machine Learning Research 17:1–54.

A NonParametric test of independenceThe Annals of Mathematical Statistics 19:546–557.https://doi.org/10.1214/aoms/1177730150

Relations between two sets of variatesBiometrika 28:321–377.https://doi.org/10.1093/biomet/28.34.321

Fast computing for distance covarianceTechnometrics 58:435–447.https://doi.org/10.1080/00401706.2015.1054435

Canonical analysis of several sets of variablesBiometrika 58:433–451.https://doi.org/10.1093/biomet/58.3.433

DeltaCon: a principled massivegraph similarity functionACM Transactions on Knowledge Discovery From Data 10:.https://doi.org/10.1145/2824443

Comparison of protein expression profiles of different stages of lymph nodes metastasis in breast cancerInternational Journal of Biological Sciences 8:353–362.https://doi.org/10.7150/ijbs.3157

Advances in Neural Information Processing SystemsMaximum likelihood estimation of intrinsic dimension, Advances in Neural Information Processing Systems, MIT Press.

Feature screening via distance correlation learningJournal of the American Statistical Association 107:1129–1139.https://doi.org/10.1080/01621459.2012.695654

Distance covariance in metric spacesThe Annals of Probability 41:3284–3305.https://doi.org/10.1214/12AOP803

The detection of disease clustering and a generalized regression approachCancer Research 27:209–220.

Kernel mean embedding of distributions: A review and beyondFoundations and Trends in Machine Learning 10:1–141.https://doi.org/10.1561/2200000060

Notes on regression and inheritance in the case of two parentsProceedings of the Royal Society of London 58:240–242.https://doi.org/10.1098/rspl.1895.0041

On quantifying dependence: a framework for developing interpretable measuresStatistical Science 28:116–130.https://doi.org/10.1214/12STS405

On measures of dependenceActa Mathematica Academiae Scientiarum Hungaricae 10:441–451.https://doi.org/10.1007/BF02024507

DISCO analysis: a nonparametric extension of analysis of varianceThe Annals of Applied Statistics 4:1034–1055.https://doi.org/10.1214/09AOAS245

Energy distanceWiley Interdisciplinary Reviews: Computational Statistics 8:27–38.https://doi.org/10.1002/wics.1375

ConferenceMIGRAINE: mri graph reliability analysis and inference for connectomicsGlobal Conference on Signal and Information Processing.

The big five default brain: functional evidenceBrain Structure and Function 219:1913–1922.https://doi.org/10.1007/s004290130610y

Multivariate TwoSample tests based on nearest neighborsJournal of the American Statistical Association 81:799–806.https://doi.org/10.1080/01621459.1986.10478337

Equivalence of distancebased and RKHSbased statistics in hypothesis testingThe Annals of Statistics 41:2263–2291.https://doi.org/10.1214/13AOS1140

Generalized canonical correlation analysis for classificationJournal of Multivariate Analysis 130:310–322.https://doi.org/10.1016/j.jmva.2014.05.011

Manifold matching using shortestpath distance and joint neighborhood selectionPattern Recognition Letters 92:41–48.https://doi.org/10.1016/j.patrec.2017.04.005

From distance correlation to multiscale graph correlationJournal of the American Statistical Association pp. 1–39.https://doi.org/10.1080/01621459.2018.1543125

The proof and measurement of association between two thingsThe American Journal of Psychology 15:72.https://doi.org/10.2307/1412159

A consistent adjacency spectral embedding for stochastic blockmodel graphsJournal of the American Statistical Association 107:1119–1128.https://doi.org/10.1080/01621459.2012.699795

ConferenceGenerative models and model criticism via optimized maximum mean discrepancyInternational Conference on Learning Representations.

Measuring and testing dependence by correlation of distancesThe Annals of Statistics 35:2769–2794.https://doi.org/10.1214/009053607000000505

A new test for multivariate normalityJournal of Multivariate Analysis 93:58–80.https://doi.org/10.1016/j.jmva.2003.12.002

Brownian distance covarianceThe Annals of Applied Statistics 3:1236–1265.https://doi.org/10.1214/09AOAS312

The distance correlation ttest of independence in high dimensionJournal of Multivariate Analysis 117:193–213.https://doi.org/10.1016/j.jmva.2013.02.012

Partial distance correlation with methods for dissimilaritiesThe Annals of Statistics 42:2382–2412.https://doi.org/10.1214/14AOS1255

A semiparametric TwoSample hypothesis testing problem for random graphsJournal of Computational and Graphical Statistics 26:344–354.https://doi.org/10.1080/10618600.2016.1193505

Regularized generalized canonical correlation analysisPsychometrika 76:257–284.https://doi.org/10.1007/s1133601192068

Conditional distance correlationJournal of the American Statistical Association 110:1726–1734.https://doi.org/10.1080/01621459.2014.993081

Penalized classification using Fisher's linear discriminantJournal of the Royal Statistical Society: Series B 73:753–772.https://doi.org/10.1111/j.14679868.2011.00783.x

Distance metric learning with application to clustering with sideinformationAdvances in Neural Information Processing Systems 15:505–512.

Serum neurogranin measurement as a biomarker of acute traumatic brain injuryClinical Biochemistry 48:843–848.https://doi.org/10.1016/j.clinbiochem.2015.05.015

A simple statistical parameter for use in evaluation and validation of high throughput screening assaysJournal of Biomolecular Screening 4:67–73.https://doi.org/10.1177/108705719900400206

Adaptive manifold learningIEEE Transactions on Pattern Analysis and Machine Intelligence 34:253–265.https://doi.org/10.1109/TPAMI.2011.115

Largescale kernel methods for independence testingStatistics and Computing 28:113–130.https://doi.org/10.1007/s1122201697217

An iterative approach to distance correlationbased sure independence screeningJournal of Statistical Computation and Simulation 85:2331–2345.https://doi.org/10.1080/00949655.2014.928820
Article and author information
Author details
Funding
Child Mind Institute (Endeavor Scientist Program)
 Joshua T Vogelstein
National Science Foundation
 Joshua T Vogelstein
Defense Advanced Research Projects Agency
 Joshua T Vogelstein
Office of Naval Research
 Joshua T Vogelstein
Air Force Office of Scientific Research
 Joshua T Vogelstein
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was partially supported by the Child Mind Institute Endeavor Scientist Program, the National Science Foundation award DMS1712947, the National Security Science and Engineering Faculty Fellowship (NSSEFF), the Johns Hopkins University Human Language Technology Center of Excellence (JHU HLT COE), the Defense Advanced Research Projects Agency’s (DARPA) SIMPLEX program through SPAWAR contract N6600115C4041, the XDATA program of DARPA administered through Air Force Research Laboratory contract FA87501220303, DARPA Lifelong Learning Machines program through contract FA86501827834, the Office of Naval Research contract N000141210601, the Air Force Office of Scientific Research contract FA95501410033. The authors thank Dr. Brett 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. And we’d like to thank Satish Palaniappan, Sambit Panda, and Junhao (Bear) Xiong for porting the code to Python.
Version history
 Received: September 3, 2018
 Accepted: January 14, 2019
 Accepted Manuscript published: January 15, 2019 (version 1)
 Version of Record published: February 22, 2019 (version 2)
Copyright
© 2019, Vogelstein et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 3,524
 views

 419
 downloads

 18
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
 Computational and Systems Biology
Bats have unique characteristics compared to other mammals, including increased longevity and higher resistance to cancer and infectious disease. While previous studies have analyzed the metabolic requirements for flight, it is still unclear how bat metabolism supports these unique features, and no study has integrated metabolomics, transcriptomics, and proteomics to characterize bat metabolism. In this work, we performed a multiomics data analysis using a computational model of metabolic fluxes to identify fundamental differences in central metabolism between primary lung fibroblast cell lines from the black flying fox fruit bat (Pteropus alecto) and human. Bat cells showed higher expression levels of Complex I components of electron transport chain (ETC), but, remarkably, a lower rate of oxygen consumption. Computational modeling interpreted these results as indicating that Complex II activity may be low or reversed, similar to an ischemic state. An ischemiclike state of bats was also supported by decreased levels of central metabolites and increased ratios of succinate to fumarate in bat cells. Ischemic states tend to produce reactive oxygen species (ROS), which would be incompatible with the longevity of bats. However, bat cells had higher antioxidant reservoirs (higher total glutathione and higher ratio of NADPH to NADP) despite higher mitochondrial ROS levels. In addition, bat cells were more resistant to glucose deprivation and had increased resistance to ferroptosis, one of the characteristics of which is oxidative stress. Thus, our studies revealed distinct differences in the ETC regulation and metabolic stress responses between human and bat cells.

 Computational and Systems Biology
 Neuroscience
Normal aging leads to myelin alterations in the rhesus monkey dorsolateral prefrontal cortex (dlPFC), which are positively correlated with degree of cognitive impairment. It is hypothesized that remyelination with shorter and thinner myelin sheaths partially compensates for myelin degradation, but computational modeling has not yet explored these two phenomena together systematically. Here, we used a twopronged modeling approach to determine how agerelated myelin changes affect a core cognitive function: spatial working memory. First, we built a multicompartment pyramidal neuron model fit to monkey dlPFC empirical data, with an axon including myelinated segments having paranodes, juxtaparanodes, internodes, and tight junctions. This model was used to quantify conduction velocity (CV) changes and action potential (AP) failures after demyelination and subsequent remyelination. Next, we incorporated the single neuron results into a spiking neural network model of working memory. While complete remyelination nearly recovered axonal transmission and network function to unperturbed levels, our models predict that biologically plausible levels of myelin dystrophy, if uncompensated by other factors, can account for substantial working memory impairment with aging. The present computational study unites empirical data from ultrastructure up to behavior during normal aging, and has broader implications for many demyelinating conditions, such as multiple sclerosis or schizophrenia.