# Correlation networks: Interdisciplinary approaches beyond thresholding

Naoki Masuda1,2, Zachary M. Boyd3, Diego Garlaschelli4,5, and Peter J. Mucha6
1Department of Mathematics, State University of New York at Buffalo
2Institute for Artificial Intelligence and Data Science, State University of
New York at Buffalo
3Department of Mathematics, Brigham Young University
4Lorentz Institute for Theoretical Physics, Leiden University, The Netherlands
5IMT School of Advanced Studies, Lucca, Italy
6Department of Mathematics, Dartmouth College

###### Abstract

Many empirical networks originate from correlational data, arising in domains
as diverse as psychology, neuroscience, genomics, microbiology, finance, and
climate science. Specialized algorithms and theory have been developed in
different application domains for working with such networks, as well as in
statistics, network science, and computer science, often with limited
communication between practitioners in different fields. This leaves
significant room for cross-pollination across disciplines. A central challenge
is that it is not always clear how to best transform correlation matrix data
into networks for the application at hand, and probably the most widespread
method, i.e., thresholding on the correlation value to create either
unweighted or weighted networks, suffers from multiple problems. In this
article, we review various methods of constructing and analyzing correlation
networks, ranging from thresholding and its improvements to weighted networks,
regularization, dynamic correlation networks, threshold-free approaches, and
more. Finally, we propose and discuss a variety of key open questions
currently confronting this field.

##  1 Introduction

Correlation matrices capture pairwise similarity of multiple, often temporally
evolving signals, and are used to describe system interactions in various
diverse disciplines of science and society, from financial economics to
psychology, bioinformatics, neuroscience, and climate science, to name a few.
Correlation analysis is often a first step in trying to understand complex
systems data [1]. Existing methods for analyzing correlation matrix data are
abundant. Very well established methods include principal component analysis
(PCA) [2] and factor analysis (FA) [3, 4], which can yield a small number of
interpretable components from correlation matrices, such as a global market
trend when applied to stock market data, or spatio-temporal patterns of air
pressure when applied to atmospheric data. Another major method for analyzing
correlation matrix data is the Markowitz’s portfolio theory in mathematical
finance, which aims to minimize the variance of financial returns while
keeping the expected return above a given threshold [5, 6]. In a related vein,
random matrix theory has been a key theoretical tool for analyzing economic
and other correlation matrix data for a couple of decades [6]. Various new
methods for analyzing correlation matrix data have also been proposed.
Examples include detrended cross-correlation analysis [7, 8, 9], correlation
dependency, defined as the difference between the partial correlation
coefficient and the Pearson correlation coefficient given three nodes [10,
11], determination of optimal paths between distant locations in correlation
matrix data [12], early warning signals for anticipating abrupt changes in
multidimensional dynamical systems including the case of networked systems
[13, 14, 15], and energy landscape analysis for multivariate time series data
particularly employed in neuroscience [16, 17].

The last two decades have also seen successful applications of tools from
network science and graph theory to correlational data. A correlation matrix
can be mapped onto a network, which we refer to here as a _correlation
network_ , where nodes represent elements and edges are informed by the
strength of the correlation between pairs of elements. Correlation network
analysis generally intends to extract useful information from data, such as
the patterns of interactions among nodes or a ranking of nodes. An ideal
correlation network analysis appropriately adapts concepts and methods
developed in network science to the case of correlation networks, generating
knowledge that standard methods for correlation matrices (such as PCA) do not
produce. Although correlation does not necessarily reflect a physical
connection or direct interaction between two nodes, correlation matrices are
conventionally used as a relatively inexpensive substitute of such direct
connections, whose data are often less available than correlation matrix data.
Correlation networks are also useful for visualization [18]. Correlation
network analysis has been used in various research disciplines, typically not
much behind wherever correlation matrix analysis is used, as we will review in
section 2. In our survey here, we focus on correlation networks, with an
emphasis on identifying different methods used to transform correlation
matrices into correlation networks.

The validity of correlation network analysis remains an outstanding question,
especially because the decisions about how to best construct network
representations from correlation matrices is far from straightforward. One of
the simplest methods is to threshold on the correlation value measured for
each pair of nodes (see section 3.2). However, while such a simple
thresholding is widely used, it introduces various problems. These problems
have led to proposals of alternative methods for generating correlation
networks, which we will cover in sections 3.3–3.6.

Before proceeding, we raise some important clarifications. First, correlation
networks as we consider here are different from network architectures that
exploit correlation in data. For example, the Progressive Spatio-Temporal
Correlation Network (PSCC-net) is an algorithm to detect and localize
manipulations in the input image data by taking advantage of spatial
correlation structure in images [19]. The superpixel group-correlation network
(SGCN) [20] and the deep correlation network (DCNet) [21] are encoder-decoder
and deep-learning network architectures, respectively, for salient object
detection in images. These “networks” are in the sense of neural network
architecture in artificial intelligence and machine learning, whereas here we
consider “networks” to denote graphs in network science.

Second, we focus on correlation networks based on the Pearson correlation
coefficient or its close variants such as the partial correlation coefficient.
In fact, there are numerous other definitions for quantifying the similarity
between data obtained from node pairs [22, 23, 24, 25, 26, 27, 28]. Examples
include similarity networks whose edges are determined using the rank
correlation coefficient [29, 30], the mutual information [31, 32, 33], and
partial mutual information [34, 35]. However, a majority of concepts and
techniques explained in our main technical section, section 3, such as the
detrimental effect of thresholding and dichotomizing the edge weight, use of
weighted networks, graphical lasso, and importance of null model, also hold
true when one constructs correlation networks using these or other alternative
methods.

Third, we do not discuss causal inference in the present paper. The inherently
symmetric nature of correlation matrices or general similarity matrices mean
that these matrices or networks in principle do not inform us of causality or
directionality between nodes. If the original data upon which one calculates
correlation matrices and networks are multivariate time series, a plethora of
methods are available for inferring causality between nodes and associated
directed networks. For these techniques, see, e.g., [36, 28, 37, 18]. In a
related vein, we do not discuss time-lagged correlation in this paper, since
these are also asymmetric in general, although many of the same considerations
we raise here also apply to lagged correlations.

##  2 Application areas

Correlation network analysis is common in many research fields. In this
section, we survey typical correlation networks and their analysis in
representative research fields.

###  2.1 Psychological networks

There are various multivariate psychological data, from which one can
construct networks [38, 18]. For example, in personality research, researchers
construct personality networks in which each node can be a personal trait or
goal such as being organized, being lazy, and wanting to stay safe. Edges
between a pair of nodes typically represent a type of conditional association
between the two nodes. Samples are frequently participants in the research
responding to various questionnaires on a numeric scale (e.g., 5-point scale
ranging from 1: strongly disagree to 5: strongly agree) corresponding to
nodes. From a cross-sectional data set, one can calculate (conditional)
correlation between pairs of nodes. Researchers are also increasingly
combining surveys with alternative data collection modalities, for example,
sensor data for daily movement or neural markers of stress [39, 40]

Another major type of psychological network is symptom networks employed in
mental health research. Symptoms of a psychological, including
psychopathological, condition, such as major depression and schizophrenia, are
interrelated. Furthermore, causality between symptoms such as fatigue,
headaches, concentration problems, and insomnia, and a psychopathological
disorder, is often unclear. It has been suggested that a disorder does not
originate from a single root cause, which motivates the study of symptom
networks [41, 42, 43, 18]. Nodes in a symptom network are symptoms, and one
can use association between pairs of symptoms calculated from the
participants’ responses to define edges. Analysis of symptom networks may help
us to predict how an individual develops psychopathology in the future,
understand comorbidity as strong connection between symptoms of two mental
disorders, and propose central nodes as possible targets of intervention [43].
Health-promoting behaviors can also be treated as nodes in these networks to
suggest key behavioral intervention points [44].

Panel data, i.e., longitudinal measurements of variables from samples, are
increasingly common for network approaches [18]. In this case, one obtains
correlation networks at multiple time points. Then, one can construct time-
varying correlation networks (see section 3.7) or within-person correlation
networks [45] that reflect temporal symptom patterns and ideally expose
individual differences and possible causal pathways in mental health patterns
related to disorders [46, 47]. However, the validity of psychological network
approach should be further studied. Research has shown that symptom networks
have poor reproducibility across samples, likely due to measurement error in
assessing symptoms among other reasons [42, 48].

###  2.2 Brain networks

Various notions of brain connectivity have been essential to better
understanding different neural functions. Studies of such brain networks
constitute a major part of a research field that is often referred to as
_network neuroscience_ [49, 50]. (See also the related material about network
representations in [51].) Multivariate time series of neuronal signals
recorded from the brain are a major source of data used in network
neuroscience research. Such data may be recorded in a resting state or when
participants are performing some task. Functional networks or functional
connectivity refer to correlation-based networks constructed from multivariate
neuronal time series data, obtained through, e.g., neuroimaging or
electroencephalography, where the term “functional” in this setting
effectively means correlational. In the case of multivariate time series data,
there are various other methods to infer directed brain networks, which is
referred to as effective connectivity, but we do not cover directed networks
in this article. Brain networks based on anatomical connectivity between brain
regions are referred to as structural networks. Functional connectivity, or an
correlation-based edge between two nodes in the brain does not imply the
presence of an edge between the same pair of nodes in the structural network.
A typical node in the functional and structural brain networks is either a
voxel (i.e., cube of side length of, e.g., 1 mm) or a spherical region of
interest (ROI), which is a sphere in the brain. See [49, 25, 52, 27] for
review of brain networks, including functional networks.

The most typical functional neuronal networks come from neuroimaging data, in
particular functional magnetic resonance imaging (fMRI) data, which are
measured using blood-oxygeneration-level-dependent (BOLD) imaging techniques
[53]. Functional connectivity between voxels or between spherical ROIs, or
other types of nodes, is calculated by a correlation between fMRI time signals
at the two nodes after one has bandpassed the fMRI time series at each node to
remove artifacts, with a frequency band of, e.g., 0.01-0.1 Hz. Functional MRI
improves on electroencephalogram (EEG) and magnetoencephalography (MEG) in
spatial resolution at the expense of temporal resolution, but functional EEG
and MEG networks are not uncommon. We also note that EEG and MEG signals are
oscillatory, so one has to calculate the functional connectivity between each
pair of nodes using methods that are aware of the oscillatory nature of the
signal, such as using phase lag index or amplitude envelope correlation [54]
rather than conventional correlation coefficients or mutual information.

Structural covariance networks are another type of correlation brain network
where the edges are defined as the correlation/covariance of the corrected
gray matter volume or cortical thickness between brain regions $i$ and $j$,
where the samples are participants [55, 56]. Morphometric similarity networks
are a variant of structural covariance networks. In morphometric similarity
networks, one uses various morphometric variables, not just a single one such
as cortical thickness, for each node (i.e., ROI) [57]. One calculates the
correlation between two nodes by regarding each morphometric variable as a
sample. Therefore, differently from structural covariance networks based on
cortical thickness, one can calculate a correlation network for each
individual.

In neuroreceptor similarity networks, an edge between two nodes, or ROIs, is
the correlation in terms of receptor density [58]. Specifically, one first
calculates a vector of neurotransmitter density for each ROI, with each entry
of the vector corresponding to one type of receptor. Then, one computes the
correlation between each pair of ROIs, called receptor similarity.

###  2.3 Gene co-expression networks

Genes do not work in isolation. Gene co-expression networks have been useful
for figuring out webs of interaction among genes using network analysis
methods [59, 60, 61, 62, 63, 28, 64, 65, 66]. They are a type of data in a
subfield of network science often referred to as _network biology_ or _network
medicine_. Gene co-expression networks are correlation networks in the
generalized sense considered here, including the case of other measures of
similarity. A typical measurement is the amount of gene expression for
different genes and samples, where a sample most commonly corresponds to a
human or animal individual. If one measures the expression of various genes
for the same set of samples, we can calculate the co-expression between each
pair of genes by calculating the sample correlation, yielding a correlation
matrix. Depending on the questions being asked in the study, it may be
important to calculate the underlying correlations with different factors to
account for the effects of heterogeneous gene frequencies [67, 68]. It is
common to transform a correlation matrix into a network and then apply various
network analysis methods, for example community detection with the aim of
estimating the group of genes that are associated with the same phenotype111A
_phenotype_ is a set of observable traits of an organism and is usually
contrasted with the underlying _genotype_ that causes (or influences) the
phenotype. such as a disease. In this manner, correlation network analysis has
been a useful tool for gene screening, which can lead to identification of
biomarkers and therapeutic targets. In addition to community detection,
identifying hub genes in co-expression networks helps finding key genes, for
example, for cancer [69].

Different ways of defining co-expression matrices and networks from gene
expression data include tissue-to-tissue co-expression (TTC) networks [70]
(also see [71, 72]). A TTC network proposed in [70] is a bipartite network,
and its node is a gene-tissue pair. An edge between two nodes, denoted by
$(\tilde{g}_{i},\tilde{t}_{i})$ and $(\tilde{g}_{j},\tilde{t}_{j})$, where
$\tilde{g}_{i}$ and $\tilde{g}_{j}$ are genes, and $\tilde{t}_{i}$ and
$\tilde{t}_{j}$ are tissues, represents the sample correlation as in
conventional co-expression networks. However, by definition, the correlation
is calculated only between node pairs belonging to different tissues, i.e.,
only for $\tilde{t}_{i}\neq\tilde{t}_{j}$. Therefore, TTC networks
characterize co-expression of genes across different tissues.

Co-expression of genes $i$ and $j$ implies that $i$ and $j$ are both expressed
at a high level in some samples (usually individuals) and both expressed at a
low level in other individuals. Co-expressed genes tend to be involved in
common biological functions. There are in fact multiple biophysical and non-
biophysical reasons for gene co-expression [63]. For example, a transcription
factor, a protein that binds to DNA, may regulate different genes $i$ and $j$
that are physically close on a chromosome. If this is the case, differential
levels of regulation by the transcription factor across individuals can create
co-expression of $i$ and $j$. Another mechanism of co-expression is that the
expression of genes $i$ and $j$, which may be located far from each other on
the chromosome or on different chromosomes, may depend on the temperature.
Then, $i$ and $j$ would be co-expressed if different individuals are sampled
from living environments with different temperatures. Variation in ages of the
individuals can similarly create co-expression among age-related genes.
Alternatively, co-expression may originate from non-biological sources, such
as technical or laboratory ones, whose exact origins are often unknown.

