Correlation networks: Interdisciplinary approaches beyond thresholding

There are various multivariate psychological data, from which one can construct networks [25, 21]. 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 [48, 49]. It is reasonable to use correlation between signals from different modalities (e.g., smartphones and brain scanners) to construct a correlation network [49].

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 [50, 51, 52, 21]. 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 [52]. Health-promoting behaviors can also be treated as nodes in these networks to suggest key behavioral intervention points [53].

Panel data, i.e., longitudinal measurements of variables from samples, are increasingly common for network approaches [21]. In this case, one obtains correlation networks at multiple time points. Then, one can construct time-varying correlation networks (see section 3.9) or within-person correlation networks [54] that reflect temporal symptom patterns and ideally expose individual differences and possible causal pathways in mental health patterns related to disorders [55, 56]. 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 [51, 57].

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 [58, 59]. (See also the related material about network representations in [60].) 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. A typical node in the 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. 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 a 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. Indeed, one does not expect a one-to-one correspondence between functional and structural brain networks because several brain states and functions continuously arise from the same anatomical structure [61]. Still, the possibility of studying structural networks in combination with functional networks on the same set of nodes is a distinguishing feature of brain networks, which can be used for an additional comparison when validating the outcome of correlation-based networks [62]. See [58, 32, 63, 64, 24, 62, 65] for reviews and comparisons of techniques for the estimation and validation of brain networks from the (partial) correlations. Examples of the use of these methods for (functional) network analysis are discussed later in this review.

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 [66]. 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 [67] 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 [68, 69]. 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) [70]. 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 [71]. 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.

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 [72, 73, 74, 75, 76, 34, 77, 78, 79]. 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 [80, 81]. 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 phenotype²²2A 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 [82].

Different ways of defining co-expression matrices and networks from gene expression data include tissue-to-tissue co-expression (TTC) networks [83] (also see [84, 85]). A TTC network proposed in [83] 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 [76]. 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 [76, 78]. 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 [86].

Metabolites are small molecules (e.g., amino acids, lipids, vitamins) that are intermediates or end products of metabolic reactions. One can also construct correlation networks from metabolomics data, or data of metabolites and their reactions [87, 88].To inform the edge, one measures pairwise correlation between the amounts of two metabolites given the samples. 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 [87]. Note that mass conservation within a moiety-conserved cycle produces negative correlation between at least one pair of metabolites involved in the reaction [89]. 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., [90]).

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 [33, 91, 92]. 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 [33, 91] 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 [93]. (Also see [33] for a historical account.) Regardless, one should keep in mind that correlation, or co-occurrence, does not immediately imply physical interaction between two taxa.

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 [94], 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 [88], are correlational in nature.

The authors of [95] 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 [95].

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 [96] (also see [97]). 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 [98]. (H₂O and ATP, for example, are excluded because they are too common.) Alternatively, in a human symptom disease network [99], 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.

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 [100]. 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 RMT [100, 101, 102, 6, 103, 104, 105, 106] and network methods such as maximum spanning trees [107, 108], community detection [109, 104, 110, 105, 106], and more advanced methods (see [23, 22] for reviews). One usually employs RMT 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 [101, 102, 103, 106]. The largest eigenvalue is usually not explained by RMT and is often called the market mode because it represents the movement of the market as a whole; moreover, other deviating eigenvalues are also present, encoding for the presence of groups of stocks that move coherently, as we discuss in section 4.3.

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 [111, 112]. One can also renormalize the covariance matrix using other indices, such as momentum [113].³³3Momentum in finance generally refers to the rate of change in price. If the prices of two assets are correlated, the momentum of one asset can be informative of the future price of the correlated asset [114]. In [113], momentum and price correlation are mixed in various ways to construct correlation-type networks that reflect collective price dynamics, and, for example, network centrality is predictive of large upcoming swings in asset prices.

Apart from microbiome studies, bibliometric and scientometric studies are another research field in which co-occurrence networks are often used [115, 116]. 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 [117]. 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 [118]. 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 [119].

Climate can be analyzed as a network of interconnected dynamical systems [120, 121, 122, 123]. 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 [124], or by choosing instead to work on a grid with (possibly only approximately) equal grid cell areas [125]. 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 equatorial 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 [120]. 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 [126]. The proposed early warning signals were effective at anticipating the collapse of Atlantic Meridional Overturning Circulation. See section 2.1 of [123] for more examples.

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 [127, 128, 129, 6] for surveys on estimation of covariance matrices. Although the primary focus of this paper is estimation of correlation networks, not covariance or correlation matrices, it is of course important to realize that correlation networks created from unreliably estimated correlation matrices are themselves unreliable. Therefore, we briefly survey a few techniques of covariance and correlation 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.⁴⁴4Note the $L-1$ in the denominator of Eq. (2) is necessary to obtain an unbiased estimator. Because Eq. (2) is a sum of $L$ outer products, 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.

