Residual spectrum: Brain functional connectivity detection beyond coherence

Functional connectivity within the human brain is important for the understanding of brain disorders such as depression, Alzheimer’s disease, autism, dyslexia, dementia, epilepsy, and other types of mental disorder. The complete map of the human brain functional connectivity, or neural pathways, is referred to as the connectome. In describing the topology of the human brain, functional connectivity is commonly calculated by characterizing the co-variation of blood oxygen level-dependent (BOLD) signals, or time series, between two distinct neural populations or regions. To investigate the functional architecture of the brain, time series from resting-state functional magnetic resonance imaging (fMRI) which measures spontaneous low-frequency fluctuations in BOLD signals are often used.

Various time and spectral domain methods have been used in the analysis of BOLD signals. Among them is the wide spread measure of coherence (see, e.g., Müller et al., (2001) and Wang et al., (2015)). Coherence is a time series index used to determine if two or more brain regions are spectrally connected, that is, if two or more brain regions have “similar neuronal oscillatory activity” (Bowyer,, 2016). See Brockwell and Davis, (1991, Section 11, p.436) for the definition. However, coherence measures the strength of linear relationships only, but fails to measure nonlinear relationships, especially, the spectral contributions of quadratic interaction terms. The reasons why nonlinearities should be incorporated into the analysis of brain connections are explained in Zhao et al., (2019, p.1572) and Cao et al., (2022, Section 5.3).

The drawback is rectified by the so called maximum residual coherence, which provides the contribution of specific interaction terms to co-variation of brain regions, over and beyond the coherence by itself. Its origin can be traced back to the orthogonal decomposition of linear and quadratic functionals in the spectral domain proposed by Kimelfeld, (1974). Based on the decomposition, Kedem-Kimelfeld, (1975) developed a selection criterion, lagged coherence, to determine the inclusion of lagged processes or interaction terms in a linear system. This criterion is a function of frequency so that “we may run into a situation where one lag maximizes the coherence over a certain frequency band while another lag maximizes it over a different band” (Kedem-Kimelfeld,, 1975). The maximum residual coherence was developed to resolve such issue and shown to be effective in both cases of continuous and binary time series (Khan et al., 2014 and Kedem, 2016). The criterion is defined using supremum with respect to frequency, while Zhang and Kedem, (2021) proposed an alternative criterion using integration with respect to frequency to choose interaction terms in a more general setting.

Although the method based on maximum residual coherence is promising, no theoretical framework for its use in testing problems exists so far. Therefore, in this paper, we propose a new spectrum, referred to as the residual spectrum to capture non-linear association, and fill the gap by proposing, as a special case, a test for the existence of the interaction terms in covariation. We demonstrate that the disparity in the test results between between healthy individual and schizophrenia patients of our residual spectrum exceeds than that of the coherence. Additionally, we found that the non-linear connectivity for schizophrenia patients is higher than that for healthy individuals, whereas it is known that the reverse is true for linear connectivity.

In addition to coherence, Mohanty et al., (2020) listed alternative measures that can capture co-variation between brain signals in different ways including wavelet coherence, mutual information, dynamic time wraping, and more. Regarding other extensions of coherence, Ombao and Van Bellegem, (2008) proposed the local band coherence to deal with non-stationary signals. Fiecas and Ombao, (2011) studied brain functional connectivity by using the partial coherence estimated by the generalized shrinkage estimator, which is based on both nonparametric and parametric spectral density estimation. Lenart, (2011) dealt with the magnitude of coherence for almost periodically correlated time series and Baruník and Kley, (2019) proposed the quantile coherence based on quantile spectra. Euan et al., (2019) proposed a cluster coherence to measure similarity between vector time series and applied to electroencephalogram (EEG) signals. See also Matsuda, (2006) for a test to construct graphical models in the frequency domain and Liu et al., (2021) for the concept of local Granger causality and its associated testing procedures.

However, our goal is to highlight a way to extend and enhance the measure of coherence in the analysis of BOLD signals by incorporating interaction terms. Cao et al., (2022) gave a comprehensive review of recent developments on statistical methods to analyze brain functional connectivity.