One is often interested in looking for differential co-expression, which
refers to the different levels of gene co-expression between two
phenotypically different sets of samples, such as a disease set versus a
control set, or in two types of tissues [63, 65]. Differential co-expression
often reveals information that one cannot obtain by examining differential
expression (as opposed to co-expression), i.e., different levels of gene
expression between the two sets of samples [73].

###  2.4 Metabolite networks

One can also construct correlation networks from metabolomics data [74, 75].
Like gene co-expression, correlation between metabolites can occur for
multiple reasons, including knock-out of a gene coding an enzyme that is
involved in a chemical reaction consuming or producing two metabolites,
different temperatures or other environmental conditions under which different
samples are obtained, or intrinsic variability owing to cellular metabolism
[74]. Note that mass conservation within a moiety-conserved cycle produces
negative correlation between at least one pair of metabolites involved in the
reaction [76]. That said, in some cases one may consider correlation or other
similarity between only a subset of metabolites that are not necessarily
associated to one another by direct chemical processes but instead draw from a
set of alternative biochemical processes (see, e.g., [77]).

###  2.5 Microbiome networks

Microbes interact with other microbe species as well as with their
environments. Understanding of microbial composition and interaction in the
human gut is expected to inform multiple diseases. Similarly, understanding
soil microbial communities may contribute to enhancing plant productivity.
Network analysis is adept at revealing, e.g., ecological community assembly
and keystone taxa, and has been increasingly contributing to these fields.

In microbiome network analysis, one collects samples from, e.g., soil, at
various time points or locations. Each sample from an environment (e.g., soil,
gut, animal corpus, or water) contains various microorganisms with different
quantities. Co-occurrence network analysis is increasingly common in this
field, aided by an increasing amount and accuracy of data [26, 78, 79]. In a
microbiome co-occurrence network, nodes are microorganisms (e.g., bacteria,
archaea, or viruses), specified at the taxa level, for example, and an edge is
defined to exist if two nodes co-occur across the samples. Therefore,
microbiome co-occurrence networks are essentially microbiome correlation
networks, and the usual correlation measures, such as Pearson correlation, can
be used to determine edge data, but more sophisticated methods to define edges
are more often used. (See [26, 78] for various co-occurrence network
construction methods.) Positively weighted edges result because of, e.g.,
cooperation between two taxa, sharing of niche requirements, or co-
colonization. Negatively edges result because of, e.g., competition for space
or resources, prey-predator relationships, or niche partitioning. A
historically famous example of negative co-occurrence in ecological community
assembly study is the checkerboard-like presence-absence patterns of two bird
species inhabiting an island, discussed by Jared Diamond [80]. (Also see [26]
for a historical account.) Regardless, one should keep in mind that
correlation, or co-occurrence, does not immediately imply physical interaction
between two taxa.

###  2.6 Disease networks

A node in a disease network is a disease phenotype. Correlation between two
diseases defines an edge, and there are various definitions of edges as we
introduce in this section. Each definition of edge creates a different type of
disease network.

Comorbidity, also called multimorbidity [81], is the simultaneous occurrence
of multiple diseases within an individual. One cause of comorbidity is that
the same gene or disease-associated protein can trigger multiple diseases.
Other causes, such as environmental factors or behaviors, such as smoking, can
also result in comorbidity. A collection of potentially comorbid diseases can
be modeled as the nodes of a network, and the edges, which are based on
comorbidity or other similarity index between diseases [75], are correlational
in nature.

The authors of [82] constructed phenotypic disease networks (PDNs) in which
nodes are disease phenotypes. The edges are sample correlation coefficients or
a variant, and the samples are patients in a hospital claim record (i.e.,
Medicare claims in the US). Note that here one uses correlation for binary
variables because each sample (i.e., patient) is either affected or not
affected by any disease $i$. The authors found that, for example, patients
tend to develop illness along the edges of the PDN [82].

Similarly, prior work constructed a human disease network when two diseases
share at least one associated gene, which is similar in principle to the
phenotypic disease network despite that the edge of the human disease network
is not a conventional correlation coefficient [83] (also see [84]). Similarly,
an edge in a metabolic disease network is defined to exist when two diseases
are either associated with the same metabolic reaction or their metabolic
reactions are adjacent to each other in the sense that they share a compound
that is not too common [85]. (H2O and ATP, for example, are excluded because
they are too common.) Alternatively, in a human symptom disease network [86],
the edge between a pair of diseases is a correlation measure in which each
sample is a symptom. In other words, roughly speaking, the edge weight is
large when two diseases share many symptoms.

###  2.7 Financial correlation networks

Stocks of different companies are interrelated, and the prices of some of them
tend to change similarly over time. A common transformation of such financial
time series before constructing correlation matrices and networks is into the
time series of logarithmic return, i.e., the successive differences of the
logarithm of the price, given by

$x_{i}(t)=\ln\frac{z_{i}(t+1)}{z_{i}(t)},$ (1)

where $z_{i}(t)$ is the price of the $i$th financial asset at time $t$, such
as the closure price of the $i$th stock on day $t$, and
$t\in\\{0,\ldots,T-1\\}$. An advantage of this method is that $x_{i}(t)$ is
not susceptible to changes in the scale of $z_{i}(t)$ over time [87]. Then,
one constructs the correlation matrix for $N$ time series
$\\{x_{i}(1),\ldots,x_{i}(T-1)\\}$, where $i\in\\{1,\ldots,N\\}$.

Financial correlation matrices have been analyzed for decades. For example,
Markowitz’s portfolio theory provides an optimal investment strategy as vector
$\bm{w}=(w_{1},\ldots,w_{N})^{\top}$, where $w_{i}$ represents the fraction of
investment in the $i$th financial asset, and ⊤ represents the transposition
[5, 6]. The theory formulates the minimizer of the variance of the return,
$\bm{w}^{\top}C\bm{w}$, where $C$ is the covariance matrix, as the solution of
a quadratic optimization problem with the constraint that the expected return,
$\bm{w}^{\top}\bm{g}$, where $\bm{g}=(g_{1},\ldots,g_{N})^{\top}$, and $g_{i}$
is the expected return for the $i$th asset, is larger than a prescribed
threshold.

Financial correlation matrices have also been extensively studied in
econophysics research since the 1990s, with successful uses of random matrix
theory [87, 88, 89, 6, 90, 91, 92, 93] and network methods such as maximum
spanning trees [94] and more advanced methods (see [24] for a review). One
usually employs random matrix theory in this context to verify that most
eigenvalues of the empirical financial correlation matrices lie in the bulk
part of the distribution of eigenvalues for random matrices. Such results
imply that most eigenvalues of the empirical correlation matrices can be
regarded as noise, and one is primarily interested in other dominant
eigenvalues of the empirical correlation matrices whose values are not
explained by random matrix theory [88, 89, 90]. The largest eigenvalue is
usually not explained by random matrix theory and is often called the market
mode.

Other types of financial data are possible. For example, correlation networks
were constructed from pairwise correlation between the daily time series of
the investor’s behavior (i.e., the net volume of Nokia stock traded or its
normalized version) for two investors [95, 96]. One can also renormalize the
covariance matrix using other indices, such as momentum [97].

###  2.8 Bibliometric networks

Apart from microbiome studies, bibliometric and scientometric studies are
another research field in which co-occurrence networks are often used [98,
99]. For example, in an academic co-authorship network, a node represents an
author, and an edge represents co-occurrence (i.e., collaboration) of two
authors in at least one paper. One can weigh the edge according to the number
of coauthored papers or its normalized variant [100]. While keeping authors as
nodes, one can also create other types of co-occurrence networks, such as co-
citation networks in which an edge connects two authors whose papers are cited
together by a later paper, and keyword-based co-occurrence networks in which
an edge connects two authors sharing keywords associated with their papers.
Nodes of co-occurrence bibliometric networks can also be journals,
institutions, research areas, and so forth. These co-occurrence networks are
mathematically close to correlation networks and have been useful for
understanding research communities and specialities, communication among
researchers, interdisciplinarities, and the structure and evolution of
science, for example.

Various other web-based information has also been analyzed as co-occurrence
networks. For example, tags annotating posts in social bookmarking systems can
be used as nodes of co-occurrence networks [101]. Two tags are defined to form
an edge if both tags appear on the same post at least once. One can also use
the number of the posts in which the two tags co-appear as the edge weight.
Another example is co-purchase networks in online marketplaces, in which a
node represents an item, and an edge represents that customers frequently
purchase the two items together [102].

###  2.9 Climate networks

Climate can be analyzed as a network of interconnected dynamical systems [103,
104, 105, 106]. In most analyses, the nodes of the network are equal-angle
latitude-longitude grid points on the globe. However, such angular partitions
lead to grid cells with geometric areas that vary with latitude, which in
particular might lead to spurious correlations in the measured quantities,
especially near the poles; such biases might be addressed either by a node
splitting scheme that aims to obtain consistent weights for the network
parameters [107], or by choosing instead to work on a grid with (possibly only
approximately) equal grid cell areas [108]. Each node has, for example, a time
series measurement of the pressure level, which represents wind circulation of
the atmosphere. The edge between a pair of nodes is based on the correlation
between the two time series. An early study showed that all nodes in tropic
regions have large degree (i.e., the number of edges that a node has)
regardless of the longitude, whereas only a small fraction of nodes in the
mid-latitude regions had large degrees [103]. Climate networks have been
further used for understanding mechanisms of climate dynamics and predicting
extreme events. For example, early warning signals were constructed from the
degree of the nodes and clustering coefficient for climate networks of the
Atlantic temperature field [109]. The proposed early warning signals were
effective at anticipating the collapse of Atlantic Meridional Overturning
Circulation. See section 2.1 of [106] for more examples.

##  3 Methods for creating networks from correlation matrices

To apply network analysis to correlation matrix data, we need to generate a
network from correlation data (usually in the form of a correlation matrix).
We call such a network a _correlation network_. Whether correlation network
analysis works or is justified depends on this process. Although there are
various methods for constructing correlation networks from data, they have
pros and cons. Furthermore, there are various unjustified practices around
correlation network generation, which may yield serious limitations on the
power of correlation network analysis. In this section, we review several
major methods.

###  3.1 Estimation of covariance matrices

How to estimate covariance matrices from observed data when the matrix size is
not small is a long-standing question in statistics and surrounding research
fields. In particular, the sample covariance matrix, a most natural candidate,
is known to be an unreliable estimate of the covariance matrix. See [110, 111,
112, 6] for surveys on estimation of covariance matrices. Although the primary
focus of this paper is estimation of correlation networks, not correlation or
covariance matrices, we briefly survey a few techniques of covariance matrix
estimation in this section, including providing the notations and
preliminaries used in the remainder of this paper. This exposition is
important in particular because correlation network analysis in non-
statistical research fields such as network science and also various
applications often ignores statistical perspectives examined in the previous
studies.

We denote by $(x_{i\ell})$ an $N\times L$ data matrix, where $N$ is the number
of nodes to observe the signal from, and $L$ is the number of samples, which
is typically the length of the time series or the number of participants in an
experiment or questionnaire. The sample covariance matrix,
$C^{\text{sam}}=(C^{\text{sam}}_{ij})$, is given by

$C^{\text{sam}}_{ij}=\frac{1}{L-1}\sum_{\ell=1}^{L}(x_{i\ell}-\overline{x}_{i})(x_{j\ell}-\overline{x}_{j}),$
(2)

where

$\overline{x}_{i}=\sum_{\ell=1}^{L}\frac{x_{i\ell}}{L}$ (3)

is the sample mean of the signal from the $i$th node.222Note the $L-1$ in the
denominator of Eq. (2) is necessary to obtain an unbiased estimator. Because
Eq. (2) is essentially a spectral decomposition of $C^{\text{sam}}$, the rank
of $C^{\text{sam}}$ is at most $L$. One can understand this fact more easily
by rewriting Eq. (2) as

$C^{\text{sam}}=\frac{1}{L-1}\sum_{\ell=1}^{L}\bm{\tilde{x}}_{\ell}\bm{\tilde{x}}_{\ell}^{\top},$
(4)

where
$\bm{\tilde{x}}_{\ell}=(x_{1\ell}-\overline{x}_{1},\ldots,x_{N\ell}-\overline{x}_{N})^{\top}$
Therefore, $C^{\text{sam}}$ is singular if $L<N$, while the converse does not
hold true in general.

Let $X_{i}$ denote a random variable for $i\in 1,\ldots,N$. The true
covariance matrix $C=(C_{ij})$ is given by

$C_{ij}=E[(X_{i}-\mu_{i})(X_{j}-\mu_{j})],$ (5)

where $E$ represents the expectation, and $\mu_{i}=E[X_{i}]$ is the ensemble
mean of $X_{i}$. In statistics, $C^{\text{sam}}$ is regarded as an estimator
of $C$, called the Pearson estimator. Equation (5) implies that a covariance
matrix is a symmetric matrix. It is also a positive semidefinite matrix.
Conversely, a symmetric positive semidefinite-matrix is always a covariance
matrix.

For any finite $L$, the sample covariance matrix obeys the so-called Wishart
distribution with $L$ degrees of freedom, denoted by $W_{N}(L,C)$, under the
assumption that the $L$ samples are i.i.d. and obey the multivariate normal
distribution whose covariance matrix is $C$ [113, 114, 3, 6]. We obtain
$E[C^{\text{sam}}]=C$. In other words, the sample covariance matrix is an
unbiased estimator of the true covariance matrix. The variance of
$C^{\text{sam}}_{ij}$ is equal to $(C_{ij}^{2}+C_{ii}C_{jj})/L$. In fact,
$C^{\text{sam}}$ is a problematic substitute of $C$, and the use of
$C^{\text{sam}}$ in applications in place of $C$ tends to fail; see [6] for an
example in portfolio optimization. An intuitive reason why $C^{\text{sam}}$ is
problematic is that, if $L$ is not much larger than $N$, which is often the
case in practice, one would need to estimate many parameters from a relatively
few observations. Specifically, the covariance and correlation matrices have
$N(N+1)/2$ and $N(N-1)/2$ unknowns to infer, respectively, whereas there are
$L$ samples of vector $(x_{1\ell},\ldots,x_{N\ell})$ available [115]. If $N/L$
is not vanishingly small (called the large dimension limit or the Kolmogorov
regime [6]), then the estimation would fail. As an extreme example, if $L<N$,
then $C^{\text{sam}}$ is singular, but $C$ is in general nonsingular. Even if
$L\geq N$, matrix $C^{\text{sam}}$ may be ill-conditioned if $L$ is not
sufficiently greater than $N$, whereas $C$ may be well-conditioned.

