Discovering and deciphering relationships across disparate data modalities

  1. Joshua T Vogelstein  Is a corresponding author
  2. Eric W Bridgeford
  3. Qing Wang
  4. Carey E Priebe
  5. Mauro Maggioni
  6. Cencheng Shen
  1. Johns Hopkins University, United States
  2. Child Mind Institute, United States
  3. University of Delaware, United States


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

eLife 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.


Identifying the existence of a relationship between a pair of properties or modalities is the critical initial step in data science investigations. Only if there is a statistically significant relationship does it make sense to try to decipher the nature of the relationship. Discovering and deciphering relationships is fundamental, for example, in high-throughput screening (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 Product-Moment 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 dependence-testing that were not addressed in Pearson’s era. First, we now desire methods that can correctly detect any kind of dependence between all kinds of data, including high-dimensional data (such as ’omics), structured data (such as images or networks), with nonlinear relationships (such as oscillators), even with very small sample sizes as is common in modern biomedical science. Second, we desire methods that are interpretable by providing insight into how or why they discovered the presence of a statistically significant relationship. Such insight can be a crucial component of designing the next computational or physical experiment.

While many statistical and machine learning approaches have been developed over the last 120 years to combat aspects of the first issue—detecting dependencies—no approach satisfactorily addressed the challenges across all data types, relationships, and dimensionalities. Hoeffding and Renyi proposed non-parametric tests to address nonlinear but univariate relationships (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 high-dimensional linear data, whereas the kernel methods and another test (Heller, Heller, and Gorfine’s test, Hhg) (Heller et al., 2013) perform well on low-dimensional nonlinear data. But no test performs particularly well on high-dimensional nonlinear data with typical sample sizes, which characterizes a large fraction of real data challenges in the current big data era.

Moreover, to our knowledge, 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 high-dimensional, nonlinear, low-sample 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). Post-hoc cross-validation 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.

Illustration of Multiscale Graph Correlation (Mgc) on simulated cloud density (xi) and grass wetness (yi).

We present two different relationships: linear (top) and nonlinear spiral (bottom; see Materials and methods for simulation details). (A) Scatterplots of the raw data using 50 pairs of samples for each scenario. Samples 1, 2, and 3 (black) are highlighted; arrows show x distances between these pairs of points while their y distances are almost 0. (B) Scatterplots of all pairs of distances comparing x and y distances. Distances are linearly correlated in the linear relationship, whereas they are not in the spiral relationship. Dcorr uses all distances (gray dots) to compute its test statistic and p-value, whereas Mgc chooses the local scale and then uses only the local distances (green dots). (C) Heatmaps characterizing the strength of the generalized correlation at all possible scales (ranging from 2 to n for both x and y). For the linear relationship, the global scale is optimal, which is the scale that Mgc selects and results in a p-value identical to Dcorr. For the nonlinear relationship, the optimal scale is local in both x and y, so Mgc achieves a far larger test statistic, and a correspondingly smaller and significant p-value. Thus, Mgc uniquely detects dependence and characterizes the geometry in both relationships.

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 high-dimensional data, where simple visualizations do not reveal relationships to the unaided human eye.

Our method, ‘Multiscale Graph Correlation’ (Mgc, pronounced ‘magic’), generalized and extends previously proposed pairwise comparison-based approaches by adaptively estimating the informative scales for any relationship — linear or nonlinear, low-dimensional or high-dimensional, unstructured or 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 hypothesis-testing and insight-providing approach that builds on recent developments in manifold and kernel learning, with complementary developments in nearest-neighbor search, and multiscale analyses.

The multiscale graph correlation procedure

Mgc is a multi-step 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):

  1. Compute two distance matrices, one consisting of distances between all pairs of one property (e.g. cloud densities, entire genomes or connectomes) and the other consisting of distances between all pairs of the other property (e.g. grass wetnesses or disease status). Then center each distance matrix (by subtracting its overall mean, the column-wise mean from each column, and the row-wise mean from each row), and denote the resulting n-by-n matrices A and B.

  2. For all possible values of k and l from 1 to n:

    1. Compute the k-nearest neighbor graphs for one property, and the l-nearest neighbor graphs for the other property. Let Gk and Hl be the adjacency matrices for the nearest neighbor graphs, so that Gk(i,j)=1 indicates that A(i,j) is within the k smallest values of the ith row of A, and similarly for Hl.

    2. Estimate the local correlations—the correlation between distances restricted to only the (k,l) neighbors—by summing the products of the above matrices, ckl=ijA(i,j)Gk(i,j)B(i,j)Hl(i,j).

    3. Normalize ckl such that the result is always between -1 and +1 by dividing byijA2(i,j)Gk(i,j)×ijB2(i,j)Hl(i,j).

  3. Estimate the optimal local correlation c* by finding the smoothed maximum of all local correlations {ckl}. Smoothing mitigates biases and provides Mgc with theoretical guarantee and better finite-sample performance.

  4. 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 p-value, rather than one p-value per scale, ensuring that it is a valid test (meaning that the false positive rate is properly controlled at the specified type I error rate).

Computing all local correlations, the test statistic, and the p-value requires O(n2logn) time, which is about the same running time complexity as other methods (Shen et al., 2018).


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):