Although covariance matrices are mathematically convenient, they are not normalized. In particular, if we multiply the data from the $i$th node by $c$ ($>0$), then $C_{ij}^{\text{sam}}$ $(=C_{ji}^{\text{sam}})$, where $j\neq i$, changes by a factor of $c$, and $C_{ii}^{\text{sam}}$ changes by a factor of $c^{2}$, whereas all the other entries of $C^{\text{sam}}$ remain the same. In practice, the data from different nodes may have different baseline fluctuation levels. For example, the price of the $i$th stock may fluctuate much more than that of the $j$th stock because the former has a larger average or the industry to which the $i$th company belongs may be subject to higher temporal variability. The correlation matrix, denoted by $\rho$, normalizes the covariance matrix such that $\rho$ is not subject to the effect of different overall amounts of fluctuations across different nodes. 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}}}.$$

(5)

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$ [130]. The set of full-rank correlation matrices for a fixed $N$ is called the elliptope, which has its own geometric structure [131, 132]. 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}}.$$

(6)

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

In the following, we refer to covariance matrices in some cases and correlation matrices in others, often interchangeably. This is because the correlation matrix, which is a normalized quantity, should be used in most data analyses, while the covariance matrix allows better mathematical analysis in most cases. This convention is not problematic because, if we standardize the given data first and then calculate the covariance matrix for the standardized data ($y_{i\ell}$), then the obtained covariance matrix is also a correlation matrix. Therefore, the mathematical results for covariance matrices also hold true for correlation matrices as long as we feed the pre-standardized data to the analysis pipeline.

With the above consideration in mind, we now stress that it is important to distinguish the sample covariance matrix, which is calculated from empirical data, from the theoretical or ‘true’ (also called population) covariance matrix. One may use the true covariance matrix to model the observed data mathematically in terms of a random process described by a (stationary) probability distribution. 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})],$$

(7)

where $E$ represents the expectation, and $\mu_{i}=E[X_{i}]$ is the ensemble mean of $X_{i}$. Equation (7) implies that a covariance matrix is a symmetric matrix. It is also a positive semidefinite matrix. Conversely, a symmetric positive semidefinite matrix $C$ is always a covariance matrix for the following reason. Any positive semidefinite matrix $M$ allows a positive semidefinite matrix, denoted by $M^{\frac{1}{2}}$, as square root. We set

$$\begin{pmatrix}X_{1}\\ \vdots\\ X_{N}\end{pmatrix}=C^{\frac{1}{2}}\begin{pmatrix}U_{1}\\ \vdots\\ U_{N}\end{pmatrix},$$

(8)

where $U_{1}$, $\ldots$, $U_{N}$ are independent random variables, with each having mean $0$ and variance $1$. Then, it is straightforward to see that $(X_{1},\ldots,X_{N})^{\top}$ has mean $0$ and covariance matrix $C$.

Having distinguished population and sample covariance matrices, we now look for a relationship between the two. If we regard $X_{i}$ as a random variable, and the observed data as a possible realization of that variable, then we must also regard the sample covariance matrix $C^{\text{sam}}$ as a random variable, and the empirical sample correlation matrix as a single realization of that variable. Then, the relevant question is the characterization of the probability distribution governing the sample covariance matrix, given the true covariance matrix $C$ which acts a set of fixed parameters (to be estimated from the data) for the distribution. Clearly, the number $L$ of observations, regarded as the number of independent draws for each random variable $X_{i}$, is another parameter (not to be estimated). 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$ [133, 134, 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, called the Pearson estimator in statistics. 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 [135]. 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 the true $C$ may be 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$ [135]. 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 so-called precision matrix (also called the concentration matrix), which is the inverse of the covariance matrix (with entries representing un-normalized partial correlation coefficients; see section 3.6). Note that a sparse precision matrix does not imply that the corresponding covariance matrix (i.e., its inverse) is sparse. Graphical lasso (see section 3.7) 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 [136] 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 [137, 138, 139] 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) [137]. Note that the shrinkage target is a biased estimator of $C$. These and other covariance selection methods balance between the estimation biases and variances.

An advanced estimation method for the entire correlation matrix, based on RMT, is the so-called optimal rotationally invariant estimator (RIE) [6, 65]. Roughly speaking, the RIE is the closest (in some spectral sense) matrix, among all those having the same eigenvectors as the sample correlation matrix, to the ‘true’ correlation matrix. It uses a certain self-averaging property to infer the spectrum of the true matrix from that of the sample matrix, and then to compute the spectral distance from the true matrix to be minimized [6]. A specific cross-validated version of the RIE has been recently shown to outperform several other estimators [65]. Since the RIE requires some notion of RMT to be properly defined, we discuss it in section 4.3.

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. 2 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.4 and 3.7.

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

First, no consensus exists regarding the method for choosing the threshold value [140, 141, 24, 142, 143, 144] despite that results of correlation network analysis are often sensitive to the threshold value [145, 141, 146, 142, 143, 144, 147]. 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 [142, 148], 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 [142, 148]; also see [149, 150, 151] 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 [148]. 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 [152, 153] 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 [154]

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

(9)