Therefore, covariance selection to reduce the number of parameters to be
estimated is a recommended practice when $L$ is not large relative to $N$
[115]. One also says that we impose some structure on the estimator of the
covariance matrix, with the mere use of the sample covariance matrix as an
estimate of the true covariance matrix corresponding to no assumed structure.

A major method of covariance selection is to impose sparsity on the covariance
matrix or the precision matrix (also called the condensation matrix), which is
the inverse covariance matrix. Note that a sparse precision matrix does not
imply that the corresponding covariance matrix (i.e., inverse of the precision
matrix) is sparse. Graphical lasso (see section 3.6) is a popular method to
estimate a sparse precision matrix. Another major method to estimate a sparse
correlation matrix is to threshold on the value of the correlation to discard
node pairs with correlation values close to $0$ (see section 3.2). Another
common method of covariance selection, apart from estimating a sparse
covariance matrix, is covariance shrinkage (see [116] for a review). With
covariance shrinkage, the estimated covariance matrix is a linear weighted sum
of the sample covariance matrix, $C^{\text{sam}}$, and a much simpler matrix,
called the shrinkage target, such as the identity matrix [117, 118, 119] or
the so-called single-index model (which is a one-factor model in factor
analysis terminology and is an approximation of $C^{\text{sam}}$ by a rank-one
matrix plus residuals) [117]. Note that the shrinkage target is a biased
estimator of $C$. These and other covariance selection methods balance between
the estimation biases and variances.

We denote the correlation matrix by $\rho$. The sample Pearson correlation
matrix, denoted by $\rho^{\text{sam}}=(\rho^{\text{sam}}_{ij})$ is defined by

$\rho^{\text{sam}}_{ij}=\frac{\sum_{\ell=1}^{L}(x_{i\ell}-\overline{x}_{i})(x_{j\ell}-\overline{x}_{j})}{\sqrt{\sum_{\ell=1}^{L}(x_{i\ell}-\overline{x}_{i})^{2}}\sqrt{\sum_{\ell=1}^{L}(x_{j\ell}-\overline{x}_{j})^{2}}}=\frac{C^{\text{sam}}_{ij}}{\sqrt{C^{\text{sam}}_{ii}C^{\text{sam}}_{jj}}}.$
(6)

Note that $\rho^{\text{sam}}_{ii}=1$ $\forall i\in\\{1,\ldots,N\\}$. Also note
that every sample correlation matrix is a sample covariance matrix of some
data but not vice versa.

A correlation matrix is characterized by positive semidefiniteness, symmetry,
range of the entries only being $[-1,1]$, and the diagonal being equal to $1$
[120]. The set of full-rank correlation matrices for a fixed $N$ is called the
elliptope, which has its own geometric structure [121, 122].

For standardized samples
$y_{i\ell}=(x_{i\ell}-\overline{x}_{i})/\sqrt{C^{\text{sam}}_{ii}}$, the
Euclidean distance between vectors $(y_{i1},\ldots,y_{iL})$ and
$(y_{j1},\ldots,y_{jL})$ is given by

$d_{ij}^{2}=\sum_{\ell=1}^{L}\left(y_{i\ell}-y_{j\ell}\right)^{2}=2-2\rho_{ij}^{\text{sam}}.$
(7)

Therefore, given the sample correlation matrix,
$d_{ij}=\sqrt{2(1-\rho_{ij}^{\text{sam}})}$ defines a Euclidean distance [94,
87].

###  3.2 Dichotomization

In this and the following subsections, we present several methods to generate
undirected networks from correlation matrices.

A simple method to generate a network from the given correlation matrix is
thresholding, which in its simplest form entails setting a threshold $\theta$,
and placing an unweighted edge $(i,j)$ if and only if the Pearson correlation
$\rho_{ij}\geq\theta$; otherwise, we do not include an edge $(i,j)$. It is
also often the case that one thresholds on $\left|\rho_{ij}\right|$. There are
mainly two choices after the thresholding. First, we may discard the weight of
the surviving edges to force it to $1$, creating an unweighted network.
Second, we may keep the weight of the surviving edge to create a weighted
network. See Fig. 1 for these two cases. The literature often use the term
thresholding in one of the two meanings without clarification. In the
remainder of this paper, we call the first case _dichotomizing_ (which can be
also called binarizing), which is, precisely speaking, a shorthand for
“thresholding followed by dichotomizing”. We discuss dichotomized networks in
this section and threshold networks without dichotomization (yielding weighted
networks) in sections 3.3 and 3.6.

Dichotomizing has commonly been used across research areas. However,
researchers have repeatedly pointed out that dichotomizing is not recommended
for multiple reasons.

Figure 1: Thresholding a correlation matrix. We set the threshold at
$\theta=0.12$. If we only threshold the correlation matrix, we obtain a
weighted network. If we further dichotomize the threshold matrix, we obtain an
unweighted network. A different $\theta$ yields a different network in
general.

First, no consensus exists regarding the method for choosing the threshold
value [123, 27, 124, 125] despite that results of correlation network analysis
are often sensitive to the threshold value [126, 127, 124, 125, 128]. In a
related vein, if a single threshold value is applied to the correlation matrix
obtained from different participants in an experiment, which is typical in
neuroimaging data analysis and referred to as an absolute threshold [124,
129], the edge density can vary greatly across participants. Since edge
density is heavily correlated with many network measures, this can be seen as
introducing a confound into subsequent analyses and casts doubt on consequent
conclusions, e.g., that sick participants tend to have less small-world brain
networks than healthy controls. (In this example, a network with a large edge
density would in general yield a small average path length and large
clustering coefficient, leading to the small-world property, so that density
differences alone could have driven the observed effect.) An alternative
method for setting the threshold is the so-called proportional thresholding,
with which one keeps a fixed fraction of the strongest (i.e., most correlated)
edges to create a network, separately for each participant [124, 129]; also
see [130, 131, 132] for an early study. In this manner, the thresholded
networks for different participants have the same density of edges. However,
while the proportional thresholding may sound reasonable, it has its own
problems [129]. First, because different participants have different
magnitudes of overall correlation coefficient values, the proportional
threshold implies that one includes relatively weakly correlated node pairs as
edges for participants with an overall low correlation coefficients. This
procedure increases the probability of including relatively spurious node
pairs, which can be regarded as type I errors (i.e., false positives),
increasing noise in the resulting network. (Also see [133, 134] for discussion
on this matter.) Second, the overall correlation strength is often predictive
of, for example, a disease in question. The proportional threshold enforces
the same edge density for the different participants’ networks. Therefore, it
gives up the possibility of using the edge density, which is a simplest
network index, to account for the group difference. If one uses the absolute
threshold, the edge density is different among participants, and one can use
it to characterize participants. The edge density in the proportional
thresholding is also an arbitrary parameter.

Second, apart from false positives due to keeping small-correlation node pairs
as edges, correlation networks at least in its original form suffer from false
positives because pairwise correlation does not differentiate between direct
effects (i.e., nodes $i$ and $j$ are correlated because they directly
interact) and indirect effects (i.e., nodes $i$ and $j$ are correlated because
$i$ and $k$ interact and $j$ and $k$ interact). In other words, correlations
are transitive. The correlation coefficient is lower-bounded by [135]

$\rho_{ij}\geq\rho_{ik}\rho_{jk}-\sqrt{(1-\rho_{ik})^{2}(1-\rho_{jk})^{2}}.$
(8)

Equation (8) implies that if $\rho_{ik}$ and $\rho_{jk}$ are large, i.e.,
sufficiently close to $1$, then $\rho_{ij}$ is positive. Furthermore, this
lower bound of $\rho_{ij}$ is usually not tight, suggesting that $\rho_{ij}$
tends to be more positive than what Eq. (8) suggests when
$\rho_{ik},\rho_{jk}>0$ [136, 137]. This false positive problem is the main
motivation behind the definition of the partial correlation networks and
related methods with which to remove such a third-party effect, i.e.,
influence of node $k$ in Eq. (8). (See section 3.5.) Instead, one may want to
suppress false positives by carefully choosing a threshold value. Let us
consider the absolute thresholding. For example, if $i$ and $j$ do not
directly interact, $i$ and $k$ do, $j$ and $k$ also do, yielding
$\rho_{ij}=0.4$, $\rho_{ik}=0.7$ and $\rho_{jk}=0.6$, then setting
$\theta=0.5$ enables us to remove the indirect effect by $k$. However, it may
be the case that $i^{\prime}$ and $j^{\prime}$ do not directly interact,
$i^{\prime}$ and $k^{\prime}$ do, $j^{\prime}$ and $k^{\prime}$ also do,
yielding $\rho_{i^{\prime}j^{\prime}}=0.2$, $\rho_{i^{\prime}k^{\prime}}=0.3$,
and $\rho_{j^{\prime}k^{\prime}}=0.4$. Then, thresholding with $\theta=0.5$
dismisses direct as well as indirect interactions (that is, it introduces
false negatives). A related artifact introduced by the combination of
thresholding and indirect effects is that thresholding tends to inflate the
abundance of triangles, as measured by the clustering coefficient, and other
short cycles [74, 137]; even correlation networks generated by thresholding
randomly and independently generated data $\\{x_{i\ell}\\}$ have high
clustering coefficients [137]. This phenomenon resembles the fact that spatial
networks tend to have high clustering just because the network is spatially
embedded [138, 139].

Third, whereas thresholding has been suggested to be able to mitigate
uncertainty on weak links (including the case of the proportional thresholding
to some extent) and enhance interpretability of the graph-theoretical results
(e.g., [27]), thresholding in fact discards the information contained in the
values of the correlation coefficient. For example, in Fig. 1, thresholding
turns a correlation of $-0.1$ and $0.1$ into the absence of an edge.
Furthermore, if we dichotomize the edges that have survived thresholding, a
correlation of $0.2$ and $0.35$ are both turned into the presence of an edge.

There are various methods to try to mitigate some of these problems. In the
remainder of this section, we cover methods related to dichotomizing.

One family of solutions is to integrate network analysis results obtained with
different threshold values [27] (but see [140] for a critical discussion). For
example, one can calculate a network index, such as the node’s degree, denoted
by $\alpha$, as a function of the threshold value, $\theta$, and fit a
functional form to the obtained function $\alpha(\theta)$ to characterize the
node [127]. Similarly, one can calculate $\alpha$ for a range of $\theta$
values and take an average of $\alpha$ [141, 142]. In the case of group-to-
group comparison, an option along this line of idea is the functional data
analysis (FDA), with which one looks at $\alpha$ as a function of $\theta$
across a range of $\theta$ values and statistically test the difference
between the obtained function for different groups by a nonparametric
permutation test [143, 144]. In these methods, how to choose the range of
$\theta$ is a nontrivial question.

A different strategy is to determine the threshold value according to an
optimization criterion. For example, a method was proposed [125] for
determining the threshold value as a solution of the optimization of the
trade-off between the efficiency of the network [145] and the density of
edges. Another method to set $\theta$ is to use the highest possible threshold
that guarantees all or most (e.g., 99%) of nodes are connected [146].

The maximal spanning tree is an easy and classical method to automatically set
the threshold and guarantee that the network is connected [94]. One adds the
largest correlation node pairs as edges one by one under the condition that
the generated network is the tree. In the end, the maximal spanning tree
contains all the $N$ nodes, and the number of edges is $N-1$. The maximal
spanning tree also allows a hierarchical tree representation, which
facilitates interpretation [94, 147, 148]. However, the generated network is
extreme in the sense that it is a most sparse network among all the connected
networks on $N$ nodes, without any triangles. A variant of the maximum
spanning tree is to sequentially add edges with the largest correlation value
under the constraint that the generated network can be embedded on a surface
of a prescribed genus value (roughly speaking, the given number of holes)
without edge intersection [149]. If the genus is constrained to be zero, the
resulting network is a planar graph, called the planar maximally filtered
graph (PMFG). The PMFG contains $3N-6$ edges. The PMFG contains more
information than the maximum spanning tree, such as some cycles and their
statistics. Another related method is to use the $k$ nearest neighbor graph of
the correlation matrix, with which each $i$th node is connected at least to
the $k$ nodes with the highest correlation with the $i$th node [134]. Yet
another choice, which is designed for general weighted networks, is the
disparity filter, with which one preserves only statistically significant
edges to generate the _network backbone_ [150]. Note that, with these methods,
some lower-correlation node pairs are retained as edges and some higher-
correlation edges are discarded.

###  3.3 Weighted networks

A strategy for avoiding the arbitrariness in the choice of the threshold value
and loss of information in dichotomizing is to use weighted networks,
retaining the pairwise correlation value as the edge weight [60, 123].
Although there are numerous settings in network science where negative edge
weights are considered, they are generally more difficult to treat. (See
section 3.4.) As such, two common methods to create positively weighted
networks are (1) using the absolute value of the correlation coefficient as
the edge weight and (2) ignoring negatively weighted edges and only using the
positively weighted edges. Both methods dismiss some information contained in
the original correlation matrix, i.e., the sign of the correlation or the
magnitude of the negative pairwise correlation. Nonetheless, these
transformations are widely used because many methods are available for
analyzing general positively weighted networks, many of which are extensions
of the corresponding methods for unweighted networks. One can also use methods
that are specifically designed for weighted networks [151].

It should be noted that weighted networks share the problem of false positives
due to indirect interaction between nodes with the unweighted networks created
by dichotomization. We also note that, in contrast to thresholding (which may
be followed by dichotomization), node pairs with any small correlation (i.e.,
correlation coefficient close to $0$) are kept as edges in the case of the
weighted network. This may increase the uncertainty of the generated network
and hence of the subsequent network analysis results.

Thresholding operations in statistics literature to increase the sparsity of
the estimated covariance matrix often produce weighted networks. This is in
contrast to the dichotomization, which produces unweighted networks. Hard
thresholding in statistics literature refers to coercing
$C^{\text{sam}}_{ij}$, with $i\neq j$, to $0$ if
$\left|C^{\text{sam}}_{ij}\right|<\theta$ and keep the original
$C^{\text{sam}}_{ij}$ if $\left|C^{\text{sam}}_{ij}\right|\geq\theta$ [152,
110, 153, 154]. Soft thresholding [152, 155, 154] transforms
$C^{\text{sam}}_{ij}$ by a continuous non-decreasing function of
$C^{\text{sam}}_{ij}$, denoted by $f(C^{\text{sam}}_{ij})$, such that