For hypothesis testing problems described by spectra, $L_{2}$- and supremum-type statistics have been studied by many authors. As for $L_{2}$-type statistics, Taniguchi et al., (1996) discussed the test based on the integrated function of spectral density matrices (see also Taniguchi and Kakizawa, 2000, Chapter 6), which can be applied to, e.g., the test for the magnitude of linear dependence and discriminant analysis (Kakizawa et al., (1998)). A closely related problem was examined in Yajima and Matsuda, (2009) under Gaussian assumption. To include other hypotheses, Eichler, (2008) considered the testing problem described by the integrated squared (Euclidean) norm of the function of spectral density matrices. In the same spirit of Eichler, (2008), tests for the equality of spectra are proposed using the bootstrap method (Dette and Paparoditis,, 2009) and the randomization method (Jentsch and Pauly,, 2015). On the other hand, supremum-type statistics are considered by Woodroofe and Van Ness, (1967), Hannan, (1970, Thorem 12, Chapter V.5), Rudzkis, (1993), and Wu and Zaffaroni, (2018). However, it is difficult to obtain asymptotics of supremum statistics for the function of a spectral density matrix in general, but it can for $L_{2}$-statistics. Therefore, in this article, we propose an $L_{2}$-type statistic.

The paper is organized as follows. Section 2 introduce the residual spectrum, which is the main objective in this paper. Section 3 delves into time series multiple regression models and assumptions, and show there exists a unique orthogonal representation of the model to derive the residual spectrum and provide an interpretation of it. In Section 4, we propose a test for the existence of the residual spectrum. The asymptotic null distribution of our test statistic and consistency of the test are derived. The consistency of estimators of unknown parameters in the model is also addressed. Section 5 provides the finite sample performance of our method. Section 6 shows the utility of our method by applying the proposed test to brain data. All proofs of Lemma and Theorems in the main article are presented in Appendix 6.