Equation (9) 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. (9) suggests when $\rho_{ik},\rho_{jk}>0$ [155, 156]. 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. (9). (See section 3.6.) 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 for dichotomized (and therefore unweighted) networks, and other short cycles [87, 156]; even correlation networks generated by dichotomizing randomly and independently generated data $\{x_{i\ell}\}$ have high clustering coefficients [156]. This phenomenon resembles the fact that spatial networks tend to have high clustering just because the network is spatially embedded [157, 158].

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., [24]), thresholding in fact discards the information contained in the values of the correlation coefficient. For example, in Fig. 2, 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 [24] (but see [159] 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 [146]. Similarly, one can calculate $\alpha$ for a range of $\theta$ values and take an average of $\alpha$ [26, 160]. 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 [161, 162]. 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 [143] for determining the threshold value as a solution of the optimization of the trade-off between the efficiency of the network [163] 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 [164].

The so-called maximal spanning tree is an easy and classical method to automatically set the threshold by guaranteeing that the network is connected [107, 100], while at the same time avoiding the creation of edges that would form loops (and are therefore unnecessary for connectedness). 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$. Thanks to the mapping from (large) sample correlation coefficients to (small) Euclidean distances established by Eq. (6), the maximal (in the sense of correlation) spanning tree is sometimes called the minimum (in the sense of distance) spanning tree [100, 107] (MST). The MST can be viewed as the graph achieving the overall minimum length among all graphs that make the $N$ nodes reachable from one another. Here, the minimum length is defined as the sum of the lengths of its realized edges, and the length of an edge is the metric distance between its endpoints. The maximal spanning tree also allows a hierarchical tree representation, which facilitates interpretation [107, 165, 166]. 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 [167]. 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 (which is in any case contained in the PMFG), such as some cycles and their statistics. Note that these methods effectively produce a threshold for the correlation value to be retained, but on the other hand preserve only some of the edges that exceed the threshold. Indeed, they are (designed to be) irreducible to the mere identification of an overall threshold value, with their merit residing in the introduction of higher-order geometric constraints guiding the dichotomization procedure.

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 [153]. 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 [168, 169]. Note that, with these methods as well, some lower-correlation node pairs are retained as edges and some higher-correlation edges are discarded.

Application example: Wang et al. compared functional networks of the brain between children with attention-deficit hyperactivity disorder (ADHD) and healthy controls using fMRI data [170]. The brain scan lasted for eight minutes in the resting state, and the brain image was acquired every two seconds. The authors used a previously published popular brain atlas defined by three-dimensional coordinates of 90 anatomical regions of interest (45 per hemisphere), each of which defines a node of the network. After various steps of preprocessing the original data, they computed the Pearson correlation coefficient between the time series at each pair of nodes, separately for each child. Aware of the problem inherent in choosing a threshold value, which we discussed earlier in this section, the authors examined dichotomized networks for a range of threshold values $\theta\in[0.05,0.5]$. As downstream analysis, they measured the small-world-ness of networks in terms of the so-called global efficiency and local efficiency [163]. A large global efficiency value and a small local efficiency value suggest that the network has a small-world property. They found that, while the networks of both ADHD and healthy control children were small-world, those of children with ADHD were somewhat less small-world than those of the controls across a wide range of $\theta$; the difference was statistically significant for the local efficiency in some range of $\theta$.

In the previous section, we discussed the idea of integrating analysis of dichotomized networks over different threshold values to mitigate the effect of selecting a single threshold value. Topological data analysis, or more specifically, persistent homology, provides a systematic and mathematically founded suite of methods to do such analyses. In topological analysis in general, one focuses on properties of shapes that do not change under continuous deformations of them, such as stretching and bending. Such topologically invariant properties include the numbers of connected components and of holes, which can be calculated through the so-called homology group. Persistent homology captures the changes in the network structure over multiple scales, or over a range of threshold values, clarifying topological features that are robust with respect to the threshold choice. In this broader perspective, networks are only a particular instance of the type of topological object under major consideration in topological data analysis, called simplicial complexes, because networks only consider pairwise interactions. One may want to consider the clique complex, in which each $k$-clique (i.e., complete subgraph with $k$ nodes) in the network is defined to be a higher-order object called a $(k-1)$-simplex and belongs to the simplicial complex of the given data. Note that the clique complex contains all the edges of the original network as well because an edge is a 2-clique by definition. For reviews of topological data analysis including persistent homology, see [171, 172, 173, 174, 175].

To analyze correlation matrix data using persistent homology, we start with a point cloud, with each point corresponding to a node in a correlation network. Then, we introduce a distance between each pair of nodes, $d_{ij}$, where $1\leq i,j\leq N$. By setting a threshold value $\theta^{\prime}$, we obtain a dichotomized network, or, e.g., a clique complex, depending on the choice, denoted by $G_{\theta^{\prime}}$. The simplicial complexes with varying $\theta^{\prime}$ values forms a filtration, i.e., with the nestedness property

$$G_{\theta^{\prime}_{1}}\subseteq G_{\theta^{\prime}_{2}}\subseteq G_{\theta^{\prime}_{3}}\subseteq\cdots,\mbox{ for }\theta^{\prime}_{1}\leq\theta^{\prime}_{2}\leq\theta^{\prime}_{3}\leq\cdots,$$

(10)

with the inclusion relationship in Eq. (10) referring to that of the edge set and that of any higher-order simplex object. The collection of clique complexes in Eq. (10) is called the Vietoris-Rips filtration. Simply put, with a larger threshold on the distance between nodes, the generated simplicial complex has more edges (and higher-order simplexes in the case of the clique complex). In practice, it suffices to consider the sequence of threshold values $\{\theta^{\prime}_{1},\theta^{\prime}_{2},\ldots\}$, such that $G_{\theta^{\prime}_{m+1}}$, contains just one more edge (and some higher-order simplexes containing that edge) than $G_{\theta^{\prime}_{m}}$; if there are multiple node pairs with exactly the same distance value, the corresponding multiple edges not in $G_{\theta^{\prime}_{m}}$ simultaneously appear in $G_{\theta^{\prime}_{m+1}}$. If $\theta^{\prime}$ is taken to vary over all possible topologies, the simplicial complex $G_{\theta^{\prime}_{1}}$ is composed of $N$ isolated nodes, while the last one in the nested sequence is the complete graph (that is, precisely, the corresponding clique complex) of $N$ nodes.

A large correlation should correspond to a small $d_{ij}$. One can realize this by setting $d_{ij}=f(\rho_{ij}^{\text{sam}})$, where $f$ is a monotonically decreasing function. Then, the network interpretation of Eq. (10) simply states the nested relationship in which the edges existing in a dichotomized network are included in dichotomized networks with smaller thresholds. However, while this practice is common, $d_{ij}$ is not mathematically guaranteed to be a distance metric, typically violating the triangle inequality if we use an arbitrary monotonically decreasing function $f$. Therefore, one often uses $f$ that makes $d_{ij}$ a distance metric, such as Eq. (6) or variants thereof. Then, the ensuing topological data analysis of correlation networks is underpinned by stronger mathematical foundations.

The next step is to calculate for each $G_{\theta^{\prime}}$ the homology groups or associated quantities, such as the $k$th Betti number. The zeroth and first Betti numbers are the numbers of connected components and essentially different cycles, respectively. These and other topological features of $G_{\theta^{\prime}}$ depend on $\theta^{\prime}$. For example, each connected component and cycle may appear (through closing loops) and disappear (through mergers with others) as one gradually increases the distance threshold $\theta^{\prime}$ from $\theta^{\prime}=0$, at which all the nodes are isolated. One can precisely visualize the occurrence of the birth and death events of each component by the persistence barcodes or persistence diagrams. For example, the persistence diagram represents each connected component, cycles, two-dimensional voids, etc. as a point $(x,y)$ in the two-dimensional space in which $x$ represents the $\theta^{\prime}$ value at which, e.g., a cycle appears and $y$ ($\geq x$) represents the $\theta^{\prime}$ value at which the same cycle disappears. If $y-x$ is large, that particular feature of the data is robust over changing scales because it is independent of the specific value of $\theta^{\prime}$ within a relatively large range.

We are often motivated to compare different correlation matrices or networks, such as in the comparison of functional brain networks between the patient group and control group. Quantitative comparison between two persistent diagrams provides a threshold-free method to effectively compare dichotomized networks. A persistence diagram consists of a set of two-dimensional points $(x,y)$. The bottleneck and Wasserstein distance, between persistence diagrams $P_{1}$ and $P_{2}$, which are commonly used, first consider the best matching between $(x,y)$ in $P_{1}$ and that in $P_{2}$. For the obtained best matching, both distance metrics measure the distance between $(x,y)$ in $P_{1}$ and that in $P_{2}$ in the Euclidean space and tally up the distance over all the points in different manners. For instance, the bottleneck distance is given by $d(P_{1},P_{2})=\inf_{\gamma}\sup_{(x,y)\in P_{1}}\left\|(x,y)-\gamma((x,y))\right\|_{\infty}$, where $\gamma$ is a matching between $P_{1}$ and $P_{2}$.

Persistent homology has been applied to correlation networks in neuronal activity data [141, 176, 177, 144], gene co-expression data [173, 178, 179], financial data [180, 181], and to co-occurrence networks of characters in literature work and in academic collaboration [182]. Scalability of persistent homology algorithms remains a concern. However, it may be of less concern for correlation network analysis because the number of nodes allowed for correlation networks is typically limited by the length of the data, not by the speed of algorithms (see section 3.1).

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 [73, 140]. Although there are numerous settings in network science where negative edge weights are considered, they are generally more difficult to treat. (See section 3.5.) 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 [183].

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$ [184, 127, 185, 186]. Soft thresholding [184, 187, 186] 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&\text{if }x\geq\theta,\\ 0&\text{if }-\theta<x<\theta,\\ x+\theta&\text{if }x\leq-\theta.\end{cases}$$

(11)

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)$ [186], 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 [186]. 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)$. An example is to use Eq. (11), but with an $(i,j)$-dependent threshold value, denoted by $\theta_{ij}$, in place of $\theta$, with $\theta_{ij}=\overline{c}\sqrt{\text{Var}_{ij}\cdot\log N/L}$, where $\overline{c}$ is a constant, and $\text{Var}_{ij}$ is an estimate of the variance of $(X_{i}-\mu_{i})(X_{j}-\mu_{j})$  [188]. This adaptive thresholding theoretically converges faster to the real covariance matrix and performs better numerically than universal thresholding schemes (i.e., using a threshold value independent of $(i,j)$).