$f(x)=\begin{cases}x-\theta&(x\geq\theta),\\\ 0&(-\theta<x<\theta),\\\
x+\theta&(x\leq-\theta).\end{cases}$ (9)

This assumption implies that, in contrast to hard thresholding, there is no
discontinuous jump in the transformed edge weight at
$C^{\text{sam}}_{ij}=\pm\theta$. Both hard and soft thresholding, as well as a
more generalized class of thresholding function $f(x)$ [154], do not imply
dichotomization and therefore generate weighted networks. In numerical
simulations, all these thresholding methods to generate weighted networks
outperformed the sample covariance matrix in estimating true sparse covariance
matrices [154]. The same study also found that there was no clear winner
between hard or soft thresholding, while combination of them tended to perform
somewhat better than other types of thresholding. Adaptive thresholding refers
to using threshold values that depend on $(i,j)$. Adaptive thresholding
performs better in terms of the approximation error and theoretically
converges faster than universal thresholding schemes (i.e., using a threshold
value independent of $(i,j)$) in numerical simulations [156].

###  3.4 Negative weights

Correlation matrices have negative entries in general. In the case of both
unweighted and weighted correlation networks, we often prohibit negative edges
either by coercing negative entries of the correlation matrix to zero or by
taking the absolute value of the pairwise correlation before transforming the
correlation matrix into a network. We prohibit negative edges for two main
reasons. First, in some research areas, it is often difficult to interpret
negative edges. In the case of multivariate financial time series, a negative
edge implies that the price of the two assets tend to move in the opposite
manner, which is not difficult to interpret. In contrast, when fMRI time
series from two brain regions are negatively correlated, it does not
necessarily imply that these regions are connected by inhibitory synapses, and
it is not straightforward to interpret negative edges in brain dynamics data
[157, 27]. Second, compared to weighted networks and directed networks, we do
not have many established tools for analyzing networks in which positive and
negative edges are mixed, i.e., signed networks. Signed network analysis is
still emerging [158]

In fact, negative edges may provide useful information. For example, they
benefit community detection because, while many positive edges should be
within a community, negative edges might ideally connect different communities
rather than lie within a community. Some community detection algorithms for
signed networks exploit this idea [158, 159]. Another strategy for analyzing
signed network data is to separately analyze the network composed of positive
edges and that composed of negative edges and then combine the information
obtained from the two analyses. For example, the modularity, an objective
function to be maximized for community detection, can be separately defined
for the positive network and the negative network originating from a single
signed network and then combined to define a composite modularity to be
maximized [160, 123]. While these methods are designed for general signed
networks, they have been applied to brain correlation networks [123, 159].

Another type of approach to signed weighted networks is nonparametric weighted
stochastic block models [161, 162], which are useful for modeling correlation
matrix data. Crucially, this method separately estimates the unweighted
network structure and the weight of each edge but in a unified Bayesian
framework. By imposing a maximum-entropy principle with a fixed mean and
variance on the edge weight, they assumed a normal distribution for the signed
edge weight. Because the edge weight in the case of correlation matrices,
i.e., the correlation coefficient, is confined between $-1$ and $1$, an ad-hoc
transformation to map $(-1,1)$ to $(-\infty,\infty)$ such as
$y=2\operatorname{arctanh}x=\ln\frac{1+x}{1-x}$ is applied before fitting the
model. One can assess the goodness of such an ad-hoc transformation by a
posteriori comparison with different forms of functions to transform $x$ to
$y$ using Bayesian model selection [162]. In this way, this stochastic block
model can handle negative correlation values. With this method, one can
determine community structure (i.e., blocks) including its number and
hierarchical structure.

###  3.5 Partial correlation

A natural method with which to avoid false positives due to indirect
interaction effects in the Pearson correlation matrix is to use the partial
correlation coefficient (as in, e.g., [163, 164, 165]). This entails measuring
the linear correlation between nodes $i$ and $j$ after partialing out the
effect of the other $N-2$ nodes. Specifically, to calculate the partial
correlation between nodes $i$ and $j$, we first compute the linear regression
of $X_{i}$ on $\\{X_{1},\ldots,X_{N}\\}\setminus\\{X_{i},X_{j}\\}$, which we
write as $X_{i}\approx\sum_{m=1;m\neq i,j}^{N}\beta_{i,m}X_{m}$, where
$\beta_{i,m}$ is the coefficient of linear regression. Similarly, we regress
$X_{j}$ on $\\{X_{1},\ldots,X_{N}\\}\setminus\\{X_{i},X_{j}\\}$, which we
write as $X_{j}\approx\sum_{m=1;m\neq i,j}^{N}\beta_{j,m}X_{m}$. The residuals
for $L$ samples are given by $\varepsilon_{i,\ell}=x_{i\ell}-\sum_{m=1;m\neq
i,j}^{N}\beta_{i,m}x_{m\ell}$ and
$\varepsilon_{j,\ell}=x_{j\ell}-\sum_{m=1;m\neq i,j}^{N}\beta_{i,m}x_{m\ell}$,
where $\ell\in\\{1,\ldots,L\\}$. The partial correlation coefficient, denoted
by $\overline{\rho}_{ij}^{\text{par}}$, is the Pearson correlation coefficient
between $\\{\varepsilon_{i,1},\ldots,\varepsilon_{i,L}\\}$ and
$\\{\varepsilon_{j,1},\ldots,\varepsilon_{j,L}\\}$.

In fact, the partial correlation coefficient between $i$ and $j$ ($j\neq i$)
is given by

$\overline{\rho}_{ij}^{\text{par}}=-\frac{\Omega_{ij}}{\sqrt{\Omega_{ii}\Omega_{jj}}},$
(10)

where $\Omega=C^{-1}$ is the precision matrix [113, 114]. Equation (10)
implies that $\Omega_{ij}=0$ is equivalent to the lack of partial correlation,
i.e., $\overline{\rho}_{ij}^{\text{par}}=0$. This conditional independence
property gives an interpretation of the precision matrix, $\Omega$.

Equation (10) also implies that the partial correlation can be calculated only
when $C$ is of full rank, whose necessary but not sufficient condition is
$L\geq N$. If $C$ is rank-deficient, a natural estimator of the $N\times N$
partial correlation matrix
$\overline{\rho}^{\text{par}}=(\overline{\rho}_{ij}^{\text{par}})$ is a
pseudoinverse of $C$. However, the standard Moore-Penrose pseudoinverse is
known to be a suboptimal estimator in terms of approximation error [166, 167],
while the pseudoinverse is useful for screening for hubs in partial
correlation networks [167]. If $C$ is of full rank, $\Omega$ as well as $C$ is
positive definite. Therefore, $\overline{\rho}^{\text{par}}$, whose diagonal
entries are forced to be $-1$, is negative definite. We can verify this by
rewriting Eq. (10) as $\overline{\rho}^{\text{par}}=-D^{-1/2}\Omega D^{-1/2}$,
where $D=\text{diag}(\Omega_{11},\ldots,\Omega_{NN})$ is the diagonal matrix
whose diagonal entries are $\Omega_{11}$, $\ldots$, $\Omega_{NN}$. If we
consider matrix $2I+\overline{\rho}^{\text{par}}$, where $I$ is the identity
matrix, as an alternative partial correlation matrix to force its diagonal
entries to $1$ instead of $-1$, the eigenvalues of the new partial correlation
matrix are upper-bounded by $2$ [168].

By thresholding the partial correlation matrix, or using an alternative
method, one obtains an unweighted or weighted partial correlation network,
depending on whether we further dichotomize the thresholded matrix. Because
the partial correlation avoids the indirect interaction affect, the network
created from random partial correlation matrices yields, for example, smaller
clustering coefficients [137] than if we had used the Pearson correlation
matrix.

While it apparently sounds reasonable to partial out the effect of the other
nodes to determine a pairwise correlation between two nodes, it is not
straightforward to determine when the partial correlation matrix is better
than the Pearson correlation one. First, Eq. (10) implies that extreme
eigenvalues of $\overline{\rho}^{\text{par}}$ are those of a normalized
precision matrix. Because the precision matrix is the inverse of the
covariance matrix $C$, extreme eigenvalues of
$\overline{\rho}^{\text{partial}}$ are derived from eigenvalues of $C$ with
small magnitudes. It is empirically known for, e.g., financial time series,
that small-magnitude eigenvalues of the covariance matrices are buried in
noise, i.e., not distinguishable from eigenvalues of random matrices [88, 89].
Therefore, the dominant eigenvalue of the precision matrix is strongly
affected by noise [169].

Second, the entries of $\overline{\rho}^{\text{par}}$ are more variable than
those of Pearson correlation matrices. Specifically, if $(x_{1},\ldots,x_{N})$
obeys a multivariate normal distribution, the Fisher-transformed partial
correlation of a sample partial correlation, i.e.,

$z_{ij}=\frac{1}{2}\ln\left(\frac{1+\overline{\rho}_{ij}^{\text{par,sam}}}{1-\overline{\rho}_{ij}^{\text{par,sam}}}\right),$
(11)

where $\overline{\rho}_{ij}^{\text{par,sam}}$ is the sample partial
correlation calculated through Eq. (10) with
$\Omega=\left(C^{\text{sam}}\right)^{-1}$, approximately obeys the normal
distribution with mean
$\frac{1}{2}\ln\left(\frac{1+\overline{\rho}_{ij}^{\text{par}}}{1-\overline{\rho}_{ij}^{\text{par}}}\right)$
and standard deviation $[L-3-(N-2)]^{-1/2}$. This result dates back to Fisher
(see e.g., [170]). In contrast, the corresponding result for the Fisher
transformation of the Pearson correlation coefficient is that the transformed
variable approximately obeys the normal distribution with mean
$\frac{1}{2}\log\left(\frac{1+\rho_{ij}^{\text{sam}}}{1-\rho_{ij}^{\text{sam}}}\right)$
and standard deviation $(L-3)^{-1/2}$ [170]. Therefore, the partial
correlation has more sampling variance than the Pearson correlation unless
$L\gg N$.

Third, partial correlation matrices typically have more negative entries and
smaller-magnitude entries than Pearson correlation matrices [171, 170].
Combined with the larger variation of the sample partial correlation than the
sample Pearson correlation discussed just above, the tendency that
$\overline{\rho}^{\text{par}}_{ij}$ has a smaller magnitude than $\rho_{ij}$
poses a challenge of statistically validating the estimated partial
correlation networks [170].

###  3.6 Graphical lasso and variants

Estimating a true correlation matrix, which contains $N(N-1)/2$ unknowns, is
an ill-founded problem unless the number of samples is sufficiently larger
than $N(N-1)/2$. A strategy to overcome this problem is to impose sparsity of
the estimated correlation network. A sparsity constraint enforces zeros on a
majority of matrix entries to suppress the number of unknowns to be estimated
relative to the number of samples. Imposing sparsity on estimated correlation
networks is a major form of covariance selection. _Structural learning_ refers
to estimation of an unknown network from data and usually assumes that the
given data obey a multivariate normal distribution and that the estimated
network is sparse. For a review with tutorials on this topic, see [38].

The Gaussian graphical model assumes that the precision matrix from the data
obeys a multivariate normal distribution and usually imposes sparsity of the
precision matrix [114]. In addition to reducing the number of unknowns to be
estimated, a motivation behind estimating a sparse precision matrix is that
$\Omega_{ij}=0$ is equivalent to the absence of conditional linear dependence
of the signals at the $i$th and $j$th nodes given all the other $N-2$
variables, which is easy to interpret. The graphical lasso is an algorithm for
learning the structure of a Gaussian graphical model [172, 173, 174, 175, 176,
177, 178]. The graphical lasso maximizes the likelihood of the multivariate
normal distribution under a lasso penalty (i.e., $\ell_{1}$ penalty), whose
simplest version is of the form
$\lambda\sum_{i,j=1}^{N}\left|\Omega_{ij}\right|$, where we recall that
$\Omega_{ij}$ is the $(i,j)$ entry of the precision matrix, and $\lambda$ is a
positive constant. This penalty term is added to the negative log likelihood
to be minimized. If $\lambda$ is large, it strongly penalizes positive
$\left|\Omega_{ij}\right|$, and the minimization of the objective function
yields many zeros of the estimated $\Omega_{ij}$. One can extend the lasso
penalty function in multiple ways, for example, by allowing $\lambda$ to
depend on $(i,j)$ and automatically determining $\lambda$ using an information
criterion. Other ways to regularize the number of nonzero elements in the
precision matrix than lasso penalty are also possible. (See e.g., [179, 180,
181].)

An alternative to the graphical lasso is to estimate sparse covariance
matrices rather than sparse precision matrices under lasso penalty [182, 183,
184, 185, 186]. With this approach, the consequence of imposing sparsity,
$C_{ij}=0$, corresponds to marginal independence between $X_{i}$ and $X_{j}$.

Most graphical lasso models and their variants do not model the estimation
problem relative to a null model correlation matrix. However, by estimating a
sparse correlation matrix that is different relative to a null model of
correlation matrix (see section 3.8 for null models), it was found that the
estimated correlation matrix gives a better description of the given financial
correlation matrix data than the graphical lasso and that the choice of the
null model also affects the performance [186]. By construction, this method
infers a set of edges that are not expected from the so-called correlation
matrix configuration model (see section 3.8 for details).

A test of significance of an edge, run on each edge, may sound like a natural
way to filter a network. However, this idea is not feasible because multiple
comparison with $N(N-1)/2$ estimates, each of whose significance would have to
be tested, is not practical given that $N(N-1)/2$ is usually large [38] and
also because the different edges are correlated with each other, particularly
if they share a node. Rather, the edges that survive the filtering method,
such as the graphical lasso, should be regarded to be sufficiently strong to
be included in the model [38].

###  3.7 Temporal correlation networks

Many empirical networks vary over time, including temporal correlation
networks [187], and many methods have been developed for analyzing time-
varying network data [188, 189, 190, 187]. If the given data is a multivariate
time series that is non-stationary, then correlation matrices computed from
the first 10% of the time points may be drastically different from that
computed from the last 10%. So, there is the possibility of greater
adaptability and better generalizability when one uses a time series of
correlation networks rather than just one. One can then apply various temporal
network analysis tools to the obtained temporal correlation networks.