H0:XandYare independentHA:XandYare not independent.

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(H0 is rejected |H0 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 distance-based test that is very powerful for detecting low-dimensional 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 trace-norm of the cross-covariance matrix, and (vii) the Cca method, which is the largest (in magnitude) singular value of the cross-covariance matrix, and can be viewed as a different generalization of Pearson in high-dimensions that is more appropriate for sparse settings (Hotelling, 1936; Witten et al., 2009; Witten and Tibshirani, 2011). Note that while we focus on high-dimensional settings, Appendix 1 shows further results in one-dimensional settings, also comparing to a number of tests that are limited to one dimension, including: (viii) Pearson’s product moment correlation, (ix) Spearman’s rank correlation (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, W-shape), 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 high-dimensional variants, to more extensively evaluate the methods; function details are in Materials and methods. The visualization of one-dimensional noise-free (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 high-dimensional relationships are more challenging because (1) they cannot be easily visualized and (2) each dimension is designed to have less and less signal, so there are many noisy dimensions. Figure 2 shows that Mgc achieves the highest (or close to the highest) power given 100 samples for each relationship and dimensionality. Figure 2—figure supplement 2 shows the same advantage in one-dimension with increasing sample size.

Figure 2 with 4 supplements see all
An extensive benchmark suite of 20 different relationships spanning polynomial, trigonometric, geometric, and other relationships demonstrates that Mgc empirically nearly dominates eight other methods across dependencies and dimensionalities ranging from 1 to 1000 (see Materials and methods and Figure 2—figure supplement 1 for details).

Each panel shows the testing power of other methods relative to the power of Mgc (e.g. power of Mcorr minus the power of Mgc) at significance level α=0.05 versus dimensionality for n=100. Any line below zero at any point indicates that that method’s power is less than Mgc’s power for the specified setting and dimensionality. Mgc achieves empirically better (or similar) power than all other methods in almost all relationships and all dimensions. For the independent relationship (#20), all methods yield power 0.05 as they should. Note that Mgc is always plotted ‘on top’ of the other methods, therefore, some lines are obscured.

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 non-monotone 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, distance-based methods including Mcorr, Dcorr, Mgc, Hhg and Hsic are slightly worse, while Mic and Mantel are the worst. Mgc’s performance is equal to linear methods on monotonic relationships. For non-monotonic relationships, traditional correlations fail to detect the existence of dependencies, Dcorr, Mcorr, and Mic, do reasonably well, but Hhg and Mgc require the fewest samples. In the high-dimensional non-monotonic relationships that motivated this work, and 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 one-dimensional settings while varying sample size. Thus, not only does Mgc empirically nearly dominate existing tests, it is a framework that one can apply to future tests to further improve their performance.

Table 1
The median sample size for each method to achieve power 85% at type one error level 0.05, grouped into monotone (type 1–5) and non-monotone relationships (type 6–19) for both one- and ten-dimensional settings, normalized by the number of samples required by Mgc.

In other words, a 2.0 indicates that the method requires double the sample size to achieve 85% power relative to Mgc. Pearson, Rv, and Cca all achieve the same performance, as do Spearman and Kendall. Mgc requires the fewest number of samples in all settings, and for high-dimensional non-monotonic relationships, all other methods require about double or triple the number of samples Mgc requires.
Dependency typeMonotoneNon-MonoAverageMonotoneNon-MonoAverage
Pearson / Rv / Cca1>10>100.8>10>10
Spearman / Kendall1>10>10n/an/an/a
Table 1—source data 1

Testing power sample size data in one dimension.
Table 1—source data 2

Testing power sample size data in high-dimensions.

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 Mgc-Map shows local correlation as a function of the scales of the two properties. More concretely, it is the matrix of ckl’s, as defined above. Thus, the Mgc-Map is an n-by-n matrix which encodes the strength of dependence for each possible scale. Figure 3 provides the Mgc-Map for all 20 different one-dimensional relationships; the optimal scale to achieve t^* is marked with a green dot. For the monotonic dependencies (1-5), the optimal scale is always the largest scale, that is the global one. For all non-monotonic dependencies (6-19), Mgc chooses smaller scales. Thus, a global optimal scale implies a close-to-linear dependency, otherwise the dependency is strongly nonlinear. In fact, this empirical observation led to the following theorem (which is proved in 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 Mgc-Map 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.

Figure 3 with 1 supplement see all
The Mgc-Map characterizes the geometry of the dependence function.

For each of the 20 panels, the abscissa and ordinate denote the number of neighbors for X and Y, respectively, and the color denotes the magnitude of each local correlation. For each simulation, the sample size is 60, and both X and Y are one-dimensional. Each dependency has a different Mgc-Map characterizing the geometry of dependence, and the optimal scale is shown in green. In linear or close-to-linear relationships (first row), the optimal scale is global, that is the green dot is in the top right corner. Otherwise the optimal scale is non-global, which holds for the remaining dependencies. Moreover, similar dependencies often share similar Mgc-Maps and similar optimal scales, such as (10) logarithmic and (11) fourth root, the trigonometric functions in (12) and (13 , 16) circle and (17) ellipse, and (14) square and (18) diamond. The Mgc-Maps for high-dimensional simulations are provided in Figure 3—figure supplement 1.

Moreover, similar dependencies have similar Mgc-Maps 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 Mgc-Maps. 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 Mgc-Maps. This indicates that the Mgc-Map partially characterizes the geometry of these relationships, differentiating different dependence structures and assisting subsequent analysis steps. Moreover, in Shen and Vogelstein, 2018, we proved that the sample Mgc-Map (which Mgc estimates) converges to the true Mgc-Map provided by the underlying joint distribution of the data. In other words, each relationship has a specific map that characterizes it based on its joint distribution, and Mgc is able to accurately estimate it via sample observations. The existence of a population level characterization of the joint distribution strongly differentiates Mgc from previously proposed multi-scale geometric or topological characterizations of data, such as persistence diagrams (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 𝒪(n4) operations, we devised a nested implementation that requires only 𝒪(n2logn) operations. Moreover, obtaining the Mgc-Map 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 𝒪(n2logn/T), where T is the number of cores (see Algorithm C1 for details). Since T is often larger than logn, in practice, Mgc can be 𝒪(n2), 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 𝒪(n3), which is orders of magnitude slower when n>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 resting-state functional magnetic resonance (rs-fMRI) activity and personality (Adelstein et al., 2011). The first row of Table 2 compares the p-value of different methods, and Figure 4A shows the Mgc-Map for the sample data. Mgc is able to yield a significant p-value (< 0.05), whereas all previously proposed global dependence tests under consideration (Mantel, Dcorr, Mcorr, or Hhg) fail to detect dependence at a significance level of 0.05. Moreover, the Mgc-Map provides a characterization of the dependence, for which the optimal scale indicates that the dependency is strongly nonlinear. Interestingly, the Mgc-Map does not look like any of the 20 images from the simulated data, suggesting that the nonlinearity characterizing this dependency is more complex or otherwise different from those we have considered so far.

Demonstration that Mgc successfully detects dependency, distinguishes linearity from nonlinearity, and identifies the most informative feature in a variety of real data experiments.

(A) The Mgc-Map for brain activity versus personality. Mgc has a large test statistic and a significant p-value at the optimal scale (13, 4), while the global counterpart is non-significant. That the optimal scale is non-global implies a strongly nonlinear relationship. (B) The Mgc-Map for brain connectivity versus creativity. The image is similar to that of a linear relationship, and the optimal scale equals the global scale, thus both Mgc and Mcorr are significant in this case. (C) For each peptide, the x-axis shows the p-value for testing dependence between pancreatic and healthy subjects by Mgc, and the y-axis shows the p-value for testing dependence between pancreatic and all other subjects by Mgc. At critical level 0.05, Mgc identifies a unique protein after multiple testing adjustment. (D) The true and false positive counts using a k-nearest neighbor (choosing the best k[1,10]) leave-one-out classification using only the significant peptides identified by each testing method. The peptide identified by Mgc achieves the best true and false positive rates, as compared to the peptides identified by Hsic or Hhg.
Table 2
The p-values for brain imaging vs mental properties.

Mgc always uncovers the existence of significant relationships and discovers the underlying optimal scales. Bold indicates significant p-value per dataset.
Testing Pairs/MethodsMgcDcorrMcorrHhgHsic
Activity vs Personality0.0430.6670.4410.0590.124
Connectivity vs Creativity0.0110.0100.0110.0310.092
Table 2—source data 1

p-value data for activity vs personality.
Table 2—source data 2

p-value data for connetivity vs creativity.

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 Mgc-Map in Figure 4B closely resembles a linear relationship where the optimal scale is global. The close-to-linear relationship is also supported from the p-value 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 p-value, implying a relatively weak and noisy relationship. A prediction of creativity via linear regression turns out to be non-significant, 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 high-dimensional 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 p-value, followed by Hhg and Hsic and then all other methods; the second dataset exhibits a close-to-linear 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 p-value of each peptide assigned by Mgc, which reveals one particular protein, neurogranin, that exhibits a strong dependency specifically with pancreatic cancer. Subsequent literature searches reveal that neurogranin is a potentially valuable biomarker for pancreatic cancer because it is exclusively expressed in brain tissue among normal tissues and has not been linked with any other cancer type (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 k-nearest-neighbor classifier to predict pancreatic cancer, and a leave-one-subject-out validation. Figure 4D shows that the peptide selected by Mgc (neurogranin) works better than any other subset of the peptides selected by Hsic or Hhg, in terms of both fewer false positives and negatives. This analysis suggests Mgc can effectively be used for screening and subsequent classification.


There are a number of connections between Mgc and other prominent statistical procedures that may be worth further exploration. First, Mgc can be thought of as a regularized or sparsified variant of distance or kernel methods. Regularization is central to high-dimensional and ill-posed problems, where dimensionality is larger than sample size. Second, Mgc can also be thought of as learning a metric because it chooses the optimal scale amongst a set of n2 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 two-dimensional parameter space to optimize the metric, deep learning searches over a much larger parameter space, sometimes including millions of dimensions. Probably neither is optimal, and somewhere between the two would be useful in many tasks. Third, energy statistics provide state of the art approaches to other problems, including goodness-of-fit (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 two-sample (or generally the K-sample) 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 K-sample testing from two-sample 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 p-values, then the optimal p-value 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 p-value. 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 ‘post-selection inference’ problem (Berk et al., 2013). The most straightforward way to address this issue is to collect new data, which is costly and time-consuming. Therefore, researchers often ignore this fact and make statistically invalid claims. Mgc circumvents this dilemma by carefully constructing its permutation test to estimate the scale in the process of estimating a p-value, rather than after. To our knowledge, Mgc is the first dependence test to take a step towards valid post-selection inference.

As a separate next theoretical extension, we could reduce the computational space and time required by Mgc. Mgc currently requires space and time that are quadratic with respect to the number of samples, which can be costly for very large data. Recent advances in related work demonstrated that one could reduce computational time of distance-based 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, semi-external 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 (Bridgeford et al., 2018; experiments archived at 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 distance-based 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 protocol

Given pairs of observations (xi,yi)Rp×Rq for i=1,,n, assume they are independently identically distributed as (X,Y)iidFXY. If the two random variables X and Y are independent, the joint distribution equals the product of the marginals, that is FXY=FXFY. 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 x’s, q is the dimension for 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 one-dimensional 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 protocol

To achieve consistent testing, most state-of-the-art dependence measures operate on pairwise comparisons, either similarities (such as kernels) or dissimilarities (such as distances).

Let 𝒳n={x1,,xn}Rp×n and 𝒴n={y1,,yn}Rq×n denote the matrices of sample observations, and dx be the distance function for x’s and dy for y’s. One can then compute two n×n distance matrices A~={a~ij} and B~={b~ij}, where a~ij=δx(xi,xj) and b~ij=δy(yi,yj). A common example of the distance function is the Euclidean metric (L2 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 A~ and B~, respectively. Any ‘generalized correlation coefficient’ (Spearman, 1904; Kendall, 1970) can be written as:

(1) c(𝒳n,𝒴n)=1zi=1nj=1naijbij,

where z is proportional to the standard deviations of A and B, that is z=n2σaσ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. Distance-based 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 A~ and B~ by certain centering schemes into A and B. Hsic chooses the Gaussian kernel and computes two kernel matrices, then transform the kernel matrices A~ and B~ by the same centering scheme as Dcorr. For Mgc, A and B are always distance matrices (or can be transformed to distances from kernels by 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 p-value can be computed by comparing the correlation of the sample data to the correlation of the permuted sample data. The independence hypothesis is rejected if the p-value is lower than a pre-determined type 1 error level, say 0.05. Then the power of the test statistic equals the probability of a correct rejection at a specific type 1 error level. Note that Hhg is the only exception that cannot be cast as a generalized correlation coefficient, but the permutation testing is similarly effective for the Hhg test statistic; also note that the iid assumption is critical for permutation test to be valid, which may not be applicable in special cases like auto-correlated time series (Guillot and Rousset, 2013).

Distance correlation (Dcorr) and the Unbiased Version (Mcorr)

Request a detailed protocol

Define the row and column means of A~ by a¯j=1ni=1na~ij and a¯i=1nj=1na~ij. Dcorr defines

aij={a~ija¯ia¯j+a¯, if ij,0, if i=j,

and similarly for bij. For distance correlation, the numerator of Equation 1 is named the distance covariance (Dcov), while sa and sb 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 FX and FY, and show that the population distance correlation equals zero if and only if X and Y are independent, for any joint distribution FXY of finite second moments and finite dimensionality. They also show that as n, the sample distance correlation converges to the population distance correlation, that is, c(𝒳n,𝒴n)c(X,Y). Thus the sample distance correlation is consistent against any dependency of finite second moments and dimensionality. Of note, the distance covariance, distance variance, and distance correlation are always non-negative. Moreover, the consistency result holds for a much larger family of metrics, those of strong negative type (Lyons, 2013).

It turns out that the sample distance correlation has a finite-sample 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. By excluding the diagonal entries and slightly modifies the off-diagonal entries of 𝒜 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(x,y) for all p,q,n, which is approximately normal even if p,q. Thus it enjoys the same theoretical consistency as Dcorr and always has zero mean under independence.

Local correlations

Request a detailed protocol

Given any matrices A and B, we can define a set of local correlations as follows. Let R(Aj,i) be the ‘rank’ of xi relative to xj, that is, R(Aj,i)=k if xi is the kth closest point (or ‘neighbor’) to xj, as determined by ranking the n-1 distances to xj. Define R(Bi,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 xi and any neighborhood size l around each yj, we define the local pairwise comparisons:

(2) a~ijk={aij,if R(Aj,i)k,0,otherwise;b~ijl={bij,if R(Bi,j)l,0,otherwise;

and then let aijk=a~ijka¯k, where a¯k is the mean of {a~ijk}, and similarly for bijl.

The local correlation coefficient at a given scale is defined to effectively exclude large distances:

(3) ckl(𝒳n,𝒴n)=1zkli,j=1naijkbijl,

where zkl=n2σakσbl, with sak and sbl is the standard deviations for the truncated pairwise comparisons. The Mgc-Map 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 aij and bij from Equation 1.

As most nonlinear relationships intrinsically exhibit a local linear structure, considering the nearest-neighbors 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 ckl 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 nearest-neighbor may significantly improve the performance over global correlation, while considering the other scenarios does not.

Mgc as the optimal local correlation

Request a detailed protocol

We 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:

R=Largest Connected Component of {(k,l) such thatckl>max(τ,cnn)}.

Then the optimal scale equals all scales within R whose local correlations are as large as c. The choice of τ is made explicit in the pseudo-code, with further discussion and justification offered in Shen and Vogelstein, 2018.

Proof for theorem 1

Request a detailed protocol

Theorem 1. When (X,Y) are linearly related (rotation, scaling, translation, reflection), the optimal scale of Mgc equals the global scale. Conversely, that. the optimal scale is local implies a nonlinear relationship.

Proof. It suffices to prove the first statement, then the second statement follows by contrapositive. When (X,Y) are linearly related, Y=WX+b for a unitary matrix W and a constant b up-to possible scaling, in which case the distances are preserved, that is yiyj=WxiWxj=xixj. It follows that Mcorr(𝒳n,𝒴n)=1, so the global scale achieves the maximum possible correlation, and the largest connected region R is empty. Thus the optimal scale is global and Mgc(𝒳n,𝒴n)=Mcorr(𝒳n,𝒴n)=1.

Computational complexity of each step

Request a detailed protocol

The distance computation takes 𝒪(n2max{p,q}), and the ranking process takes 𝒪(n2logn). Once the distance and ranking are completed, computing one local generalized correlation requires 𝒪(n2) (see Algorithm C4). Thus, a naive approach to compute all local generalized correlations requires at least 𝒪(n2max{n2,p,q}) by going through all possible scales, meaning possibly 𝒪(n4) which would be computationally prohibitive. However, given the distance and ranking information, we devised an algorithm that iteratively computes all local correlations in 𝒪(n2) by re-using adjacent smaller local generalized correlations (see Algorithm C5). Therefore, when including the distance computation and ranking overheads, the MGC statistic is computed in 𝒪(n2max{logn,p,q})), which has the same running time as the Hhg statistic, and the same running time up to a factor of logn as global correlations like Dcorr and Mcorr, which require 𝒪(n2max{p,q}) time. By utilizing a multi-core architecture, Mgc can be computed in 𝒪(n2max{logn,p,q}/T) instead. As T=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 n2 term, so computing the p-value requires 𝒪(n2max{logn,p,q,r}/T).

Mgc algorithms and testing procedures

Request a detailed protocol

Six 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 p-value 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 𝒪(n2max(logn,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 (xi,yi) pairs, an integer r for the number of random permutations.
Output: (i) MGC statistic c*, (ii) the optimal scale (k,l), (iii) the p-value p(c),
     function MG((xi,yi), for i[n])
     (1) Calculate all pairwise distances:
for i,j:=1,,n do
aij=δx(xi,xj)dx is the distance between pairs of x samples
bij=δy(yi,yj)dy is the distance between pairs of y samples
end for
Let A={aij} and B={bij}.
     (2) Calculate Multiscale Correlation Map 𝒞 & Mgc Test Statistic:
[c,𝒞,k,l]=MGCSAMPLESTAT(A,B)Algorithm C2
     (3) Calculate the p-value
pval(c)=PERMUTATIONTEST(A,B,r,c)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 𝒪(n2).
Input: A pair of distance matrices (A,B)Rn×n×Rn×n.
Output: The Mgc statistic cR, all local statistics 𝒞Rn×n, and the corresponding local scale (k,l)N×N.
1:function MGCSampleStat(A,B)
2: 𝒞=MGCALLLOCAL(A,B)All local correlations
3: τ=THRESHOLDING(𝒞)find a threshold to determine large local correlations
4: for i,j:=1,,ndo rijI(cij>τ)end for identify all scales with large correlation
5: {rij:i,j=1,,n}binary map encoding scales with large correlation
6: =CONNECTED()largest connected component of the binary matrix
7: c𝒞(n,n)use the global correlation by default
8: kn,ln
9: if (i,jrij)2n thenproceed when the significant region is sufficiently large
10: [c,k,l]max(𝒞)find the smoothed maximum and the respective scale
11: end if
12:end Function
Input: 𝒞Rn×n.
Output: A threshold t to identify large correlations.
13:function Thresholding 𝒞
14: τcij<0(cij)2/cij<01variance of all negative local generalized correlations
15: τmax{0.01,τ}×3.5threshold based on negative correlations
16: τmax{τ,2/n,cnn}
17:end Function
Pseudocode C3 Permutation Test. This algorithm uses the random permutation test with r random permutations for the p-value, requiring 𝒪(rn2logn) for Mgc. In the real-data experiment, we always set r=10,000. Note that the p-value computation for any other global generalized correlation coefficient follows from the same algorithm by replacing Mgc with the respective test statistic.
Input: A pair of distance matrices (A,B)Rn×n×Rn×n, the number of permutations r, and Mgc statistic c* for the observed data.
Output: The p-value pval[0,1].
1:function PermutationTest(A, B, r, c*)
2: for t:=1,,r do
3: π=RANDPERM(n)generate a random permutation of size n
4: c0[t]=MGCSAMPLESTAT(A,B(π,π))calculate the permuted Mgc statistic
5: end for
6: pval(c)1tt=1rI(cc0[t])compute p-value of Mgc
7:end function
Pseudocode C4 Compute local test statistic at a given scale. This algorithm runs in 𝒪(n2) once the rank information is provided, which is suitable for Mgc computation if an optimal scale is already estimated. But it would take 𝒪(n4) if used to compute all local generalized correlations. Note that for the default Mgc implementation uses single centering, the centering function centers A by column and B by row, and the sorting function sorts A within column and B within row. By utilizing T=log(n) cores, the sorting function can be easily parallelized to take 𝒪(n2log(n)/T)=𝒪(n2).
Input: A pair of distance matrices (A,B)Rn×n×Rn×n, and a local scale (k,l)N×N.
Output: The local generalized correlation coefficient ckl[1,1].
1:function LocalGenCorr(A, B, k, l)
2:  for Z:=A,B do Z=SORT(Z) end for parallelized sorting
3: for Z:=A,B do Z=CENTER(Z) end forcenter distance matrices
4: c~kltr((AA)T×(B(B)T))un-normalized local distance covariance
5: vAtr((AA)T×(A(A)T))local distance variances
6: vBtr((BB)T×(B(B)T))
7: eAi,j=1n(AA)ijsample means
8: eBi,j=1n(BB)ij
9: ckl(c~kleAeB/n2)/(vA(eA/n)2)(vB(eB/n)2)center and normalize
10:end function
Pseudocode C5 Compute the multiscale correlation map (i.e., all local generalized correlations) in 𝒪(n2logn/T). Once the distances are sorted, the remaining algorithm runs in 𝒪(n2). An important observation is that each product aijbij is included in ckl if and only if (k,l) satisfies kR(Aj,i) and lR(Bj,i), so it suffices to iterate through aijbij for i,j:=1,,n, and add the product simultaneously to all ckl whose scales are no more than (R(Aj,i),R(Bj,i)). To achieve the above, we iterate through each product, add it to ckl at (kl)=(R(Aj,i),R(Bj,i)) only (so only one local scale is accessed for each operation); then add up adjacent ckl for k,l=1,,n. The same applies to all local covariances, variances, and expectations.
Input: A pair of distance matrices (A,B)Rn×n×Rn×n.
Output: The multiscale correlation map 𝒞[1,1]n×nfor k,l=1,,n.
1:function MGCAllLocal(A, B)
2:  for Z:=A,B do Z=SORT(Z) end for
3: for Z:=A,B do Z=CENTER(Z)end for
4: for i,j:=1,,n doiterate through all local scales to calculate each term
5: kijZ
6: lijZ
7: c~klc~kl+aijbij
8: vkAvkA+aij2
9: vlBvlB+bij2
10: ekAekA+aij
11: elBelB+bij
12: end for
13: for k:=1,,n-1 doiterate through each scale again and add up adjacent terms
14: c~1,k+1c~1,k+c~1,k+1
15: c~k+1,1c~k+1,1+c~k+1,1
16:  for Z:=A,B do vk+1ZvkZ+vk+1Z end for
17:  for Z:=A,B do ek+1ZekZ+ek+1Z end for
18: end for
19:  for k,l:=1,,n-1 do
20: c~k+1,l+1c~k+1,l+c~k,l+1+c~k+1,l+1c~k,l
21: end for
22: for k,l:=1,,n do
23: ckl(c~klekAelB/n2)/(vkAekA2/n2)(vlBelB2/n2)
24: end for
25:end function
Pseudocode C6 Power computation of Mgc against a given distribution. By repeatedly sampling from the joint distribution FXY, sample data of size n under the null and the alternative are generated for r Monte-Carlo replicates. The power of Mgc follows by computing the test statistic under the null and the alternative using Algorithm 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 FXY, the sample size n, the number of MC replicates r, and the type 1 error level a.
Output: The power ß of Mgc.
1:function MGCPower(FXY, n, r, a)
2:  for t:=1,,r do
3: for i:=[n] do
4: xi0iidFX,yi0iidFYsample from null
5: (xi1,yi1)iidFXY,sample from alternative
6: end for
7: for i,j:=1,,n do
8: aij0=δx(xi0,xj0), bij0=δy(yi0,yj0)pairwise distances under the null
9: aij1=δx(xi1,xj1), bij1=δy(yi1,yj1)pairwise distances under the alternative
10: end for
11: c0[t]=MGCSAMPLESTAT(A0,B0)Mgc statistic under the null
12: c1[t]=MGCSAMPLESTAT(A1,B1)Mgc statistic under the alternative
13: end for
14: ωαCDF1α(c0[t],t[r])the critical value of Mgc under the null
15: βt=1r(c1[t]>ωα)/rcompute power by the alternative distribution
16:end function

Simulation dependence functions

Request a detailed protocol

This section provides the 20 different dependency functions used in the simulations. We used essentially the exact same relationships as previous publications to ensure a fair comparison (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 xRp, we denote x[d],d=1,,p as the dth dimension of the vector x. For the purpose of high-dimensional simulations, wRp is a decaying vector with w[d]=1/d for each d, such that wTx is a weighted summation of all dimensions of x. Furthermore, 𝒰(a,b) denotes the uniform distribution on the interval (a,b), (p) denotes the Bernoulli distribution with probability p, 𝒩(μ,Σ) denotes the normal distribution with mean µ and covariance S, U and V represent some auxiliary random variables, κ is a scalar constant to control the noise level (which equals 1 for one-dimensional simulations and 0 otherwise), and ϵ is a white noise from independent standard normal distribution unless mentioned otherwise.

For all the below equations, (X,Y)iidFXY=FY|XFX. For each relationship, we provide the space of (X,Y), and define FY|X and FX, as well as any additional auxiliary distributions.

1. Linear (X,Y)Rp×R,


2. Exponential (X,Y)Rp×R:


3. Cubic (X,Y)Rp×R:


4. Joint normal (X,Y)Rp×Rp: Let ρ=1/2p, Ip be the identity matrix of size p×p, Jp be the matrix of ones of size p×p, and Σ=[IpρJpρJp(1+0.5κ)Ip]. Then


5. Step Function (X,Y)Rp×R


where I is the indicator function, that is I(z) is unity whenever z true, and zero otherwise.

6. Quadratic (X,Y)Rp×R:


7. W Shape (X,Y)Rp×R:U𝒰(1,1)p,


8. Spiral (X,Y)Rp×R:U𝒰(0,5), ϵ𝒩(0,1)


9. Uncorrelated Bernoulli (X,Y)Rp×R:U(0.5)ϵ1𝒩(0,Ip),ϵ2𝒩(0,1),


10. Logarithmic (X,Y)Rp×Rp:ϵ𝒩(0,Ip)


11. Fourth Root (X,Y)Rp×Rp:


12. Sine Period 4π(X,Y)Rp×Rp:U𝒰(1,1),V𝒩(0,1)p,θ=4π,


13. Sine Period 16π(X,Y)Rp×Rp: Same as above except θ=16π and the noise on Y is changed to 0.5κϵ.

14. Square (X,Y)Rp×Rp: Let U𝒰(1,1),V𝒰(1,1),ϵ𝒩(0,1)p,θ=π8. Then


for d=1,,p.

15. Two Parabolas (X,Y)Rp×R: ϵ𝒰(0,1),U(0.5),


16. Circle (X,Y)Rp×R:U𝒰(1,1)p,ϵ𝒩(0,Ip),r=1,


17. Ellipse (X,Y)Rp×R: Same as above except r=5.

18. Diamond (X,Y)Rp×Rp: Same as 'Square' except θ=π4.

19. Multiplicative Noise (x,y)Rp×Rp:u𝒩(0,Ip),


20. Multimodal Independence (X,Y)Rp×Rp:LetU𝒩(0,Ip),V𝒩(0,Ip),U(0.5)p,V(0.5)p. Then


For each distribution, X and Y are dependent except (20); for some relationships (8,14,16-18) they are independent upon conditioning on the respective auxiliary variables, while for others they are 'directly' dependent. A visualization of each dependency with D=Dy=1 is shown in Figure 2—figure supplement 1.

For the increasing dimension simulation in the main paper, we always set κ=0 and n=100, with p increasing. Note that q=p for types 4, 10, 12, 13, 14, 18, 19, 20,, otherwise q=1. The decaying vector w is utilized for p>1 to make the high-dimensional relationships more difficult (otherwise, additional dimensions only add more signal). For the one-dimensional simulations, we always set p=q=1, κ=1 and n=100.

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 Inventory-Revised the characterized personality along five dimensions (Costa and McCrae, 1992).

This dataset consists of 42 subjects, each with 197 time-steps of resting-state functional magnetic resonance activity (rs-fMRI) activity, as well as the subject’s five-dimensional 'personality'. Adelstein et al. (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 five-factor personality modality, we used the Euclidean distance. For the brain activity modality, we derived the following comparison function. For each scan, (i) run Configurable Pipeline for the Analysis of Connectomes pipeline (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 Kullback-Leibler divergence across regions to obtain a similarity matrix across comparing all regions. Then, use the normalized Hellinger distance to compute distances between each subject.

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 semi-parametric 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 follow-up experiment.

In this study, we performed the following five sets of tests on those data:

  1. Ovarian vs. normal for all proteins,

  2. Ovarian vs. normal for each individual protein,

  3. Pancreas vs. normal for all proteins,

  4. Pancreas vs. all others for each individual protein,

  5. Pancreas vs. normal for each individual protein.

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 p-values; we used the Benjamini-Hochberg 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
Results for cancer peptide screening.

The first two rows report the p-values for the tests of interest based on all peptides. The next four rows report the number of significant proteins from individual peptide tests; the Benjamini-Hochberg procedure is used to locate the significant peptides by controlling the false discovery rate at 0.05.
Testing pairs / MethodsMgcMantelDcorrMcorrHhg
1Ovar vs. Norm: p-value0.00010.00010.00010.00010.0001
2Ovar vs. Norm: # peptides218190186178225
3Pancr vs. Norm: p-value0.00820.06850.06690.01920.0328
4Panc vs. Norm: # peptides976711
5Panc vs. All: # peptides10003
6# peptides unique to Panc10002
7# false positives for Panc0n/an/an/a1
Appendix 1—table 1—source data 1

Ovarian testing results.
Appendix 1—table 1—source data 2

Pancreatic testing results.
Appendix 1—table 1—source data 3

Peptide screening results for pancreatic.

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 p-value ≈ 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 p-value (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 p-value for each peptide.

Appendix 1—table 2
For each of Mgc, Dcorr, Mcorr, Hhg, Hsic, Mantel, Pearson, and Mic, list the top four peptides identified for Panc vs All and the respective corrected p-value using Benjamini-Hochberg.

Bold indicates a significant peptide at type 1 error level 0.05. The top candidates are very much alike except Mic. In particular, neurogranin is consistently among the top candidates for all methods, but is only significant while using Mgc, Hsic, and Hhg; there are two other significant proteins from Hsic and Hhg, but they do not further improve the classification performance comparing to just using neurogranin. Note that the p-values from Mantel and Pearson are always 1 after Benjamini-Hochberg correction, so their respective top peptides are identified using raw p-values without correction.
methodTop four identified peptides
Mgcneurograninfibrinogen protein 1tropomyosin alpha-3ras suppressor protein 1
Dcorrneurograninfibrinogen protein 1kinase 6twinfilin-2
Mcorrneurograninfibrinogen protein 1kinase 6tropomyosin alpha-3
Hsicneurogranintropomyosin alpha-3kinase 6tripeptidyl-peptidase 2
Hhgneurograninfibrinogen protein 1tropomyosin alpha-3platelet basic protein
Mantelneurograninadenylyl cyclasetropomyosin alpha-3alpha-actinin-1
Pearsonneurograninadenylyl cyclasetropomyosin alpha-3alpha-actinin-1
Mickinase BS100-A9ERF3Athymidine

All three methods reveal the same unique protein for pancreas: neurogranin. Hsic identifies another peptide (tropomyosin alpha-3 chain isoform 4), and Hhg identifies a third peptide (fibrinogen-like protein 1 precursor). However, fibrinogen-like protein 1 precursor is not significant for p-value 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 p-values 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 time-varying stimuli? Any region that is selected as significant is a false positive by construction. By testing each brain region separately, Mgc provides a distribution of false positive rates. If Mgc is valid, the resulting distribution should be centered around the significance level, which is set at 0.05 for these experiments.

We considered 25 resting state fMRI experiments from the 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 pre-processing 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.

Appendix 1—figure 1
We demonstrate that Mgc is a valid test that does not inflate the false positives in screening and variable selection.

This figure shows the density estimate for the false positive rates of applying Mgc to select the 'falsely significant' brain regions versus independent noise experiments; dots indicate the false positive rate of each experiment. The mean ± standard deviation is 0.0538 ± 0.0394.

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 I7-6850K 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(n3) and thus significantly slower.

Appendix 1—table 3
The actual testing time (in seconds) on real data.

Data availability

To facilitate reproducibility, we make all datasets available from: (copy archived at


  1. Book
    1. Bracewell RN
    2. Bracewell RN
    The Fourier Transform and Its Applications
    New York: McGraw-Hill.
    1. Coifman RR
    2. Maggioni M
    (2006) Diffusion wavelets
    Applied and Computational Harmonic Analysis 21:53–94.
  2. Book
    1. Costa PT
    2. McCrae RR
    Neo PI-R Professional Manual, 396
     PAR Psychological Assessment R.
  3. Conference
    1. Da Zheng DM
    2. Burns R
    3. Vogelstein JT
    4. Priebe CE
    5. Szalay AS
    FlashGraph: processing Billion-Node graphs on an array of commodity SSDs
    USENIX Conference on File and Storage Technologies.
  4. Book
    1. Daubechies I
    Ten Lectures on Wavelets
  5. Book
    1. Edelsbrunner H
    2. Harer JL
    Computational Topology: An Introduction
    American Mathematical Society.
  6. Book
    1. Good P
    Permutation, Parametric, and Bootstrap Tests of Hypotheses
    1. Goodfellow I
    2. Pouget-Abadie J
    3. Mirza M
    4. Xu B
    5. Warde-Farley D
    6. Ozair S
    7. Courville A
    8. Bengio Y
    Advances in Neural Information Processing System
    2672–2680, Generative adversarial nets, Advances in Neural Information Processing System, MIT Press.
  7. Report
    1. Gorfine M
    2. Heller R
    3. Heller Y
    Comment on Detecting Novel Associations in Large Data Sets
    Israel Institute of Technology.
    1. Gretton A
    2. Borgwardt KM
    3. Rasch M
    4. Schölkopf B
    5. Smola AJ
    Advances in Neural Information Processing Systems
    513–520, A kernel method for the two-sample-problem, Advances in Neural Information Processing Systems, MIT Press.
    1. Gretton A
    2. Gyorfi L
    Consistent nonparametric tests of independence
    Journal of Machine Learning Research 11:1391–1423.
  8. Book
    1. Hastie T
    2. Tibshirani R
    3. Friedman JH
    Elements of Statistical Learning
    New York: Springer.
    1. Heller R
    2. Heller Y
    3. Kaufman S
    4. Brill B
    5. Gorfine M
    Consistent distribution-free -sample and independence tests for univariate random variables
    Journal of Machine Learning Research 17:1–54.
  9. Book
    1. Kendall MG
    Rank Correlation Methods
    London: Griffin.
  10. Book
    1. Lee JA
    2. Verleysen M
    Nonlinear Dimensionality Reduction
    Springer Science & Business Media.
    1. Levina E
    2. Bickel PJ
    Advances in Neural Information Processing Systems
    Maximum likelihood estimation of intrinsic dimension, Advances in Neural Information Processing Systems,  MIT Press.
    1. Mantel N
    The detection of disease clustering and a generalized regression approach
    Cancer Research 27:209–220.
  11. Book
    1. Pearl J
    Causality: Models, Reasoning, and Inference
    Cambridge University Press.
    1. Rényi A
    (1959) On measures of dependence
    Acta Mathematica Academiae Scientiarum Hungaricae 10:441–451.
    1. Rizzo ML
    2. Székely GJ
    (2016) Energy distance
    Wiley Interdisciplinary Reviews: Computational Statistics 8:27–38.
  12. Conference
    1. Roncal WG
    2. Koterba ZH
    3. Mhembere D
    4. Kleissas DM
    5. Vogelstein JT
    6. Burns R
    7. Bowles AR
    8. Donavos DK
    9. Ryman S
    10. Jung RE
    11. Wu L
    12. Calhoun VD
    13. Jacob Vogelstein R
    MIGRAINE: mri graph reliability analysis and inference for connectomics
    Global Conference on Signal and Information Processing.
  13. Conference
    1. Sutherland DJ
    2. Tung H-Y
    3. Strathmann H
    4. De S
    5. Ramdas A
    6. Smola A
    7. Gretton A
    Generative models and model criticism via optimized maximum mean discrepancy
    International Conference on Learning Representations.
    1. Szekely GJ
    2. Rizzo ML
    Testing for equal distributions in high dimension
    InterStat 10:.
    1. Xing EP
    2. Ng AY
    3. Jordan MI
    4. Russell S
    Distance metric learning with application to clustering with side-information
    Advances in Neural Information Processing Systems 15:505–512.
    1. Zhang Z
    2. Wang J
    3. Zha H
    (2012) Adaptive manifold learning
    IEEE Transactions on Pattern Analysis and Machine Intelligence 34:253–265.

Article and author information

Author details

  1. Joshua T Vogelstein

    1. Johns Hopkins University, Baltimore, United States
    2. Child Mind Institute, New York, United States
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2487-6237
  2. Eric W Bridgeford

    Johns Hopkins University, Baltimore, United States
    Competing interests
    No competing interests declared
  3. Qing Wang

    Johns Hopkins University, Baltimore, United States
    Data curation
    Competing interests
    No competing interests declared
  4. Carey E Priebe

    Johns Hopkins University, Baltimore, United States
    Competing interests
    No competing interests declared
  5. Mauro Maggioni

    Johns Hopkins University, Baltimore, United States
    Competing interests
    No competing interests declared
  6. Cencheng Shen

    University of Delaware, Delaware, United States
    Conceptualization, Formal analysis, Investigation, Visualization, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1030-1432


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.


This work was partially supported by the Child Mind Institute Endeavor Scientist Program, the National Science Foundation award DMS-1712947, 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 N66001-15-C-4041, the XDATA program of DARPA administered through Air Force Research Laboratory contract FA8750-12-2-0303, DARPA Lifelong Learning Machines program through contract FA8650-18-2-7834, the Office of Naval Research contract N00014-12-1-0601, the Air Force Office of Scientific Research contract FA9550-14-1-0033. 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

  1. Received: September 3, 2018
  2. Accepted: January 14, 2019
  3. Accepted Manuscript published: January 15, 2019 (version 1)
  4. Version of Record published: February 22, 2019 (version 2)


© 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.


  • 3,467
    Page views
  • 417
  • 13

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Joshua T Vogelstein
  2. Eric W Bridgeford
  3. Qing Wang
  4. Carey E Priebe
  5. Mauro Maggioni
  6. Cencheng Shen
Discovering and deciphering relationships across disparate data modalities
eLife 8:e41690.

Share this article

Further reading

    1. Computational and Systems Biology
    2. Neuroscience
    Domingos Leite de Castro, Miguel Aroso ... Paulo Aguiar
    Research Article Updated

    Closed-loop neuronal stimulation has a strong therapeutic potential for neurological disorders such as Parkinson’s disease. However, at the moment, standard stimulation protocols rely on continuous open-loop stimulation and the design of adaptive controllers is an active field of research. Delayed feedback control (DFC), a popular method used to control chaotic systems, has been proposed as a closed-loop technique for desynchronisation of neuronal populations but, so far, was only tested in computational studies. We implement DFC for the first time in neuronal populations and access its efficacy in disrupting unwanted neuronal oscillations. To analyse in detail the performance of this activity control algorithm, we used specialised in vitro platforms with high spatiotemporal monitoring/stimulating capabilities. We show that the conventional DFC in fact worsens the neuronal population oscillatory behaviour, which was never reported before. Conversely, we present an improved control algorithm, adaptive DFC (aDFC), which monitors the ongoing oscillation periodicity and self-tunes accordingly. aDFC effectively disrupts collective neuronal oscillations restoring a more physiological state. Overall, these results support aDFC as a better candidate for therapeutic closed-loop brain stimulation.

    1. Cancer Biology
    2. Computational and Systems Biology
    Sara Latini, Veronica Venafra ... Francesca Sacco
    Research Article

    Currently, the identification of patient-specific therapies in cancer is mainly informed by personalized genomic analysis. In the setting of acute myeloid leukemia (AML), patient-drug treatment matching fails in a subset of patients harboring atypical internal tandem duplications (ITDs) in the tyrosine kinase domain of the FLT3 gene. To address this unmet medical need, here we develop a systems-based strategy that integrates multiparametric analysis of crucial signaling pathways, and patient-specific genomic and transcriptomic data with a prior knowledge signaling network using a Boolean-based formalism. By this approach, we derive personalized predictive models describing the signaling landscape of AML FLT3-ITD positive cell lines and patients. These models enable us to derive mechanistic insight into drug resistance mechanisms and suggest novel opportunities for combinatorial treatments. Interestingly, our analysis reveals that the JNK kinase pathway plays a crucial role in the tyrosine kinase inhibitor response of FLT3-ITD cells through cell cycle regulation. Finally, our work shows that patient-specific logic models have the potential to inform precision medicine approaches.