Tapering estimators also threshold and suppress $C^{\text{sam}}_{ij}$ in an $(i,j)$-dependent manner. A tapering estimator for the $(i,j)$ entry of the covariance matrix is given by $\overline{f}(i,j)\cdot\frac{L-1}{L}C^{\text{sam}}_{ij}$, where $\overline{f}(i,j)=1$ if $|i-j|\leq k/2$, $\overline{f}(i,j)=2-2|i-j|/k$ if $k/2<|i-j|<k$, and $\overline{f}(i,j)=0$ if $|i-j|\geq k$. Here, tapering parameter $k$ is an even integer, and $\frac{L-1}{L}C^{\text{sam}}_{ij}$ is the maximum likelihood estimator of $C_{ij}$. The tapering estimator more strongly suppresses the $(i,j)$ entries that are farther from the diagonal respecting the sparsity of the estimated covariance matrix. This estimator is optimal in terms of the rate of convergence to the true covariance matrix, and the value of $k$ realizing the optimal rate differs between the two major matrix norms with which the estimation error is measured [189].

Application example: Chen et al. compared gene co-expression networks between people with schizophrenia and non-schizophrenic controls [190]. They used the Weighted Correlation Network Analysis (WGCNA) software (see section 6), which is frequently used in gene co-expression network analysis, to construct weighted networks of genes. In short, it uses $\left|\rho_{ij}\right|^{\beta}$ as the weight of edge $(i,j)$, where $\beta$ is a parameter. Then, they compared gene modules, or communities of the network, between the schizophrenia and control groups. The module detection was carried out by an algorithm implemented in WGCNA. There was no significant difference between the schizophrenia and control groups in terms of the module structure (i.e., which genes are in each module). However, the eigengene — that is, the first principal component of the co-expression matrix within a module — of each of two modules was significantly associated with schizophrenia compared to the control. The hub genes of each of these two modules, NOTCH2 and MT1X, were the largest contributor to the respective eigengenes. The authors further carried out biological analyses of these two genes to clarify where their expressions are upregulated and the functions in which these genes may be involved. The results were similar for a biopolar disorder data set.

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 [191, 24]. 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 [192]

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 [192, 193]. 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 [194, 140]. While these methods are designed for general signed networks, they have been applied to brain correlation networks [140, 193].