A simple method to create dynamic correlation networks from multivariate time
series data is sliding-window correlation [191] (also called rolling-window
correlation in the finance literature; see e.g. [192]). With this method, one
considers time windows within the entire observation time horizon,
$t=\\{1,\ldots,t_{\max}\\}$. These time windows may be overlapping or non-
overlapping. Then, within each time window, one calculates the correlation
matrix and network. If there are $100$ time windows within $[1,t_{\max}]$,
then this method creates a temporal network of length $100$. Reliably
computing a single correlation matrix and a static correlation network from
multivariate time series requires a reasonable length (i.e., the number of
time points) of a time window. Generation of a reliable dynamic correlation
network requires longer data because one needs a multitude of such reasonably
long time windows. A limitation of sliding-window correlations is that they
are susceptible to large variability if the size of the time window is small,
whereas a large window size sacrifices sensitivity to temporal changes [193].

Early seminal reports analyzed temporal correlation networks of stock markets
by tracking financial indices and central nodes of the static correlation
network over more than a decade [148, 194, 195]. However, methods of dynamic
correlation networks have been particularly developed in brain network
analysis. In neuroimaging studies, in particular in fMRI studies, dynamic
correlation networks are known as dynamic (or time-varying) functional
connectivity networks [196, 197, 198, 191, 199]. Temporal changes in
functional (i.e., correlational) brain networks may represent neural
signaling, behavioral performance, or changes in the cognitive state, for
example. Patterns of time-varying functional networks may alter under a
disease. One can also analyze stability of community structure of temporal
correlation networks over time [200, 201, 202]. (See also [190] for temporal
community detection in financial correlation networks.) Many variations of the
method, such as on how to create sliding windows and how to cluster time
windows to define discrete-state transition dynamics, and change-point
detection are available in the field (e.g., see [199]). Many of these methods
should be applicable to multivariate time series data to generate temporal
correlation networks in other domains.

There are also methods for estimating dynamic precision matrices, proposed
outside neuroscience [203, 204, 205]. For example, the time-varying graphical
lasso (TVGL) formulates the inference of discrete-time time-varying precision
matrix as a convex optimization problem [205]. The objective function is
composed of maximization of the log likelihood under a lasso sparsity
constraint and the constraint that the precision matrix at adjacent time
points does not change too drastically, enforcing the temporal consistency.

###  3.8 Null models

A good practice to verify any structural or dynamical measurement $\alpha$ is
to compare the value of $\alpha$ in the given network with the value of
$\alpha$ achieved in random networks. This allows one to determine whether the
$\alpha$ value measured for the given network data is explainable purely by
the gross structural properties (e.g. edge density, degree distributions) of
the random graph family (when the $\alpha$ value is similar between the given
and random networks) or it is the result of other distinctive features of the
data (when the $\alpha$ value is statistically different between the given and
random networks). Many discoveries in network science owe to the fact that key
analyses have prudently implemented this practice by using or inventing
appropriate null models of networks. Already in one of the earliest seminal
papers on small-world networks, Watts and Strogatz compared the average path
length and the clustering coefficient of empirical networks with those of the
Erdős-Rényi random graph having the same number of nodes and edges [206].

One of the most popular null models for networks is the configuration model,
i.e., a uniform ensemble of networks under the constraint that the degree of
each node is conserved [207, 208], either exactly or in expectation. By
comparing given networks against a configuration model, one can reliably
quantify and discuss various network properties such as network motifs [209],
community structure [210], rich clubs [211], and core-periphery structure
[212]. The rationale behind the use of the configuration model is that the
node’s degree is heterogeneously distributed for many empirical networks and
that one generally wants to explore structural, dynamical, or other properties
of networks that are not immediate outcomes of the heterogeneous degree
distribution. One can extend the configuration model by imposing additional
constraints that the network under discussion is supposed to satisfy, such as
spatiality, i.e., the constraint that the nodes are embedded in a metric space
and the probability of an edge depends on the distance between the two nodes
[213]. See [214, 215, 207] for reviews of configuration models and best
practices for generating random realizations from such models.

We should similarly test properties found for given correlation networks
against appropriate null models. However, the usual configuration models are
not appropriate as null models of correlation networks because they are
significantly different from correlation networks derived from purely random
data [74, 137, 90, 190]. The expectation $\langle A_{ij}\rangle$ of the
$(i,j)$ entry of the adjacency matrix of the configuration model conditioned
on the degrees of all nodes, at least in the idealized and unrealistic regime
of weak heterogeneity of the degrees [214, 215], is equal to

$\langle A_{ij}\rangle=\frac{k_{i}k_{j}}{N\overline{k}},$ (12)

where $k_{i}=\sum_{j=1}^{N}A_{ij}$ is the degree of the $i$th node in the
original network, $A_{ij}$ is the entry of the _empirical_ adjacency matrix of
the original network (i.e., $A_{ij}=1$ if there is an edge between the $i$th
and $j$th nodes, and $A_{ij}=0$ otherwise), and $\overline{k}$ is the average
degree over all nodes. The above expected value also represents the
probability of independently connecting the $i$th and $j$th nodes in
realizations of the configuration model for networks. Note that, one can
indeed realize graphs in the configuration model by sampling edges
independently (at least if the constraint on the degree is ‘soft’, i.e.
realized only as an ensemble average over realizations [215]), correlation
matrices cannot be generated with independent entries, even in the null model
of independent signals. This is because, even under the null hypothesis of
independent realizations of the original time series, the correlation matrix
constructed from such time series still obeys the ‘metric’ (or triangular)
inequality in (8). We will elaborate more on this point below.

To see what Eq. (12) yields by merely replacing an empirical adjacency matrix
with an empirical correlation matrix $\rho_{ij}$, we proceed as follows [90].
We assume that each empirical signal is standardized in advance such that
$\text{Var}(X_{i})=1$, $\forall i\in\\{1,\ldots,N\\}$. In this way, we do not
need to distinguish between the correlation (i.e., $\rho_{ij}$) and covariance
(i.e., $C_{ij}$) matrices. We express the ‘degree’ as

$k_{i}=\sum_{j=1}^{N}\rho_{ij}=\sum_{j=1}^{N}C_{ij}=\sum_{j=1}^{N}\text{Cov}(X_{i},X_{j})=\text{Cov}(X_{i},X_{\text{tot}}),$
(13)

where Cov represents the covariance and $X_{\text{tot}}=\sum_{i=1}^{N}X_{i}$
is the ‘total’ signal. Then, we obtain

$N\overline{k}=\sum_{i=1}^{N}k_{i}=\text{Cov}(X_{\text{tot}},X_{\text{tot}})=\text{Var}(X_{\text{tot}}),$
(14)

where $\text{Var}(X_{\text{tot}})$ is the variance of $X_{\text{tot}}$. By
inserting the above quantities into (12) with $A_{ij}$ replaced by
$\rho_{ij}$, we obtain for the expected correlation matrix

$\langle\rho_{ij}\rangle=\frac{\text{Cov}(X_{i},X_{\text{tot}})\text{Cov}(X_{j},X_{\text{tot}})}{\text{Var}(X_{\text{tot}})}=\frac{\text{Cov}(X_{i},X_{\text{tot}})}{\sqrt{\text{Var}(X_{\text{tot}})}\sqrt{\text{Var}(X_{i})}}\cdot\frac{\text{Cov}(X_{j},X_{\text{tot}})}{\sqrt{\text{Var}(X_{\text{tot}})}\sqrt{\text{Var}(X_{j})}}=\rho(X_{i},X_{\text{tot}})\rho(X_{j},X_{\text{tot}}).$
(15)

Technically, this expected matrix is a covariance matrix because it is
symmetric, rank 1, and the only nonzero eigenvalue is positive [90, 216]. To
interpret the meaning of the above expression for $\langle\rho_{ij}\rangle$,
we recall the definition of the conditional three-way partial Pearson
correlation coefficient [113, 3]:

$\rho(X_{i},X_{j}\mid
X_{\text{tot}})=\frac{\rho(X_{i},X_{j})-\rho(X_{i},X_{\text{tot}})\rho(X_{j},X_{\text{tot}})}{\sqrt{1-\rho(X_{i},X_{\text{tot}})^{2}}\sqrt{1-\rho(X_{j},X_{\text{tot}})^{2}}}.$
(16)

We therefore conclude that the expected correlation matrix in (15) is a
correlation matrix of $N$ signals that satisfy the conditional independence
relationship

$\rho(X_{i},X_{j}\mid X_{\text{tot}})=0$ (17)

$\forall i,j(\neq i)\in\\{1,\ldots,N\\}$ [190].

However, when one generates a correlation network from the configuration
model, i.e. from a correlation matrix obeying (17), the generated network is
far from a typical correlation network generated by random data, due to the
triangular inequality mentioned above. To see an example of this, let us
revisit the example briefly explained in section 3.2. Let us consider purely
random data in which we generate each sample $x_{i\ell}$,
$i\in\\{1,\ldots,N\\}$, $\ell\in\\{1,\ldots,L\\}$ (where we recall that $L$ is
the number of samples for each node) as i.i.d. random numbers obeying a given
normal distribution. Then, we calculate the correlation matrix for
$\\{x_{i\ell}\\}$ and then a correlation network. For a broad class of
methods, including thresholding, the generated correlation network has a high
clustering coefficient [137], precisely because of the inequality (8).
Therefore, high clustering coefficients in correlation networks should not
come as a surprise. In contrast, networks generated by the ordinary
configuration model yield low clustering coefficients [217], disqualifying it
as a null model for correlation networks. If we use the usual configuration
model as the null model, we would incorrectly conclude that a given
correlation network has high clustering even if the network does not have
particularly high clustering among correlation networks. The configuration
model as null model also underperforms the simpler null model, called the
uniform null model (which is analogous to the random regular graph and to the
Hirschberger-Qi-Steuer (H-Q-S) algorithm explained below) on benchmark
problems of community detection when communities of different sizes are
present in single networks and those communities are detected by modularity
maximization [190].

Several papers have noted the need for null models specifically for
correlation networks [26, 137, 90, 218, 208]. We should use correlation
networks derived from a random correlation matrix as a null model. We stress
once again that a random correlation matrix is different from a random network
model (e.g., Erdős-Rényi model or configuration model), because of the
dependencies between entries. Similarly, many classes of random matrices are
not appropriate null models for correlation networks, either. For example, a
symmetric matrix whose all on-diagonal entries are $1$ and off-diagonal
entries are i.i.d. uniformly on $(-1,1)$ is almost never a correlation matrix
unless $N$ is small [219]. In this section, we introduce several null models
of correlation matrices. All the null models presented give distributions over
correlation matrices. Then, using any method introduced in previous sections
(e.g., by thresholding in various ways), one can define corresponding null
models over correlation networks.

A straightforward and traditional null model consists in shuffling the
original data, $\\{x_{i\ell}\\}$, based on which the correlation matrix is
calculated. This method is especially typical for multivariate time series
data. For example, one preserves the power spectrum of the time series at each
node while the time series is otherwise randomized [220, 221, 137, 208]. More
specifically, one Fourier transforms the time series at the $i$th node,
randomize the phase, and carry out the inverse Fourier transform. Then, for
the randomized multivariate time series, one calculates the correlation
matrix, which we use as control. One can impose additional constraints on the
randomization of time series depending on the properties other than the power
spectrum that one wants to have the null model preserve [222].

Another null model is based on the expected spectral properties of a
correlation matrix, rather than on the expected matrix itself. Note that the
expected correlation matrix, under the null hypothesis of independent signals,
is the identity matrix whose entries we denote as
$\rho^{\text{MG1}}_{ij}=\delta_{ij}$ (where $\delta_{ij}$ is the Kronecker
delta symbol) [90, 190]. This null model corresponds to an expectation under
white noise signals $\\{x_{i\ell}\\}$ that are independent for different nodes
$i\in\\{1,\ldots,N\\}$ and samples $\ell\in\\{1,\ldots,L\\}$. However the
sample pairwise correlation $\rho^{\text{sam}}_{ij}$ measured from the
signals, even under the null hypothesis of independence, will be different
from the identity matrix, unless $L\to\infty$. To take into account the
effects of finite $L$, it is convenient to look at the spectrum of the sample
correlation matrix. If we assume a standardized independent normal
distribution for $x_{i}$ $\forall i$, then $\rho^{\text{sam}}$ obeys the
Wishart distribution $W_{N}(I,L)$, where we recall that $I$ is the $N\times N$
identity matrix, the covariance matrix of a white noise signal. Note that the
expectation of $W_{N}(I,L)$ is $I$. Early studies of financial multivariate
time series data found that only a small number of the largest eigenvalues of
the covariance matrix deviate from the distribution of eigenvalues predicted
by random matrix theory for the Wishart ensemble [88, 89]. In other words,
only the leading eigenvalues and the associated eigenvectors outside the
prediction of random matrix theory, often called the market modes, are
significant relative to the random matrix null model. Random matrix theory has
provided useful null models for covariance/correlation matrix data, in
particular in financial time series data analysis; see [6] for a review. Also
see [120] for a review of random correlation as opposed to covariance
matrices.

Among many possible specific choices of null models for correlation matrices
based on random matrix theory, here we consider two models, which we denote by
$\rho^{\rm MG2}$ and $\rho^{\rm MG3}$, proposed in [90]. A given correlation
matrix is symmetric and positive semidefinite and therefore can be decomposed
as

$\rho^{\text{sam}}=\sum_{i=1}^{N}\lambda_{i}\bm{u}_{(i)}\bm{u}_{(i)}^{\top},$
(18)

where $\lambda_{i}(\geq 0)$ is the $i$th eigenvalue of $\rho^{\text{sam}}$,
and $\bm{u}_{(i)}$ is the associated normalized right eigenvector. The null
model correlation matrix $\rho^{\rm MG2}$ preserves the contribution of small
eigenvalues to Eq. (18), which are regarded to be noisy and described by
random matrix theory, and is given by

$\rho^{\rm
MG2}=\sum_{i=1;\lambda_{i}\leq\lambda_{+}}^{N}\lambda_{i}\bm{u}_{(i)}\bm{u}_{(i)}^{\top},$
(19)

where

$\lambda_{+}=\left(1+\sqrt{\frac{N}{L}}\right)^{2}.$ (20)