In this section, we introduce a new spectrum, termed the residual spectrum of order $j$, for the vector stationary process $(X_{0}(t),X_{1}(t),\ldots,X_{K}(t))$, designed to capture both linear and non-linear terms of $(X_{1}(t),\ldots,X_{K}(t)$ in explaining $X_{0}(t)$. For the sake of brevity in notation, denote the spectral density matrix for $(X_{0}(t),X_{1}(t),\ldots,X_{K}(t))$ and $(X_{1}(t),\ldots,X_{K}(t))$ by $\bm{f}(\lambda):=\left(f_{ij}(\lambda)\right)_{i,j=0,\ldots,K}$ and $\bm{f}_{K}(\lambda):=\left(f_{ij}(\lambda)\right)_{i,j=1,\ldots,K}$, respectively.

Under the regularity conditions, as derived in the subsequent section, the residual spectrum $f_{G_{j}G_{j}}$ of order $j$ is defined as follows: for $j=1$,

where, for $j\in\{2,\ldots,K\}$,

and

The interpretation of our spectrum, as elucidated in the subsequent section, is as follows:

1. 1.

   $f_{G_{1}G_{1}}(\lambda)$ quantifies the strength of the linear association between $X_{0}(t)$ and $X_{1}(t)$.

2. 2.

   $f_{G_{2}G_{2}}(\lambda)$ measures the strength of the linear relationship between $X_{0}(t)$ and $X_{2}(t)$ after eliminating the linear effect of $X_{1}(t)$.

3. 3.

   $f_{G_{3}G_{3}}(\lambda)$ indicates the strength of the linear association between $X_{0}(t)$ and $X_{3}(t)$ after removing the linear effects of $X_{1}(t)$ and $X_{2}(t)$.

4. 4.

   Similarly, $f_{G_{j}G_{j}}(\lambda)$ represents the strength of the linear relationship between $X_{0}(t)$ and $X_{j}(t)$ subsequent to the removal of the linear effects of $X_{1}(t),\ldots,X_{j-1}(t)$.

While one might initially presume that $f_{G_{j}G_{j}}(\lambda)$ is incapable of capturing non-linear relationships, it does indeed possess this capability. We shall elaborate on the integration of non-linear relationships through our residual spectrum. Let $X_{0}(t)$ and $X_{1}(t)$ be stationary processes and consider taking $X_{2}(t)$ as a non-linear process of $X_{1}(t)$, for instance, $X_{2}(t)=X_{1}(t)X_{1}(t+u)$ for some predetermined time lag $u$. In this scenario, $f_{G_{2}G_{2}}(\lambda)$ accounts for the contribution of the non-linear terms $X_{2}(t)$ which cannot be measured by the linear term $X_{1}(t)$. The quantity $f_{G_{1}G_{1}}(\lambda)+f_{G_{2}G_{2}}(\lambda)$ can be interpreted as the spectral density, encompassing not solely the linear association of $X_{1}(t)$ but also the non-linear association of $X_{2}(t)$. This extension can be generalized to higher-order scenarios in a similar fashion. In the following section, we shall provide a derivation of our residual spectrum and elucidate the rationale behind its interpretation as aforementioned.

We consider the model

where $X_{0}(t)$ is a response variable, $\zeta$ is intercept, $\bm{X}(t):=(X_{1}(t),\ldots,X_{K}(t))^{\top}$ is a covariate process with autocovariance matrix $\Gamma_{\bm{X}}(u):={\rm E}\bm{X}(t)\bm{X}({t-u})$ such that, for any $i=1,\ldots,K$, $X_{i}(t):=X_{i}^{\prime}-{\rm E}X_{i}^{\prime}$ where $(X_{1}^{\prime}(t),\ldots,X_{K}^{\prime}(t))^{\top}$ a $s$-th order stationary process for any $s\in\mathbb{N}$, $\epsilon(t)$ is an i.i.d. centered disturbance process independent of $\{\bm{X}{(t)};t\in\mathbb{Z}\}$ with moments of all orders, and $\sum_{k=-\infty}^{\infty}k^{2}|b_{i}(k)|<\infty$ for all $i=1,\ldots,K$. The continuous version of this model was considered by Jenkins and Watts, (1968, Section 11.4.1, p.485).

For any random vector $\bm{W}(t):=(W_{1}(t),\cdots,W_{\ell}(t))^{\top}$ and $t_{1},\ldots,t_{\ell}\in\mathbb{Z}$, the cumulant of order $\ell$ of $(W_{1}(t_{1}),\cdots,W_{\ell}(t_{\ell}))$ is defined as

where the summation $\sum_{(\nu_{1},\ldots,\nu_{p})}$ extends over all partitions $(\nu_{1},\ldots,\nu_{p})$ of $\{1,2,\cdots,\ell\}$ (see Brillinger,, 1981, p.19). In this paper, we make the following assumption.

In conjunction with $\sum_{k=-\infty}^{\infty}k^{2}|b_{i}(k)|<\infty$ for all $i=1\ldots,K$ and the existence of moments of all orders for $\epsilon(t)$, Assumption 3.1 implies the cumulant summability (2) for any integer $\ell\geq 2$ and $(i_{1},\ldots,i_{\ell})\in\{0,1,\ldots,K\}^{\ell}$. The condition about the summability of cumulant is standard in this context, e.g., Brillinger, (1981, Assumption 2.6.2), Taniguchi et al., (1996, Assumption 2), Eichler, (2008, Assumption 3.1), Shao, (2010, Section 4), Jentsch and Pauly, (2015, Assumption 2.1), and Aue and Van Delft, (2020, Section 4). The sufficient condition for the summability of the cumulant is often discussed. The summability of the fourth-order cumulant under $\alpha$-mixing condition was shown by Andrews, (1991, Lemma 1). For univariate alpha-mixing processes, Neumann, (1996, Remark 3.1), Lee and Subba Rao, (2017, Remark 3.1), and Bücher et al., (2020, Lemma 2) showed the summability condition holds for higher-order cases. The essential tool of the proof is Theorem 3 of Statulevicius and Jakimavicius, (1988), which provides the upper bound of higher-order cumulant in terms of alpha-mixing coefficients and moments for univariate processes. Thus, the theorem is not applicable for multivariate processes. Besides, Panaretos and Tavakoli, 2013a (, Proposition 4.1) showed the MA($\infty$) process meets the summability condition. On the other hand, Kley, (2014, Lemmas 4.1 and 4.2) gave a sufficient condition for the summability of cumulant for indicator functions of the processes in a beautiful way (see also Kley et al., (2016, Proposition 3.1)). Employing his idea, we obtain the following lemma.

Therefore, the geometrically $\alpha$-mixing strictly stationary vector process with the appropriate moment condition holds under Assumption 3.1. We present two examples.

Under Assumption 3.1, the spectral density matrix $\bm{f}(\lambda):=\left(f_{ij}(\lambda)\right)_{i,j=0,\ldots,K}$ for the process $(X_{0}(t),X_{1}(t),\ldots,X_{K}(t))$ exists and twice continuously differentiable. For $i\in\{0,\ldots,K\}$, $f_{ii}(\lambda)$ is an auto-spectrum for the process $X_{i}(t)$ and, for $i_{1},i_{2}(\neq i_{2})\in\{0,\ldots,K\}$, $f_{i_{1}i_{2}}(\lambda)$ is a cross-spectrum for the processes $X_{i_{1}}(t)$ and $X_{i_{2}}(t)$.

The next theorem shows the model (1) has the following orthogonal representation.

From Theorem 3.1, the spectral density for $G_{j}$ is given by

By construction, $G_{j}$ can be interpreted as a process after eliminating the effect of the processes $G_{1},\ldots,G_{j-1}$. Thus, we refer to $f_{G_{j}G_{j}}$ as a residual spectral density of order $j$. These spectra are related to the important indexes.

When $K=1$, the squared coherence $\mathcal{C}_{1}(\lambda)$ is defined, for $\lambda\in[-\pi,\pi]$, as

which is the time series extension of the correlation, measures linear dependence between $\{X_{0}(t)\}$ and $\{X_{1}(t)\}$ (see, e.g., Brockwell and Davis, 1991, p.436).

In the case of $K=2$, $G_{2}$ has the spectral density of the form

In particular, when $X_{2}(t)$ is defined as the lagged process $X_{2}(t):=X_{1}(t)X_{1}(t-u)$, the lagged spectral coherence $\mathcal{C}_{2}(\lambda,u)$ is defined, for $\lambda\in[-\pi,\pi]$, as

The lagged coherence $\mathcal{C}_{2}(\lambda,u)$ measures the nonlinear (quadratic) dependence between $\{X_{0}(t)\}$ and $\{X_{1}(t)\}$. See details in Kedem-Kimelfeld, (1975). In this case, we also refer to $f_{G_{2}G_{2}}$ as a lagged spectral density.

As a higher-order extension of squared coherence, it is natural to define a squared coherence of order $d$ as

Since $f_{00}(\lambda)=\sum_{i=1}^{K}f_{G_{i}G_{i}}(\lambda)+{{\rm Var}(\epsilon(t))}/{(2\pi)}$, $\mathcal{C}_{d}(\lambda)\in(0,1)$ for any $\lambda\in[-\pi,\pi]$.

The residual spectral density is closely related to the partial cross spectrum, which was introduced by Dahlhaus, (2000), and is the cross spectrum after removing the linear effect of a certain process (see also Fiecas and Ombao, (2011, Section 2) and Eichler, (2008, p.969)). The partial cross spectrum of $X_{0}$ and $X_{2}$ given $X_{1}$ is defined by

and we can see the relationship between $f_{G_{2}G_{2}}$ and $f_{02|1}$:

Thus, we can interpret $f_{G_{2}G_{2}}$ as the rescaled squared partial cross spectrum of $X_{0}$ and $X_{2}$ given $X_{1}$.

We are interested in the existence of $X_{K}$ in the model (1), that is, $b_{K}(k)=a_{KK}(k)=0$ for any $k\in\mathbb{Z}$ under the presence of $X_{1},\ldots,X_{K-1}$. To this end, we consider the following hypothesis testing problem: the null hypothesis is

and the alternative hypothesis is $K_{0}:H_{0}$ does not hold. The null hypothesis $H_{0}$ holds true if and only if that $a_{KK}(k)=0$ for all $k\in\mathbb{Z}$ since $a_{KK}(k):=\int_{-\pi}^{\pi}A_{KK}\left(e^{-\mathrm{i}\lambda}\right)e^{-\mathrm{i}k\lambda}{\rm d}\lambda/(2\pi)$ and the proof of Theorem 3.1 shows $A_{KK}\left(e^{-\mathrm{i}\lambda}\right)=0$ if and only if $f_{G_{K}G_{K}}(\lambda)=0.$

Suppose the observed stretch $\{X_{j}(t);t=1,\ldots,n;j=0,\ldots,K\}$ is available. To estimate $\bm{f}(\lambda)$, we define the kernel spectral density estimator $\hat{\bm{f}}(\lambda):=\left(\hat{f}_{ij}(\lambda)\right)_{i,j=0,1,\ldots,K}$ and

where $M_{n}$ is a bandwidth parameter satisfies Assumption 4.1 (A1), $\omega(\cdot)$ is the lag window function defined as $\omega(x):=\int_{-\infty}^{\infty}W(t)e^{\mathrm{i}xt}\textrm{ d}t$, $W(\cdot)$ is the kernel function satisfy Assumption 4.1 (A2), for $h\in\{0,\ldots,n-1-u\}$,

and, for $h\in\{-n+1+u,\ldots,-1\}$,

with ${{X}_{i}{(.)}}=\sum_{t=1}^{n-u}{X}_{i}({t})/{{(n-u)}}$ for $i=0,1,\ldots,K$.

The condition (A1) is slight stronger than Eichler, (2008, Assumption 3.3 (iii)), that is, $n/M_{n}^{9/2}\to 0$ and $n/M_{n}^{2}\to\infty$ as $n\to\infty$. This is because that he dealt with the zero mean process. On the other hand, we need to estimate ${\rm E}X_{i}^{\prime}$ in the kernel density estimator. In the same reason, we assume the condition $\sum_{h=-n+1}^{n-1}\omega\left({h}/{M_{n}}\right)=O(M_{n})$ in (A2).

From the fact $A_{KK}\left(e^{-\mathrm{i}\lambda}\right)=0$ if and only if $f_{G_{K}G_{K}}(\lambda)=0$, together with the non-singularity of $\bm{f}_{K}(\lambda)$ and (7), the hypothesis (9) is equivalent to

and subsequently equivalent to

where

Our proposed test statistic is defined as

where