Another type of approach to signed weighted networks is nonparametric weighted stochastic block models [195, 196], 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 [196]. 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.

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., [197, 198, 199]). 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_{j,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}}},$$

(12)

where $\Omega=C^{-1}$ is the precision matrix [133, 134]. Equation (12) 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 (12) 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 [200, 201], while the pseudoinverse is useful for screening for hubs in partial correlation networks [201]. If $C$ is of full rank, $\Omega$ as well as $C$ is positive definite. Therefore, although Eq. (12) only holds true for $i\neq j$, if we denote the matrix defined by the right-hand side of Eq. (12) including the diagonal entries by $\tilde{\rho}$, then the diagonal entries of $\tilde{\rho}$ are equal to $-1$, and $\tilde{\rho}$ is negative definite. We can verify this by rewriting Eq. (12) as $\tilde{\rho}=-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 $\overline{\rho}^{\text{par}}=2I+\tilde{\rho}$, where $I$ is the identity matrix, as a partial correlation matrix to force its diagonal entries to $1$ instead of $-1$, the eigenvalues of $\overline{\rho}^{\text{par}}$ are upper-bounded by $2$ [202].

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 [156] 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. (12) 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 [101, 102], as we discuss later in section 4. Therefore, the dominant eigenvalue of the precision matrix is strongly affected by noise [203].

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),$$

(13)

where $\overline{\rho}_{ij}^{\text{par,sam}}$ is the sample partial correlation calculated through Eq. (12) 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., [204]). 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}$ [204]. 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 [205, 204]. 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 [204].

Studies in neuroscience have compared partial correlation networks with simple correlation networks and/or with the corresponding underlying structural networks [198, 206, 32, 207, 208, 209, 62] (see section 2.2). One study found that the similarity between partial correlation networks and structural networks is higher than that between correlation networks and structural networks [62].