The boundary $\lambda_{+}$ originates from the Marchenko-Pastur distribution
(also known as the Sengupta-Mitra distribution) and represents the expected
largest eigenvalue under the null hypothesis of independent signals. Matrix
$\rho^{\rm MG2}$ is not a correlation matrix because its diagonal elements are
not equal to 1. However, this does not affect most network analysis because we
usually ignore diagonal elements or self-loops in correlation networks. Matrix
$\rho^{\rm MG2}$ represents a null model constructed only from the eigenvalues
of the empirical correlation matrix that are deemed to be noisy. Therefore
comparing the empirical matrix against $\rho^{\rm MG2}$ singles out properties
that cannot be traced back to noise [90]. By contrast, matrix $\rho^{\rm MG3}$
also preserves the contribution of the largest eigenvalue (if all the entries
of the associated eigenvector is positive) in addition to that of the noisy
eigenvalues and is given by

$\rho^{\rm
MG3}=\lambda_{\max}\bm{u}_{(\max)}\bm{u}_{(\max)}^{\top}+\sum_{i=1;\lambda_{i}\leq\lambda_{+}}^{N}\lambda_{i}\bm{u}_{(i)}\bm{u}_{(i)}^{\top},$
(21)

where $\max$ is the index for the dominant eigenvalue of $\rho^{\text{sam}}$.
All the entries of the dominant eigenvector, $\bm{u}_{(\max)}$, is positive if
there is a sufficiently strong common trend affecting all the $N$ signals. If
this common trend is so strong that all the entries of the correlation matrix
are positive, the Perron-Frobenius theorem ensures the positivity of the
dominant eigenvector. When this happens, the common trend obscures all the
mutual correlations among the signals. Matrix $\rho^{\rm MG3}$ deliberately
removes this global trend in addition to the noise, to reveal the underlying
structure. Therefore, properties of correlation matrices or correlation
networks that are not expected from $\rho^{\text{MG3}}$ represent those not
anticipated by the simultaneous presence of local noise and global trends. In
other words, they reflect the presence of mesoscopic correlated structures
such as correlation-induced communities [90]. As we will discuss in section
4.2, one can indeed use matrix $\rho^{\rm MG3}$ to successfully detect
communities of correlated signals.

As another null model, the H-Q-S algorithm, invented by Hirschberger, Qi, and
Steuer [223] is an equivalent of the Erdős-Rényi random graph in general
network analysis. Specifically, given the covariance matrix, $C^{\text{sam}}$,
the H-Q-S algorithm generates random covariance matrices, $C^{\rm HQS}$, under
the following constraints. First, the expectation of each on-diagonal entry of
$C^{\rm HQS}$ is equal to the average of the $N$ on-diagonal entries of
$C^{\text{sam}}$, denoted by $\mu_{\rm
on}\equiv(1/N)\times\sum_{i=1}^{N}C^{\text{sam}}_{ii}$. Second, the
expectation and variance of each off-diagonal entry of $C^{\rm HQS}$ are equal
to those of $C^{\text{sam}}$ calculated on the basis of all the off-diagonal
entries, denoted by $\mu_{\rm off}$ and $\sigma_{\rm off}^{2}$, respectively.
Optionally, one can also constrain the variance of the on-diagonal entries of
$C^{\rm HQS}$ [223] or use a fine-tuned heuristic variant of the algorithm
[137]. To implement the most basic H-Q-S algorithm without constraining the
variance of the on-diagonal entries of $C^{\rm HQS}$, we set

$L^{\rm HQS}\equiv\max\left(2,\lfloor\left(\mu_{\rm on}^{2}-\mu_{\rm
off}^{2}\right)/\sigma_{\rm off}^{2}\rfloor\right),$ (22)

where $\lfloor\cdot\rfloor$ is the largest integer that is not greater than
the argument. Then, we draw $N\times L^{\rm HQS}$ variables, denoted by
$x_{i\ell}$ (with $i\in\\{1,\ldots,N\\}$ and $\ell\in\\{1,\ldots,L^{\rm
HQS}\\}$), independently from the identical normal distribution with mean
$\sqrt{\mu_{\rm off}/L^{\rm HQS}}$ and variance $-\mu_{\rm off}/L^{\rm
HQS}+\sqrt{\mu_{\rm off}^{2}/(L^{\rm HQS})^{2}+\sigma_{\rm off}^{2}/L^{\rm
HQS}}$. Then, the H-Q-S algorithm sets the covariance matrix by

$C^{\rm HQS}_{ij}=\sum_{\ell=1}^{L^{\rm HQS}}x_{i\ell}x_{j\ell}\quad
i,j\in\\{1,\ldots,N\\}.$ (23)

It is known that $\langle C^{\rm HQS}\rangle_{ij}=\delta_{ij}\mu_{\rm
on}+(1-\delta_{ij})\mu_{\rm off}$. Therefore, the expectation of the
correlation matrix, $\rho^{\text{HQS}}$, is approximately given by
$\langle\rho^{\rm HQS}\rangle_{ij}=\delta_{ij}+(1-\delta_{ij})\mu_{\rm
off}/\mu_{\rm on}$. The degree of correlation networks is usually
heterogeneously distributed. This property is not captured by the H-Q-S
algorithm [224].

A configuration model for covariance matrices that preserves the expectation
of each row sum excluding the diagonal entry, which is equivalent to each
node’s degree in the case of the adjacency matrix of a conventional network,
can be defined as follows [218]. This algorithm, which we refer to as the
correlation matrix configuration model, preserves the expectation of each
diagonal entry of $C^{\text{sam}}$, or the variance of each variable, and the
expectation of each row sum excluding the diagonal entry, i.e.,
$\sum_{j=1;j\neq i}^{N}C^{\text{sam}}_{ij}$ $\forall i$, corresponding to the
degree of the $i$th node. Under these constraints, the correlation matrix
configuration model uses the distribution of $x_{i\ell}$, determined from the
maximum entropy principle. In fact, each
$(x_{1\ell},\ldots,x_{N\ell})^{\top}$, $\ell\in\\{1,\ldots,L\\}$,
independently obeys an identical multivariate normal distribution whose mean
is the zero matrix and covariance matrix is denoted by $C^{\text{cm}}$.
Therefore, the correlation matrix configuration model is the Wishart
distribution $W_{N}(C^{\text{cm}},L)$. Matrix $C^{\text{cm}}$ is of the form
$\left[(C^{\text{cm}})^{-1}\right]_{ij}=-\left[\delta_{ij}\cdot
2\alpha_{i}+(1-\delta_{ij})(\beta_{i}+\beta_{j})\right]$, where $\alpha_{i}$
and $\beta_{i}$ are parameters to be determined. One determines the values of
$\alpha_{i}$ and $\beta_{i}$ using a gradient descent algorithm [218] or by
reformulating the problem as a convex optimization problem and solving it
[186].

##  4 Network-analysis inspired analysis directly applied to correlation
matrices

As we already mentioned, a straightforward way to use null models for
correlation matrices as controls of correlation networks is to generate
correlation networks from the correlation matrix generated by the null model
or its distribution. In this section, we showcase another usage of null models
for correlation matrices, which is to conduct analysis inspired by network
analysis but directly on correlation matrix data with the help of null model
correlation matrices. Crucially, this scenario does not involve transformation
of a given correlation matrix into a correlation network. To explain the idea,
consider financial time series analysis using correlation matrices. Portfolio
optimization and random matrix theory directly applied to correlation matrix
data are among powerful techniques to analyze such data [6]. These methods do
not suffer from difficulties in transforming correlation matrices into
correlation networks because they do not carry out such a transformation. In
contrast, a motivation behind carrying out such a transformation is that one
can then use various network analysis methods. A strategy to take advantage of
both approaches is to adapt network analysis methods for conventional networks
to the case of correlation matrix data.

###  4.1 Degree

Many empirical networks show heterogeneous degree distributions such as a
power-law-like distribution [225, 217]; such networks are called scale-free
networks. The same holds true for the weighted degree of many networks [226].
Correlation networks are no exception, not much depending on how one
constructs a network from correlation matrix data [194, 227, 228, 229, 218].

If we do not transform the given correlation matrix into a network, the node’s
weighted degree represents how the node’s signal, $X_{i}$, is close to the
signal averaged over all the nodes, $X_{\text{total}}$, as shown in Eq. (13).
Previous research showed that the weighted degree calculated in this manner is
heterogeneously distributed for some empirical data, while the right tail of
the distribution is not as fat as typical degree distributions for
conventional empirical networks [218]. The results are qualitatively similar
when one calculates the weighted degree of the $i$th node as $\sum_{j=1;j\neq
i}^{N}|C_{ij}|$ or $\sum_{j=1;j\neq i;C_{ij}>0}^{N}C_{ij}$. Therefore,
heterogeneous degree distributions of the correlation network are not an
artifact of the thresholding or other operations for creating networks from
correlation data, at least to some extent.

###  4.2 Community detection

Community structure in networks is a partition of the nodes into (generally
non-overlapping) groups that are internally well connected and sparsely
connected across. One can detect communities with many different algorithms,
and one popular family of methods, despite some shortcomings [128], is
modularity maximization [230, 210], which aims at placing a higher-than-
expected number of edges connecting nodes within the same community.
Modularity with a resolution parameter $\gamma$ is defined by

$Q=\frac{1}{N\overline{k}}\sum_{i,j=1}^{N}\left(A_{ij}-\gamma\frac{k_{i}k_{j}}{N\overline{k}}\right)\delta_{g_{i},g_{j}},$
(24)

where we remind that $A_{ij}$ is the entry of the empirical adjacency matrix,
$g_{i}$ is the community in which the $i$th node is placed by the current
partition, and $\delta$ is again the Kronecker delta. Note the presence of the
term $k_{i}k_{j}/N\overline{k}$ coming from Eq. (12), which signifies the use
of the configuration model as standard null model in the ordinary modularity
for networks. Approximate maximization of $Q$ by varying $g_{1}$, $\ldots$,
$g_{N}$ for a given network identifies its community structure.

Community detection, in particular modularity maximization, is desirable for
correlation network data, too. Given that the original correlation matrix has
both positive and negative entries in general, a possible variant of
modularity maximization for correlation networks is to maximize a modularity
designed for signed networks. The modularity for signed networks may be
defined as a weighted difference of the modularity calculated for the positive
network (i.e., the weighted network only containing positively weighted edges)
and the modularity calculated for the negative network (i.e., the weighted
network only containing negatively weighted edges with the edge weight being
the absolute value of the correlation) [123]. However, this procedure assumes
that, in the null model, edges can be thought as independent (as in the model
described by Eq. (12)) and that positive edges can be randomized independently
of negative edges. We have seen that both assumptions are clearly not valid
for correlation matrices.

One can bypass the analogy with networks by directly computing and maximizing
a modularity function that is appropriately defined for correlation matrices
[90]. A viable redefinition of modularity for correlation matrices is given by

$Q^{\text{cor}}=\frac{1}{\mathcal{N}}\sum_{i,j=1}^{N}\left(\rho^{\text{sam}}_{ij}-\langle\rho_{ij}\rangle\right)\delta_{g_{i},g_{j}},$
(25)

where $\mathcal{N}=\sum_{i,j=1}^{N}\rho^{\text{sam}}_{ij}$ is a normalization
constant. Matrix $\langle\rho\rangle$ is a proper null model for the
correlation matrix, which one can scale by a resolution parameter $\gamma$ if
desired. Approximate maximization of $Q^{\text{cor}}$ provides communities in
correlation matrices. The crucial ingredient is the choice of the null model,
$\langle\rho_{ij}\rangle$. Depending on what features of the original data one
desires to preserve, the use of any of the models described in section 3.8 is
in principle legitimate, e.g., the white-noise (identity matrix) model
$\rho^{\text{MG1}}_{ij}=\delta_{ij}$, the noise-only model
$\rho^{\text{MG2}}_{ij}$, the noise+global model $\rho^{\text{MG3}}_{ij}$, or
the correlation matrix configuration model $C^{\text{cm}}_{ij}$. However, note
that in the summation in Eq. (25) some terms must be positive and some must be
negative in order to identify a nontrivial community structure, i.e., one
different from a single community enclosing all nodes. For example, if all the
entries of the empirical correlation matrix, $\rho^{\text{sam}}$, are
positive, then the null model $\rho^{\text{MG1}}$ will keep all terms in the
summation in Eq. (25) non-negative, and the resulting optimal partition will
be a single community [90]. By contrast, the use of $\rho^{\text{MG2}}$ yields
negative terms in the same situation unless a global trend is present. If a
global trend is present, one obtains the best results by using
$\rho^{\text{MG3}}$, which reveals group-specific correlations [90].
Modularity maximization using $\rho^{\text{MG3}}$ has successfully revealed
nontrivial community structure in time series of financial stock prices [90,
91, 92], credit default swaps [93], brain rhythms [231], and in single-cell
gene expression profiles [66].

###  4.3 Clustering coefficient

Clustering coefficients measure the abundance of triangles in a network [217].
A dominant definition of clustering coefficient for unweighted networks,
denoted by $\tilde{C}$, is given by [206]

$\tilde{C}=\frac{1}{N}\sum_{i=1}^{N}\tilde{C}_{i},$ (26)

where $\tilde{C}_{i}$ is the local clustering coefficient at the $i$th node
given by

$\tilde{C}_{i}=\frac{(\text{number of triangles involving the }i\text{th
node})}{k_{i}(k_{i}-1)/2}.$ (27)

The denominator of the right-hand side of Eq. (27) is a normalization constant
to ensure that $0\leq\tilde{C}_{i}\leq 1$.

One often measures clustering coefficients for both unweighted and weighted
correlation networks. There are various definitions of weighted clustering
coefficients [232, 233]. One definition [228] is given by $\tilde{C}^{\rm
wei,Z}=N^{-1}\sum_{i=1}^{N}\tilde{C}_{i}^{\rm wei,Z}$, where

$\tilde{C}_{i}^{\rm
wei,Z}=\frac{1}{\max_{i^{\prime}j^{\prime}}w_{i^{\prime}j^{\prime}}}\frac{\sum_{\begin{subarray}{c}1\leq
j,\ell\leq N\\\ j,\ell\neq
i\end{subarray}}w_{ij}w_{i\ell}w_{j\ell}}{\sum_{\begin{subarray}{c}1\leq
j,\ell\leq N\\\ j,\ell\neq i;j\neq\ell\end{subarray}}w_{ij}w_{i\ell}},$ (28)