Application example: Wang, Xie, and Stanley analyzed correlation networks composed of stock market indices from 2005 to 2014 from 57 countries [210], widely covering the continents of the world, with each country corresponding to a node. They computed Pearson and partial correlation coefficients for the time series of the logarithmic returns, given by Eq. (1), between each pair of countries, converted the correlation coefficient value into a distance (see Eq. (6)), and constructed the minimal spanning trees, called the MST-Pearson and MST-Partial networks. These networks appeared to be scale-free (i.e., with a power-law-like degree distribution) trees. Among other things, they compared clusters and top centrality nodes (i.e., countries) between the MST-Pearson and MST-Partial. They observed that the results from the MST-Partial networks are more reasonable than those from the MST-Pearson construction in light of our general understanding of world economics.

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 reviews with tutorials and examples on this topic, see [211, 25, 212].

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 [134]. 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 [213, 214, 215, 216, 217, 218]. 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}$. The $(i,j)$ pairs for which the estimated $\Omega_{ij}$ is nonzero form edges of the network; if there is no edge between $i$ and $j$, they are conditionally independent, i.e., conditioned on the other $N-2$ variables. This conditional independence is also referred to as the pairwise Markov property because the distribution of $X_{i}$ only depends on $\{X_{j}:(i,j)\text{ is an edge of the network}\}$. The pairwise Markov property is a special case of the global Markov property. The global Markov property dictates that $\{X_{i}:i\in A\}$ and $\{X_{i}:i\in B\}$ are conditionally independent given $\{X_{i}:i\in S\}$ if $S$ is a cutset of $A$ and $B$, that is, any path in the graph connecting a node in $A$ and a node in $B$ passes through $S$. The pairwise and global Markov properties are major cases of the Markov random field, which is defined by a set of random variables having Markov property specified by an undirected graph [219, 134, 220, 221, 222]. The network is sparse by design. 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., [223, 224, 225].)

A neighborhood selection method is another algorithm to estimate a sparse Gaussian graphical model [226]. With this method, one first carries out lasso regression for each $i$th variable (i.e., node) to determine a tentative small set of $i$’s neighbors. Second, if $j$ is a tentative neighbor of $i$, and $i$ is a tentative neighbor of $j$, then they are connected by an undirected edge $(i,j)$. The method named “space” (Sparse PArtial Correlation Estimation) advances the aforementioned neighborhood selection method by proposing to minimize a single composite penalized loss function, not separately minimizing the penalized loss function for each variable [227]. The loss function of “space” is a weighted sum of the regression error (i.e., error in estimating $X_{i}$ in terms of a linear combination of other $X_{j}$’s) over all $i$’s plus the usual $\ell_{1}$ penalty term. These regression-based methods as well the graphical lasso algorithms permit the case in which the number of observations, $L$, is smaller than the number of variables, $N$.

Bayesian variants of graphical lasso provide the level of undertainty in the estimated model. One such Bayesian approach assumes that we know node pairs that are not adjacent to each other, which is equivalent to imposing $\Omega_{ij}=0$ for a given set of node pairs $(i,j)$ [228]. Such a situation is possible when both the correlation matrix and structural network are available, as is common in MRI experiments in the brain. The Bayesian method uses the G-Wishart distribution, which is the Wishart distribution constrained by $\Omega_{ij}=0$ for non-adjacent node pairs in the given graph (i.e., structural network) $G$, as the prior distribution of the precision matrix, $\Omega$ [229, 230]. If $G$ is the complete graph, then the G-Wishart distribution is the Wishart distribution. Then, it uses data $(x_{i\ell})\in\mathbb{R}^{N\times L}$ to update the distribution using Bayes’ rule, obtaining the posterior distribution of $\Omega$, which is again a G-Wishart distribution but with updated parameter values. In other words, the G-Wishart distribution is a conjugate prior, giving us good intuitive pictures of how the prior distribution is transformed to the posterior distribution and mitigating the computational burden of the Bayesian method when $N$ is not small.

Gaussian graphical models assume normality. One approach to relax this assumption is to use a so-called nonparanormal distribution [231]. The idea is to assume that the transformed random variables $\left(f_{1}(X_{1}),\ldots,f_{N}(X_{N})\right)$ obey a multivariate normal distribution, where $f_{1}$, $\ldots$, $f_{N}$ are monotone and differentiable functions. In this case, we say that the original variables $(X_{1},\ldots,X_{N})$ have a nonparanormal distribution or a Gaussian copula (given by $f_{1}$, $\ldots$, $f_{N}$, and the mean vector and the covariance matrix of the normal distribution after the transformation). The nonparanormal distribution is nonidentifiable. For example, if we replace $f_{1}(X_{1})$ by a new function $\tilde{f}_{1}(X_{1})\equiv 2f_{1}(X_{1})$ and appropriately scale the first row and column of the covariance matrix and the first entry of mean vector used for the multivariate normal distribution that the transformed $N$ variables obey, then one gets the same distribution of $(X_{1},\ldots,X_{N})$. Therefore, we further impose that the transformations $f_{1}$, $\ldots$, $f_{N}$ conserve the mean and variance of the original $X_{1}$, $\ldots$, $X_{N}$. Then, as in the case of Gaussian graphical models, $\Omega_{ij}=0$ implies that $X_{i}$ and $X_{j}$ are conditionally independent of each other given all the other $N-2$ variables, in addition to that $f(X_{i})$ and $f(X_{j})$ are conditionally independent of each other. There are methods to estimate from data the functions $f_{1}$, $\ldots$, $f_{N}$ as well as a sparse covariance matrix $\Omega$, assuming a lasso penalty. Naturally, the nonparanormal distribution outperforms the graphical lasso when the true distribution of the data is not multivaraite normal [231]. Later work proposed algorithms for estimating the nonparanormal distribution with optimal convergence rate [232, 233]; the idea behind the improved algorithms is to avoid explicitly calculating $f_{1}$, $\ldots$, $f_{N}$ and to exploit statistics of rank correlation coefficient estimators. There are also other methods for estimating non-normal graphical models without relying on copula [234, 235].

Extreme cases of non-Gaussian distribution of the random variables are discrete multivariate distributions, in particular when each $X_{i}$ is binary (i.e., $\in\{-1,1\}$.). In this case, the form of the distribution of ($X_{1}$, $\ldots$, $X_{N}$) proportional to $e^{\bm{x}\Omega\bm{x}^{\top}}$, where $\bm{x}=(X_{1}-\mu_{1},\ldots,X_{N}-\mu_{N})$, defines the Ising model. Note that the same form of distribution for continuous variables is the multivariate normal distribution, which we discussed above for graphical lasso. There are many methods for inferring the Ising model from observed multivariate binary data, which is the task referred to as the inverse Ising model and Boltzmann machine learning [236, 237, 238, 239, 240]. Although exact likelihood maximization is available, it is practical only for small $N$. Therefore, various methods aim to overcome this limitation by allowing some approximation. Similar to graphical models, inference algorithms for the Ising model respecting the sparsity of the underlying network have also been investigated. For example, $\ell_{1}$-regularized methods based on pseudo-likelihood maximization were developed for estimating a sparse Ising Markov random field, or a graph in which the absence of the edge signifies the conditional independence between two nodes [241, 242, 243]. Decimation, or to recursively set the small edge weights to zero, combined with a stopping criterion and other heuristics, improves the performance of psuedo-likelihood maximization [244].