and $w_{ij}$ $(=w_{ji}\geq 0)$ is the weight of edge $(i,j)$.

Many empirical networks show large unweighted or weighted clustering
coefficient values, and correlation networks are no exception. However, as we
pointed out in section 3.2, a high clustering coefficient of the correlation
network is at least partly due to pseudo correlation.

Given this background, clustering coefficients for correlation matrices were
proposed using a similar idea to the case of modularity directly defined for
correlation matrices [234]. Because correlation matrices are naturally
clustered if we threshold on the Pearson correlation matrix, the authors used
the three-way partial correlation coefficient or partial mutual information to
partial out the effect of a common neighbor of nodes $j$ and $\ell$ (say $i$)
to quantify partial connection between $j$ and $\ell$. In other words, we
measure the connectivity between neighbors of $i$ by, for example, the partial
correlation coefficient $\rho(X_{j},X_{\ell}\mid X_{i})$, which we abbreviate
as $\rho_{j\ell\mid i}$; the partial correlation coefficient is defined by Eq.
(16). Because there is no clear notion of neighborhood for correlation
matrices, we need to consider all triplets of different nodes, ($i$, $j$,
$\ell$). Then, as for the definition of the original clustering coefficient
for networks, they took the average of $\rho_{j\ell\mid i}$ over the $i$th
node’s neighbors $j$ and $\ell$ to define a local clustering coefficient for
$i$. For example, we define a local clustering coefficient for node $i$ as a
weighted average by

$C_{i}^{\rm cor,A}=\frac{\sum_{\begin{subarray}{c}1\leq j<\ell\leq N\\\
j,\ell\neq i\end{subarray}}\left|\rho_{ij}\rho_{i\ell}\rho_{j\ell\mid
i}\right|}{\sum_{\begin{subarray}{c}1\leq j<\ell\leq N\\\ j,\ell\neq
i\end{subarray}}\left|\rho_{ij}\rho_{i\ell}\right|}.$ (29)

Finally, as in the case of the clustering coefficient for networks, we define
the global clustering coefficient by $C^{\rm cor,A}=\sum_{i=1}^{N}C_{i}^{\rm
cor,A}/N$. This method borrows the idea of clustering coefficient from complex
network studies and tailors it for correlation matrix data. Clustering
coefficients $C_{i}^{\rm cor,A}$ and $C^{\rm cor,A}$ already partial out the
effect of pseudo correlation between $X_{j}$ and $X_{\ell}$ due to $X_{i}$.
However, we can still compare the observed clustering coefficient values
against those for null models to validate whether or not the clustering
coefficient values for the given data are significantly different from those
for the null model [218].

##  5 Software

In this section, we introduce freely available code useful for analyzing
correlation networks. Obviously, one can apply software for analyzing general
unweighted or weighted networks after thresholding (and optionally further
dichotomizing) the given correlation matrix data. There are numerous packages
for unweighted and weighted network analysis, which we do not mention here
without a few exceptions.

In gene correlation network analysis, Weighted Correlation Network Analysis
(WGCNA) is a popular freely available R software package [60]. WGCNA provides
various outputs, such as community structure, weighted clustering
coefficients, and visualization. WGCNA transforms the original edge weight,
denoted by $w_{ij}$ for edge $(i,j)$, by a so-called soft thresholding
transformation, i.e., by $\left|w_{ij}\right|^{\beta}$, where $\beta\geq 1$ is
a parameter333The soft thresholding here, coined in [60], is different from
the same term defined in section 3.3. It also does not belong to thresholding
in the sense used in this article (defined in section 3.2)., such that one
obtains an unsigned weighted network. Phiclust [66] is another community-
detection tool for correlations from single-cell gene expression data, derived
from random matrix theory (i.e., Wishart ensemble as described in section
3.8). It can be used to identify cell clusters with non-random substructure,
possibly leading to the discovery of previously overlooked phenotypes.
Phiclust is written in R and is freely available at Github [235] and Zenodo
under GNU General Public License V3.0 [236].

The Brain Connectivity Toolbox is a MATLAB toolbox for analyzing networks
[237]. It is also implemented in Python. Despite the name, many of the
functions provided by the Brain Connectivity Toolbox can be applied to general
network data, and many people outside neuroscience also use it. In relation to
correlation networks, this toolbox is particularly good at handling weighted
and signed networks, such as their degree, clustering coefficients, and
community structure.

The graph-tool module in Python provides powerful network analysis and
visualization functionality [238]. In relation to correlation networks, graph-
tool is particularly strong at network inference based on stochastic block
models.

The bootnet package in R can be used for estimating psychological networks
using graphical lasso, with a unique feature of being able to assess the
accuracy of the estimated network [38]. This package can be used for other
types of correlation matrix data as well. Also see [239, 38] for various R
packages for sparse and related estimation of the precision and covariance
matrices. For example, the qgraph package can also generate networks using the
graphical lasso and visualize Pearson and partial correlation matrices [240],
with beginner-friendly tutorials (e.g., [241]).

Graphical lasso is also implemented in Python, through the GraphicalLassoCV
function in the scikit-learn package [242, 243]. The sklearn.covariance
package also contains other functions for covariance estimation such as
covariance shrinkage. Table 2 of [18] lists other code packages for estimating
graphical models as well as different models.

Papers [90, 208] contain references to Python, MATLAB, and R codes for
generating null models of correlation matrices discussed in section 3.8. For
example, the “spatiotemporal modeling tools for Python” contains functions to
generate null model correlation matrices such as the H-Q-S model (named
Zalesky matching model in their package) and methods through generating
surrogate time series [222]. Another package in the list is the Scola, in
Python, which generates the H-Q-S model and the correlation matrix
configuration model [186]. Finally, a MATLAB package on MathWorks [244]
implements the null models based on random matrix theory, $\rho^{\rm MG2}$ and
$\rho^{\rm MG3}$, derived in [90] from the Wishart ensemble.

##  6 Outlook

We have reviewed various techniques to obtain and analyze networks generated
from correlation matrix data, which naturally arise across many domains.
Sections 2 and 3 emphasize for readers that there is not a single dominant
method. We also highlighted good practices and pitfalls of individual methods.
Naïve applications of network generation and analysis methods to correlation
matrix data can easily yield flawed results. We should be careful about both
how to generate correlation networks and how to analyze them. Specifically, in
resonance with previous reports, we explained that a simplistic dichotomizing,
which is widely used, is problematic for multiple reasons (see section 3.2).
The problems, i.e., lack of an established method to set the threshold value,
false positives, and losing the values of subthreshold correlation values, are
shared by weighted networks generated by thresholding (without dichotomizing).
Our recommendation is to resort to other methods, such as partial correlation
networks, graphical lasso (which is a partial correlation network method) and
its variants, and covariance shrinkage, which avoid some of these problems.
Another option is to still use the simplistic thresholding, possibly followed
by dichotomization, but combine it with a proper null model to avoid pitfalls.
Although thresholded correlation networks have high clustering no matter what
data we use, a proper null model (e.g., shuffling of $\\{x_{i\ell}\\}$) will
also produce thresholded correlation networks with high clustering. Therefore,
by comparing the results with those for the null model, one can avoid wrong
conclusions such as that almost all empirical correlation networks have high
clustering. We recall that null models for networks, most famously the
configuration model, are not a proper null model for correlation networks
because they generally do not originate from correlation matrices and do not
generally match key properties of typical correlation matrices. See section
3.8 for proper choices. We conclude this article by stating some promising
directions for future research.

Reproducibility of correlation networks arising from common approaches is a
major practical issue, as has been pointed out in the literature in psychology
and psychotherapy (e.g. [42, 48]) and microbiome analysis [79]. In these
research fields and others, it is often the case that only relatively few
samples are available given the size of matrix and network, $N$, we wish to
analyze. Especially if the number of samples, $L$, is smaller than or close to
$N$, the partial correlation matrix and spectra of random matrices carry large
variation, in part because they are inherently rank deficient. From the
statistical inference point of view, it is not sound to try to infer many
parameters, such as entries of the covariance matrix, from relatively few
observations. The recognition of such phenomena led to the idea of covariance
selection and various other methods. The amount of data needed for reliably
estimating correlation networks, i.e., power analysis in statistics, should be
further pursued for various correlation matrix data [38]. The development of
methods to help practitioners use correlation networks better (e.g., by
providing uncertainty quantification or clarifying the various noise trade-
offs) can be transformational. Despite these challenges, there is a pressing
need to understand complex systems of a very high dimension (i.e., $N\gg L$)
with correlational data. One approach to this problem is to formulate
estimation of large correlation networks as a computational rather than
statistical challenge, as a problem to be solved under runtime and memory
constraints, and to search feasible solutions in combination with machine
learning [1]. How to reconsile the statistical and computational types of
approach and deepen usage of machine learning and artificial intelligence
techniques to correlation network analysis may be a beneficial research
direction.

A key outstanding question is the treatment of low-correlation edges. On one
hand, we have surveyed attempts to remove “noise” edges (for example by
thresholding or graphical lasso), which is supposed to improve the overall
signal-to-noise ratio of the graph representation. Sparser models are more
parsimonious, easier to process quickly and with a lower memory footprint, and
amenable to a range of network science analysis tools. On the other hand, one
can argue that getting rid of low-correlation edges risks losing valuable
information (see section 3.2). In fact, it has been shown in the neuroscience
context that the low-correlation edges alone can have substantial predictive
power [143, 245, 246].

A related question is how to use _a priori_ domain knowledge to choose
appropriate preprocessing steps, such as what threshold to apply and whether
or not to dichotomize. For example, dichotomizing may be more appropriate when
the _a priori_ belief is that nodes are either coordinating or not, with no
appreciable difference in the degree of coordination when two nodes
coordinate. As another example, one could use random matrix theory on domain-
relevant distributions to compare the dominant eigenvectors or eigenvalues
before and after thresholding. This exercise may provide guidance on when
thresholding is unlikely to have adverse effects.

Another viable alternative to the current focus on trying to recover and
analyze the most accurate possible correlation matrix may be to treat the
constructed correlation network as inherently uncertain and to regard it as a
sample from a distribution of possible matrices as part of the analysis. For
example, when assessing community structure, it may make sense to focus on
structures that appear consistently across samples of correlation networks
drawn from the estimated distribution, even when exact correlation values
(perhaps especially the weaker correlations) considerably vary from sample to
sample. Although there are established ways to do this for general networks
[247], we are not aware of any work either applying these approaches to
correlation networks or developing theory and methods tailored to correlation
networks. Such approaches may leverage existing work on null models for
correlation networks, for example, as priors when forming a posterior
distribution to sample from. On the other hand, some studies have documented
the stability of the detected correlation-induced communities across time and
their robustness under change of temporal resolution [90, 93].

There are many multilayer network data sets, including multilayer correlation
matrix data sets, and various data-analysis methods for them [248, 249, 190,
250]. Examples include brain activity, where different layers of correlation
matrices correspond to, for example, different frequency bands [251, 252, 253,
254], or brain activity during different tasks [255]. In gene co-expression
networks, different layers correspond to, for example, different levels of co-
expression between gene pairs [256] or different tissue types such as
different organs [257]. Overlaying different methods to construct correlation
networks from one data set in each layer is another method to construct
multilayer correlation matrices [258]. There are methods for analyzing
multilayer correlation matrices and networks such as multilayer community
detection algorithms [190, 257]. However, methods that exploit the
mathematical nature of correlation matrices data are still scarce and left for
future work. Furthermore, multilayer networks are a major representation of
temporal network data, where each layer corresponds to one time window [259,
190, 187]. Therefore, methods for analyzing multilayer correlation network
data will also contribute to analysis of time-varying correlation network
data.

Similar to other studies, we emphasize the potentially negative effects of
thresholding, motivating our explanation of other methods for constructing
correlation networks. However, thresholding also has positive effects such as
reducing false positives by discarding edges with small weights including the
case of correlation networks. Such positive effects of thresholding may be
manifested in multilayer data. For example, aggregating layers in a multilayer
network and dichotomizing the aggregated edge weight can enhance detectability
of multilayer communities compared to no thresholding under certain conditions
[260, 261]. Furthermore, some layers in a multilayer network may be more
informative than other layers. While these arguments are for general
multilayer networks, many of them may directly apply to multilayer correlation
matrices.

In section 4, we showcased some methods for analyzing correlation matrix data
without transforming them into correlation networks but yet by benefiting from
ideas and formulations of network analysis. Although we only presented the
degree, modularity maximization and clustering coefficients, this approach is
generalizable to other structural properties and analysis methods formulated
for networks. Examples include various node centrality measures, motifs,
community detection methods apart from modularity maximization, rich clubs,
fractals, and simplicial complexes.

For a given $N$, the set of covariance matrices constitutes the positive
semidefinite cone, which is convex. Similarly, the set of full-rank
correlation matrices, which is a strict subset of full-rank covariance
matrices, is called the elliptope [121, 122]. Positive semidefinite cones and
elliptopes are manifolds and have their own geometric structure, which have
been suggested to be useful for measuring the similarity between pairs of
covariance or correlation matrices. Quantitative comparison of two covariance
and correlation matrices is useful for various tasks such as fingerprinting of
individuals, anomaly detection, and change-point detection in multivariate
time series data. A straightforward way to measure the distance between two
covariance/correlation matrices is to use a common matrix norm such as the
Frobenius norm (i.e.,
$\sqrt{\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\rho^{(1)}_{ij}-\rho^{(2)}_{ij}\right|^{2}}$
in the case of correlation matrices, where $\rho^{(1)}$ and $\rho^{(2)}$ are
two correlation matrices). However, research (in particular in neuroimaging
studies) has suggested that geodesic distance measures respecting the
structure of these manifolds is better at, for example, fingerprinting of
individuals from fMRI data [262, 263, 264]. In these geodesic distance
measures, one considers the tangent space at a given point $x$ on the
manifold, which corresponds to a correlation/covariance matrix. The so-called
exponential map provides a one-to-one mapping from the straight line segment
on the tangent space, which is essentially the Euclidean space, to the
geodesic from $x$ to $y$ on the manifold. The logarithm map is the inverse of
the exponential map. The geodesic distance between $x$ and $y$ is the length
of the geodesic and has a practical matrix algebraic formula. Multiple
reasonable definitions of such geodesic distances exist [263]. See [265, 266,
263] for mathematical expositions. Although these techniques are not for
correlation networks but for matrices, they may potentially benefit
understanding and algorithms for correlation networks. For example, can we
understand null models of correlation matrices as a projection onto a
submanifold of the entire elliptope? What are geometric meanings of
thresholding, dichotomizing, and other operations to create correlation
networks? Do we benefit by measuring distances between correlation networks
rather than between correlation matrices?

One complication that has not received enough attention is that many in-
practice comparisons involve ensembles of observed networks rather than single
networks. This is the case in most fMRI studies where networks are used. When
working with an ensemble of networks, one must make various decisions, such as
whether or not to ensure that edge density is constant across networks (see
section 3.2 for the absolute versus proportional threshold). The development
of mathematical theories for how to construct correlation-based networks for
ensembles may be helpful because most null models and other tools are only
oriented toward single networks. Multilayer approach and geometric approaches
to correlation networks and matrices, both of which cater to between-
network/matrix comparisons, are promising paths towards this goal.

A graphon is a symmetric measurable function $W:[0,1]^{2}\to[0,1]$. Given $W$,
we generate dense graphs in which there is an edge between the $i$th and $j$th
nodes with probability $W(u_{i},u_{j})$, where each $u_{i}$ $\forall i$
independently obeys the uniform density on $[0,1]$ [267]. In network science,
assigning a node weight, either from a probability distribution or empirical
data, and connecting two nodes probabilistically as a function of the two node
weights has been a major method to generate networks [268, 269, 270, 271, 272,
273, 274]. Basic correlation networks are equivalent to an extension of this
class of network models where $u_{i}$ is an $L$-dimensional vector of the
$i$th feature from $L$ samples, and $W$ is a criterion with which to determine
the edge. In fact, similar to the construction of correlation networks by
dichotomizing, dichotomizing functions have been used as $W$ to generate
networks with power-law degree distributions from scalar node weights that do
not have to obey long-tailed distributions [271, 273, 275, 276, 277].
Importing mathematical frameworks and methods from graphon-related models,
such as the limit of a sequence of dense networks, to correlation network
analysis may be an interesting idea. Those frameworks may be able to provide
null models for correlation networks or give more mathematical foundations of
correlation networks.

We have seen that there are various research fields in which they collect and
analyze correlation networks. We have also seen that some particular analysis
techniques are heavily studied in one field, but others are preferred in
different fields. For example, random matrices for correlation matrices have
particularly been used in financial data studies, including econophysics. As
another example, a majority of studies of temporal correlation networks has
been done in network neuroscience under the name of dynamic functional
connectivity/networks. However, very often, methods for analyzing correlation
networks developed in one research field do not rely on particularities of the
field and are therefore transferable to other research fields. While such
cross-fertilization has been ongoing and advocated [18], we emphasize that
much more of it will be useful for furthering correlation network analysis and
applications.

## Acknowledgements

N.M. acknowledges support from National Institute of General Medical Sciences
(under grant no. R01GM148973), the Japan Science and Technology Agency (JST)
Moonshot R&D (under grant no. JPMJMS2021), the National Science Foundation
(under grant nos. 2052720 and 2204936), and JSPS KAKENHI (under grant nos. JP
21H04595 and 23H03414). P.J.M. and Z.M.B. acknowledge support from the Army
Research Office (under MURI award W911NF-18-1-0244). P.J.M. also acknowledges
support from the National Institute of Diabetes and Digestive and Kidney
Diseases (under grant no. R01DK125860) and from the National Science
Foundation (under grant no. 2140024). Z.M.B. also acknowledges support from
the National Science Foundation (under grant no. 2137511). D.G. acknowledges
support by the European Union - NextGenerationEU - National Recovery and
Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR), project
‘SoBigData.it - Strengthening the Italian RI for Social Mining and Big Data
Analytics’ - Grant IR0000013 (n. 3264, 28/12/2021)
(https://pnrr.sobigdata.it/). The content is solely the responsibility of the
authors and does not necessarily represent the official views of any agency
supporting this work.

## References

  * [1] M. Becker, H. Nassar, C. Espinosa, I. A. Stelzer, D. Feyaerts, E. Berson, Neda H. Bidoki, A. L. Chang, G. Saarunya, A. Culos, D. De Francesco, R. Fallahzadeh, Q. Liu, Y. Kim, I. Marić, S. J. Mataraso, S. N. Payrovnaziri, T. Phongpreecha, N. G. Ravindra, N. Stanley, S. Shome, Y. Tan, M. Thuraiappah, M. Xenochristou, L. Xue, G. Shaw, D. Stevenson, M. S. Angst, B. Gaudilliere, and N. Aghaeepour.  Large-scale correlation network construction for unraveling the coordination of complex biological systems.  Nat. Comput. Sci., 3:346–359, 2023.
  * [2] I. T. Jolliffe.  Principal Component Analysis.  Springer, New York, NY, second edition, 2002.
  * [3] T. W. Anderson.  An Introduction to Multivariate Statistical Analysis.  John Wiley & Sons, Hoboken, NJ, third edition, 2003.
  * [4] S. A. Mulaik.  Foundations of Factor Analysis.  Taylor & Francis, Boca Raton, FL, second edition, 2010.
  * [5] H. Markowitz.  Portfolio selection.  J. Finance, 7:77–91, 1952.
  * [6] J. Bun, J. P. Bouchaud, and M. Potters.  Cleaning large correlation matrices: Tools from Random Matrix Theory.  Phys. Rep., 666:1–109, 2017.
  * [7] B. Podobnik and H. E. Stanley.  Detrended cross-correlation analysis: A new method for analyzing two nonstationary time series.  Phys. Rev. Lett., 100:084102, 2008.
  * [8] X.-Y. Qian, Y.-M. Liu, Z.-Q. Jiang, B. Podobnik, W.-X. Zhou, and H. E. Stanley.  Detrended partial cross-correlation analysis of two nonstationary time series influenced by common external forces.  Phys. Rev. E, 91:062816, 2015.
  * [9] N. Yuan, Z. Fu, H. Zhang, L. Piao, E. Xoplaki, and J. Luterbacher.  Detrended partial-cross-correlation analysis: A new method for analyzing correlations in complex system.  Sci. Rep., 5:8143, 2015.
  * [10] D. Y. Kenett, M. Tumminello, A. Madi, G. Gur-Gershgoren, R. N. Mantegna, and E. Ben-Jacob.  Dominating clasp of the financial sector revealed by partial correlation analysis of the stock market.  PLoS ONE, 5:e15032, 2010.
  * [11] D. Y. Kenett, T. Preis, G. Gur-Gershgoren, and E. Ben-Jacob.  Dependency network and node influence: Application to the study of financial markets.  Int. J. Bifu. Chaos, 22:1250181, 2012.
  * [12] D. Zhou, A. Gozolchiani, Y. Ashkenazy, and S. Havlin.  Teleconnection paths via climate network direct link detection.  Phys. Rev. Lett., 115:268501, 2015.
  * [13] L. Chen, R. Liu, Z.-P. Liu, M. Li, and K. Aihara.  Detecting early-warning signals for sudden deterioration of complex diseases by dynamical network biomarkers.  Sci. Rep., 2:342, 2012.
  * [14] S. Chen, E. B. O’Dea, J. M. Drake, and B. I. Epureanu.  Eigenvalues of the covariance matrix as early warning signals for critical transitions in ecological systems.  Sci. Rep., 9:2572, 2019.
  * [15] N. G. MacLaren, P. Kundu, and N. Masuda.  Early warnings for multi-stage transitions in dynamics on networks.  J. R. Soc. Interface, 20:20220743, 2023.
  * [16] T. Watanabe, N. Masuda, F. Megumi, R. Kanai, and G. Rees.  Energy landscape and dynamics of brain activity during human bistable perception.  Nat. Comm., 5:4765, 2014.
  * [17] T. Ezaki, T. Watanabe, M. Ohzeki, and N. Masuda.  Energy landscape analysis of neuroimaging data.  Phil. Trans. R. Soc. A, 375:20160287, 2017.
  * [18] D. Borsboom, M. K. Deserno, M. Rhemtulla, S. Epskamp, E. I. Fried, R. J. McNally, D. J. Robinaugh, M. Perugini, J. Dalege, G. Costantini, A.-M. Isvoranu, A. C. Wysocki, C. D. van Borkulo, R. van Bork, and L. J. Waldorp.  Network analysis of multivariate data in psychological science.  Nat. Rev. Methods Primers, 1:58, 2021.
  * [19] X. Liu, Y. Liu, J. Chen, and X. Liu.  PSCC-Net: Progressive spatio-channel correlation network for image manipulation detection and localization.  IEEE Trans. Circ. Syst. Video Technol., 32:7505–7517, 2022.
  * [20] M. Lee, C. Park, S. Cho, and S. Lee.  Superpixel group-correlation network for co-saliency detection.  In Proceedings of the International Conference on Image Processing (ICIP), pages 806–810, 2022.
  * [21] Z. Tu, Z. Li, C. Li, and J. Tang.  Weakly alignment-free RGBT salient object detection with deep correlation network.  IEEE Trans. Image Proc., 31:3752–3764, 2022.
  * [22] U. von Luxburg.  A tutorial on spectral clustering.  Stat. Comput., 17:395–416, 2007.
  * [23] M. Tumminello, C. Coronnello, F. Lillo, S. Micciche, and R. N. Mantegna.  Spanning trees and bootstrap reliability estimation in correlation-based networks.  Int. J. Bifu. Chaos, 17:2319–2329, 2007.
  * [24] M. Tumminello, F. Lillo, and R. N. Mantegna.  Correlation, hierarchies, and networks in financial markets.  J. Econ. Behav. Org., 75:40–58, 2010.
  * [25] S. M. Smith, K. L. Miller, G. Salimi-Khorshidi, M. Webster, C. F. Beckmann, T. E. Nichols, J. D. Ramsey, and M. W. Woolrich.  Network modelling methods for FMRI.  NeuroImage, 54:875–891, 2011.
  * [26] K. Faust and J. Raes.  Microbial interactions: From networks to models.  Nat. Rev. Microbiol., 10:538–550, 2012.
  * [27] F. De Vico Fallani, J. Richiardi, M. Chavez, and S. Achard.  Graph analysis of functional brain networks: Practical issues in translational neuroscience.  Phil. Trans. R. Soc. B, 369:20130521, 2014.
  * [28] Y. X. R. Wang and H. Huang.  Review on statistical methods for gene network reconstruction using expression data.  J. Theor. Biol., 362:53–61, 2014.
  * [29] A. Fukushima, M. Kusano, H. Redestig, M. Arita, and K. Saito.  Metabolomic correlation-network modules in Arabidopsis based on a graph-clustering approach.  BMC Syst. Biol., 5:1, 2011.
  * [30] T. Millington and M. Niranjan.  Construction of minimum spanning trees from financial returns using rank correlation.  Physica A, 566:125605, 2021.
  * [31] A. J. Butte and I. S. Kohane.  Mutual information relevance networks: Functional genomic clustering using pairwise entropy measurements.  In Pacific Symposium on Biocomputing, pages 418–429, Singapore, 2000. World Scientific.
  * [32] P. E. Meyer, F. Lafitte, and G. Bontempi.  minet: A r/bioconductor package for inferring large transcriptional networks using mutual information.  BMC Bioinformatics, 9:461, 2008.
  * [33] X. Guo, H. Zhang, and T. Tian.  Development of stock correlation networks using mutual information and financial big data.  PLoS ONE, 13:e0195941, 2018.
  * [34] R. Salvador, J. Suckling, C. Schwarzbauer, and E. Bullmore.  Undirected graphs of frequency-dependent functional connectivity in whole brain networks.  Phil. Trans. R. Soc. B, 360:937–946, 2005.
  * [35] P. Fiedor.  Partial mutual information analysis of financial networks.  Acta Physica Polonica A, 127:863–867, 2015.
  * [36] J. Pearl.  Causal inference.  Journal of Machine Learning Research Workshop and Conference Proceedings, 6:39–58, 2010.
  * [37] J. Runge, S. Bathiany, E. Bollt, G. Camps-Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M. D. Mahecha, J. Muñoz-Marí, E. H. van Nes, J. Peters, R. Quax, M. Reichstein, M. Scheffer, B. Schölkopf, P. Spirtes, G. Sugihara, J. Sun, K. Zhang, and J. Zscheischler.  Inferring causation from time series in Earth system sciences.  Nat. Commun., 10:2553, 2019.
  * [38] S. Epskamp, D. Borsboom, and E. I. Fried.  Estimating psychological networks and their accuracy: A tutorial paper.  Behav. Res., 50:195–212, 2018.
  * [39] G. M. Harari, N. D. Lane, R. Wang, B. S. Crosier, A. T. Campbell, and S. D. Gosling.  Using smartphones to collect behavioral data in psychological science: Opportunities, practical considerations, and challenges.  Perspectives on Psychological Science, 11:838–854, 2016.
  * [40] A. L. McGowan, F. Sayed, Z. M. Boyd, M. Jovanova, Y. Kang, M. E. Speer, D. Cosme, P. J. Mucha, K. N. Ochsner, D. S. Bassett, E. B. Falk, and D. M. Lydon-Staley.  Dense sampling approaches for psychiatry research: Combining scanners and smartphones.  Biol. Psychiatry, 93:681–689, 2023.
  * [41] D. Borsboom and A. O. J. Cramer.  Network analysis: An integrative approach to the structure of psychopathology.  Annu. Rev. Clin. Psychol., 9:91–121, 2013.
  * [42] E. I. Fried and A. O. J. Cramer.  Moving forward: Challenges and directions for psychopathological network theory and methodology.  Pers. Psychol. Sci., 12:999–1020, 2017.
  * [43] E. I. Fried, C. D. van Borkulo, A. O. J. Cramer, L. Boschloo, R. A. Schoevers, and D. Borsboom.  Mental disorders as networks of problems: A review of recent insights.  Soc. Psychiatry Psychiatr. Epidemiol., 52:1–10, 2017.
  * [44] A. L. McGowan, Z. M. Boyd, Y. Kang, L. Bennett, P. J. Mucha, K. N. Ochsner, D. S. Bassett, E. B. Falk, and D. M. Lydon-Staley.  Within-person temporal associations among self-reported physical activity, sleep, and well-being in college students.  Psychosomatic Medicine, 85:141–153, 2023.