An alternative to the graphical lasso is to estimate sparse covariance matrices rather than sparse precision matrices under lasso penalty [245, 246, 247, 248, 249, 250]. With this approach, the consequence of imposing sparsity, $C_{ij}=0$, corresponds to marginal independence between $X_{i}$ and $X_{j}$. Similar to the case of the graphical lasso, one regards $(i,j)$ pairs for which the estimated $C_{ij}$ is nonzero as edges of the network.

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 4 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 [250]. By construction, this method infers a set of edges that are not expected from the so-called correlation matrix configuration model (see section 4 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 easily feasible because multiple comparisons 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 [25]. If we require a significance level of $0.05$, then Bonferroni correction for multiple comparisons implies that the edge statistic has to be tested to be significant with the $p$ value less than $0.05/[N(N-1)/2]$ for the edge to be actually significant at the $0.05$ significance level. Unless $N$ is small, this condition is usually too harsh. Furthermore, the different edges are correlated with each other, particularly if they share a node. In contrast, the Bonferroni correction assumes that the different tests are independent. Extensions of the Bonferroni corrections, such as the Šidák and Holm corrections, do not resolve these issues.

The false discovery rate approach, more specifically, the Benjamini-Yekutieli procedure [251], provides a solution to these problems. This test is less stringent than family-wise error rate controls including the Bonferroni correction and allows dependency between different tests. To run this procedure to generate a correlation network, we first calculate the $p$ value for each node pair $(i,j)$. Second, we arrange all the $p$ values in ascending order, which we denote by $p_{(1)}$, $p_{(2)}$, $\ldots$, $p_{(N(N-1)/2)}$. Then, we reject the null hypothesis (i.e., regard that the edge with the corresponding $p$ value is significant) for $p_{(1)}$, $\ldots$, $p_{(M)}$, where

$$M=\max\left\{i:p_{(i)}<\frac{0.05i}{\frac{N(N-1)}{2}\sum_{i^{\prime}=1}^{N(N-1)/2}\frac{1}{i^{\prime}}}\right\}.$$

(14)

This procedure is a thresholding method with an automatically set threshold value implied by the number of edges, $M$, determined by Eq. (14). We note that the type of correlation matrix (e.g., Pearson or partial correlation) does not matter once $p$ values have been obtained.

In the above procedure, what does it mean to calculate a $p$ value for a node pair? The answer varies. If we are given just one correlation matrix, which we assume in most of this article, the $p$ value can be that of the Pearson correlation coefficient in a standard textbook, if the Pearson correlation matrix is given as input [252] (see [253, 254, 255] for different calculations of the $p$ value when a single matrix is given as input). If we are given multiple correlation matrices forming one group, the one-sample $t$-test provides a $p$ value for each $(i,j)$ pair, quantifying whether the $(i,j)$ correlation for the group is different from $0$ [198]. If correlation matrices from two groups need to be compared, the two-sample $t$-test provides a $p$ value for each $(i,j)$. In this case, the generated network will be a difference network, in which the edges represent node pairs for which the correlation is significantly different between the two groups.

In neuroimaging research, the network-based statistic (NBS) is popularly used for controlling for the same multiple comparison problem [256]. In the NBS, one first calculates the $p$ value for each $(i,j)$, as in the Benjamini-Yekutieli procedure. Then, we keep $(i,j)$ whose $p$ values are smaller than an arbitrary threshold. If we just use those $(i,j)$ pairs to form a network, we would suffer from the aforementioned problems (i.e., multiple comparisons and dependencies between edges). Instead, the NBS focuses on the connected components induced by the surviving edges, which are small in number in general, and tests the significance of the sizes (i.e., the number of nodes) of the connected components using a permutation test. The NBS has been extended to a threshold-free method [257]. NBS and many of its extensions are applicable to correlation matrix data beyond neuroscience. Note that the NBS is not a method to estimate a correlation network; it tactically avoids network estimation and the problem of multiple comparisons, while providing a statistically controlled downstream network analysis (i.e., test on the size of connected components). In this sense, the NBS casts a key question: why do we need to estimate a network first of all? We will discuss this topic in section 7.1.

Covariance selection methods, such as the graphical lasso, do not explicitly test the significance of individual edges. The edges that survive these types of filtering methods should be regarded to be sufficiently strong to be included in the model [25].

Many empirical networks vary over time, including temporal correlation networks [258], and many methods have been developed for analyzing time-varying network data [259, 260, 110, 258]. 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 [261] (also called rolling-window correlation in the finance literature; see e.g. [262]). 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 [263].

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 [166, 264, 265]. 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 [266, 267, 268, 261, 269]. 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 [252, 270, 271]. (See also [110, 103, 104, 105, 106] for temporal stability analyses of 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 [269]). 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 [272, 273, 274]. For example, the time-varying graphical lasso (TVGL) formulates the inference of discrete-time time-varying precision matrix as a convex optimization problem [274]. 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.

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. In the simplest case, one randomizes all entries independently within each time series, thereby destroying all the cross-correlations while preserving the original values for each time series separately.

As a more constrained option, one preserves the power spectrum of the time series at each node while the time series is otherwise randomized [279, 280, 156, 277]. 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 is used as control.

Another method that preserves the full autocorrelation structure within each single time series, while randomizing cross-correlations among the $N$ time series, has been proposed in [281, 282]. The method is called the rotational random shuffling (RRS) model because it first imposes periodic boundary conditions (i.e., ‘gluing’ the last timestep to the first one) to turn each time series into a ‘ring’, and then randomly rotates the $N$ rings with respect to each other while keeping each ring internally intact.

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 [283].

Other null models are explicitly inspired by null models that are routinely used for networks. For instance, the H-Q-S algorithm, invented by Hirschberger, Qi, and Steuer [284] 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}$ [284] or use a fine-tuned heuristic variant of the algorithm [156]. 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),$$

(15)

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\}.$$

(16)

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}$. This highlights the completely homogeneous nature of the null model. For instance, while the degree of empirical correlation networks is usually heterogeneously distributed [108], this property is not captured by the H-Q-S algorithm [285].

One of the most popular heterogeneous 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  [286, 277], 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 [287], community structure [288], rich clubs [289], and core-periphery structure [290]. 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 [291]. See [292, 293, 286] 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 [87, 156, 103, 110]. 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 [292, 293], is equal to

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

(17)

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 [293]), 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 (9). We will elaborate more on this point below.

To see what Eq. (17) yields by merely replacing an empirical adjacency matrix with an empirical correlation matrix $\rho_{ij}$, we proceed as follows [103]. 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}}),$$

(18)

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}}),$$

(19)

where $\text{Var}(X_{\text{tot}})$ is the variance of $X_{\text{tot}}$. By inserting the above quantities into (17) 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}}).$$

(20)

Technically, this expected matrix is a covariance matrix because it is symmetric, rank 1, and the only nonzero eigenvalue is positive [103, 294]. 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 [133, 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}}}.$$

(21)

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

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

(22)

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

However, when one generates a correlation network from the configuration model, i.e. from a correlation matrix obeying (22), 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 distribution. Then, we calculate the sample correlation matrix for $\{x_{i\ell}\}$ and then a sample correlation network. This procedure immediately establishes a connection between null models for sample correlation matrices and RMT [7, 8, 9], some elements of which are discussed in section 4.3. For a broad class of methods, including dichotomizing, the generated correlation network has a high clustering coefficient  [156], precisely because of the inequality (9). 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 [295], disqualifying it as a null model for correlation networks [103]. 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 above) on benchmark problems of community detection when communities of different sizes are present in single networks and those communities are detected by modularity maximization [110].

A different configuration model, specifically designed for covariance matrices, can be defined as follows. This model 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 [276]. 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,L)$ that we have introduced in section 3.1, with $C=C^{\text{cm}}$. The 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 [276] or by reformulating the problem as a convex optimization problem and solving it [250].