Estimating Time-Varying Effective Connectivity in High-Dimensional fMRI Data Using Regime-Switching Factor Models

TL;DR

This paper introduces a regime-switching factor model for high-dimensional fMRI data that captures both smooth and abrupt changes in brain connectivity, enabling more accurate detection of dynamic states and effective connectivity patterns.

Contribution

It proposes a novel three-step estimation method combining PCA, switching VAR models, and subspace estimates to analyze dynamic brain connectivity with high-dimensional data.

Findings

Outperforms K-means clustering in regime detection accuracy

Identifies distinct brain states with unique connectivity patterns

Reveals modular organization in resting-state networks

Abstract

Recent studies on analyzing dynamic brain connectivity rely on sliding-window analysis or time-varying coefficient models which are unable to capture both smooth and abrupt changes simultaneously. Emerging evidence suggests state-related changes in brain connectivity where dependence structure alternates between a finite number of latent states or regimes. Another challenge is inference of full-brain networks with large number of nodes. We employ a Markov-switching dynamic factor model in which the state-driven time-varying connectivity regimes of high-dimensional fMRI data are characterized by lower-dimensional common latent factors, following a regime-switching process. It enables a reliable, data-adaptive estimation of change-points of connectivity regimes and the massive dependencies associated with each regime. We consider the switching VAR to quantity the dynamic effective connectivity. We propose a three-step estimation procedure: (1) extracting the factors using principal component analysis (PCA) and (2) identifying dynamic connectivity states using the factor-based switching vector autoregressive (VAR) models in a state-space formulation using Kalman filter and expectation-maximization (EM) algorithm, and (3) constructing the high-dimensional connectivity metrics for each state based on subspace estimates. Simulation results show that our proposed estimator outperforms the K-means clustering of time-windowed coefficients, providing more accurate estimation of regime dynamics and connectivity metrics in high-dimensional settings. Applications to analyzing resting-state fMRI data identify dynamic changes in brain states during rest, and reveal distinct directed connectivity patterns and modular organization in resting-state networks across different states.

Equations

$${\mathbf{y}}_t=\mathbf{Q}{\mathbf{f}}_t+{\mathbf{ϵ}}_t$$

$${\mathbf{f}}_t={\mathbf{\Phi }}_{\mathbf{f}}(1){\mathbf{f}}_{t-1}+\dots +{\mathbf{\Phi }}_{\mathbf{f}}(P){\mathbf{f}}_{t-P}+{\mathbf{\eta }}_t,{\mathbf{\eta }}_t\sim N(0,{\mathbf{\Sigma }}_{\mathbf{\eta }})$$

$${\mathbf{y}}_t=\sum _{\ell =1}^P{\mathbf{\Phi }}_{\mathbf{y}}(\ell ){\mathbf{y}}_{t-\ell }+\mathbf{Q}{\mathbf{\eta }}_t$$

$${\mathbf{f}}_t=\sum _{\ell =1}^P{\mathbf{\Phi }}_{\mathbf{f}}^{[S_t]}(\ell ){\mathbf{f}}_{t-\ell }+{\mathbf{\eta }}_t,{\mathbf{\eta }}_t\sim N(0,{\mathbf{\Sigma }}_{\mathbf{\eta }}^{[S_t]})$$

$$z_{ij}=P(S_t=j|S_{t-1}=i)$$

$${\mathbf{y}}_t$$

$$=$$

$${\mathbf{y}}_t=\sum _{\ell =1}^P{\mathbf{\Phi }}_{\mathbf{y}}^{[S_t]}(\ell ){\mathbf{y}}_{t-\ell }+{\mathbf{\nu }}_t$$

$${\mathbf{F}}_t$$

$${\mathbf{y}}_t$$

$${\bf A}_{\bf F}^{[S_{t}]}=\left(\begin{array}[]{ccccc}{\bf\Phi}_{\bf f}^{[S_{t}]}(1)&{\bf\Phi}_{\bf f}^{[S_{t}]}(2)&\ldots&{\bf\Phi}_{\bf f}^{[S_{t}]}(P-1)&{\bf\Phi}_{\bf f}^{[S_{t}]}(P)\\ {\bf I}_{r}&{\bf 0}&\ldots&{\bf 0}&{\bf 0}\\ {\bf 0}&{\bf I}_{r}&\ldots&{\bf 0}&{\bf 0}\\ \vdots&&\ddots&&\vdots\\ {\bf 0}&{\bf 0}&\ldots&{\bf I}_{r}&{\bf 0}\\ \end{array}\right).$$

$$\hat{r}=\underset{\{1,\dots ,L_r\}}{\mathrm{argmax}}\left\{\ln \left(\frac{1}{NT}\sum _{t=1}^T\Vert {\widehat{\mathbf{ϵ}}}_t(r){\Vert }_2^2\right)+r\left(\frac{N+T}{NT}\right)\ln \left(\frac{NT}{N+T}\right)\right\}$$

$${\mathbf{F}}_{t|t}^j$$

$$V_{t|t}^j$$

$${\mathbf{X}}_{t|T}^j$$

$$V_{t|T}^j$$

$$V_{t,t-1|T}^j$$

$${\mathbf{M}}_{t|t}^j$$

$${\mathbf{M}}_{t|T}^j$$

$$P_t$$

$$P_{t,t-1}$$

$${\hat{\mathbf{A}}}_{\mathbf{F}}^{[j]}$$

$${\mathbf{\Sigma }}_{\mathbf{w}}^{[j]}$$

$${\hat{z}}_{ij}$$

$$\sqrt{T}(\tilde{\mathbf{b}}-\mathbf{b})\stackrel{D}{\to }N(0,\mathbf{G})$$

Full text

Keywords: Regime-switching models, Large VAR models, Factor analysis; Principal components analysis; Dynamic Brain Connectivity.

1 Introduction

Most analyses of functional connectivity (FC) using functional magnetic resonance imaging (fMRI) data implicitly assumed that relationships between distinct brain regions are static (stationary) across time. Time-invariant FC metrics such as correlations between fMRI time series are computed over the entire period of recording. Recent years have seen increased interest in investigating dynamic changes in FC patterns over time, often referred to as dynamic (time-varying) functional connectivity [Hutchison et al., 2013a, Calhoun et al., 2014]. Several studies have reported temporal fluctuations in FC at time-scales of seconds to minutes, in both strength and directionality of the connections, even during resting state [Chang and Glover, 2010, Allen et al., 2012, Leonardi et al., 2013, Hutchison et al., 2013b, Zalesky et al., 2014].

The simplest and most common approach to examining the dynamic behavior in connectivity is the sliding-window correlation, which involves computing locally stationary correlations over consecutive windowed short-time segments of data to produce time-varying FC metrics [Chang and Glover, 2010, Allen et al., 2012, Hutchison et al., 2013a, Zalesky et al., 2014]. However, this approach is limited by the choice of optimal window length: a long window has low statistical power to detect abrupt and highly localized changes, while a short window produces noisy estimates for smooth changes. An alternative strategy is the model-based approach, which can provide a unified, parsimonious framework to characterize the dynamic connectivity structure based on the time-dependent model parameters. For example, the time-varying multivariate volatility models [Lindquist et al., 2014] and the time-varying vector autoregressive (VAR) models [Havlicek et al., 2010, Samdin et al., 2015] have been used to capture effectively instantaneous temporal changes in fMRI-based functional and effective connectivity (a more specific cross-dependence with directionality, in a sense that it measures the causal influence of one brain region on another).

Recent evidence from fMRI studies suggested state-related types of dynamics in FC: time-varying but reoccurring connectivity patterns which switches according to a few discrete underlying quasi-stable brain states (regimes). This non-stationarity is characterized by rapid transitions between regimes and smooth changes within a regime. Various analytical approaches have been used to identify these replicable dynamic ‘connectivity states’. These include K-mean clustering of the windowed correlations [Allen et al., 2012, Hutchison et al., 2013a], which, however, ignores information about the temporal order of the dynamics, hidden Markov models producing the state-time alignments [Baker et al., 2014], or algorithms to detect change points in connectivity [Cribben et al., 2012, Jeong et al., 2016]. However, these studies focused on evaluating the un-directed connectivity. Our recent work [Samdin et al., 2016] proposed a more general method based on the switching VAR (SVAR) models to infer dynamic states of effective connectivity in fMRI and EEG data. The next challenge is to estimate the high-dimensional connectivity states for a large number of brain regions, where traditional analyses based on pair-wise correlations or using average signals from parcellated regions of interest (ROIs) might produce sub-optimal results.

In this paper, we propose a new approach based on regime-switching factor models for estimating temporal changes in effective connectivity states in high-dimensional fMRI data for a whole-brain network analysis. Precisely, our approach is to first employ a factor analysis model to characterize the large fMRI data via a small number of common, latent (unobservable) factor components, and then identify the dynamic connectivity regimes based on these low-dimensional summary signals. We develop a non-stationary factor model which takes into account the time-variation of the underlying serial cross-correlation structure of the high-dimensional data, by introducing regime-switching in the factor dynamics, By specifying the factors to evolve as a Markov-switching VAR process, we derive a factor SVAR model for the observation space, which is an extension of the SVAR model used in [Samdin et al., 2016] to the large-dimensional case. Such formulation implies projection of the high-dimensional directed connectivity matrix onto a lower-dimensional subspace (small VAR coefficient matrix of factors, spanned by the factor loadings), and thus allows us to capture the changes in connectivity regimes in these subspaces driven by the few factors. It enables a reliable and computationally-efficient estimation of the regime change-points and the massive dependence measures associated with each regime.

We develop a three-step estimation procedure. The first step is initial estimation of connectivity subspace shared across regimes based on a stationary factor model. The number of factors and a common factor loading (specifying the dimension and the span of the subspace) are estimated by applying the principal component analysis (PCA) on the data. The second step is the dynamic regime segmentation based on the factor SVAR model formulated in a state-space form. The change-points between connectivity states are identified via switching Kalman filter and switching Kalman smoother (SKF and SKS). Simultaneously, the regime-dependent parameters of the latent switching factor process are updated by the expectation-maximization (EM) algorithm, with the common factor loadings initialized and fixed using the estimates from the first step. Then, the fMRI signals are partitioned according to the estimated states, and fitted with a separate factor model for each regime to obtain regime-dependent factor loadings. The third step is estimation of within-regime connectivity metrics, where the estimates of the high-dimensional VAR connectivity graph/matrix for each state are constructed from the low-dimensional factor subspace parameters estimated from the previous two steps. We evaluated the performance of our method via simulations by comparing with the K-mean clustering approach. Application to the resting-state fMRI data reveals switching states of the resting state connectivity networks, with the modular organization changes across different states.

2 Regime-switching Factor VAR Models

In this section, we first describe the stationary factor model with an autoregressive factor process. Then, we introduce a non-stationary generalization with regime-switching in the factor dynamics.

2.1 The Factor Model

Let ${\bf y}_{t}=[y_{1t},\ldots,y_{Nt}]^{\prime}$ be a $N\times 1$ observed vector of non-stationary time series of fMRI at time points $t=1,\ldots,T$. The cross-section dimension of the time series $N$ can be comparable to or even larger than the sample size $T$ (or length of time series). We suppose the high-dimensional time series is driven by a small number of latent factors. Specifically, we consider a factor model defined by

[TABLE]

where ${\bf f}_{t}=[f_{t1},\ldots,f_{tr}]^{\prime}$ is a $r\times 1$ vector of unobserved common factors with mean zero and covariance matrix $\boldsymbol{\Sigma}_{\bf f}$, ${\bf Q}=[{\bf q}_{1},\ldots,{\bf q}_{r}]$ is a $N\times r$ constant factor loading matrix assumed to be orthonormal, i.e. ${\bf Q}^{\prime}{\bf Q}={\bf I}_{r}$ where ${\bf I}_{r}$ denotes a $r\times r$ identity matrix, $r$ is the number of factors satisfying $r<<N$, and ${\boldsymbol{\epsilon}}_{t}=[\epsilon_{t1},\ldots,\epsilon_{tN}]^{\prime}$ is $N\times 1$ vector of noise components with mean zero and covariance matrix $\boldsymbol{\Sigma}_{\boldsymbol{\epsilon}}={\mbox{diag}}({\sigma}_{\epsilon_{1}}^{2},\ldots,{\sigma}_{\epsilon_{N}}^{2})$, assuming the error terms are cross-sectionally uncorrelated. The model captures the correlation between the time series via the mixing of some common factors ${\bf f}_{t}$ by ${\bf Q}$. The model (1) allows for dimension-reduction in the sense that the serial and cross-correlation in the high-dimensional observational process $\{{\bf y}_{t}\}$ is driven by the much lower-dimensional factor process $\{{\bf f}_{t}\}$ and mixing matrix ${\bf Q}$.

The evolution of the latent factor dynamics in $\{{\bf f}_{t}\}$ can be modeled by a stationary vector autoregressive (VAR) process of order $P$, VAR($P$)

[TABLE]

where ${\bf\Phi}_{\bf f}(\ell)$ is the $r\times r$ AR coefficients matrix at lag $\ell$ for $l=1,\ldots,P$ and ${\boldsymbol{\eta}}_{t}$ is a $r\times 1$ Gaussian white noise process with mean zero and covariance matrix $\boldsymbol{\Sigma}_{\boldsymbol{\eta}}$. Both processes $\{{\bf y}_{t}\}$ and $\{{\bf f}_{t}\}$ are non-stationary where the factor loadings ${\bf Q}$ and the AR coefficients matrices for factors ${\bf\Phi}_{\bf f}(\ell)$ are time-constant.

Factor VAR Model: The temporal inter-dependence in the high-dimensional observation process ${\bf y}(t)$ can be characterized by the much lower-dimensional VAR process of ${\bf f}_{t}$ in (2). This forms the basic idea of our recent works [Ting et al., 2014, Wang et al., 2016], where we developed a factor-based VAR (f-VAR) model for the observations ${\bf y}(t)$ by substituting (2) into (1) and assuming ${\boldsymbol{\epsilon}}_{t}$ approximately zero, which gives

[TABLE]

where ${\bf\Phi}_{\bf y}(\ell)={\bf Q}{\bf\Phi}_{\bf f}(\ell){\bf Q}^{\prime}$ are high-dimensional $N\times N$ coefficients matrices for ${\bf y}_{t}$, an orthogonal projection of the smaller matrices ${\bf\Phi}_{\bf f}(\ell)$ on to lower-dimensional subspace spanned by the columns of ${\bf Q}$. It provides a low-rank approximation for the dependence structure in ${\bf y}(t)$. The model subspace can be learned by using the principal component analysis (PCA) where the estimator for ${\bf Q}$ are defined by eigenvectors corresponding to the $r$ largest eigenvalues of the sample covariance matrix of ${\bf y}(t)$. It leads to substantially improved consistency and computational efficiency in estimating VAR models under high-dimensional settings, compared to the traditional least-squares estimator as shown in [Ting et al., 2014].

The coefficients matrix ${\bf\Phi}_{\bf y}(\ell)$ can quantify directed interactions in a network with large number of nodes (e.g. a large-scale network of brain regions) at time lag $\ell$. There exists a directed influence in the Granger-causality sense with direction from node $j$ to node $i$ for any connection strength $|{\bf\Phi}_{ij}|>0$, where ${\bf\Phi}_{ij}$ is the ($i,j$)-th element of ${\bf\Phi}_{\bf y}$. When applied to identify effective brain connectivity networks with a large number of nodes from resting-state fMRI data [Ting et al., 2014], the estimates provided more reliable interpretation and capable of revealing the modular, hierarchical structure of the brain networks during rest, by varying the subspace dimension $r$.

2.2 Regime-switching in Factor Dynamics

Factor Model with Regime-Switching: We now generalize the stationary factor model in (1) to allow for time-variation in the serial interdependence structure of the latent factors, by introducing regime-switching in the coefficient matrices of the VAR factor specification in (2). In this respect, we propose a non-stationary factor model with regime-switching factor dynamics. More precisely, we assume the factor loadings remain stationary but the factors to follow a Markov-switching VAR (SVAR) process of order $P$, SVAR($P$). This class of models has been applied for modeling of econometric data [Krolzig, 2013]. The SVAR is a quasi-stationary model consisting of a set of $K$ independent VAR models, each indexed by a hidden random indicator $S_{t}$

[TABLE]

here $\{S_{t}\in\{j=1,\ldots,K\},\ t=1,\ldots,T\}$ is a sequence of state/regime variables, which is time-dependent and take values in a discrete space $j=1,\ldots,K$; and $\{{\bf\Phi}_{\bf f}^{[j]}(\ell),\ell=1,\ldots,P\}$ are coefficient matrices for state $j$. This is a generalized version of (2) which allows for structural changes in the VAR coefficients. The AR coefficients matrices ${\bf\Phi}_{\bf f}^{[S_{t}]}(\ell)$ are piecewise constant function of the discrete state $S_{t}$, i.e. constant within time-blocks belong to a same regime but change across different regimes. This renders the factor process piecewise stationary, a special form of departure from stationarity. However, the proposed model differs from the classical piecewise constant processes (e.g., piecewise VAR processes) primarily because in the classical piecewise processes the future blocks are not at all related to previous blocks. Our proposed model permits recurring regimes where future blocks could be related to past blocks if they were both indexed by the same state. This has important implications in estimation and inference because we can pool together different time-blocks of the same regime (that are indexed by the same state) thus producing more accurate and more efficient estimates

We assume that $S_{t}$ follows a $K$-state first-order Markovian process with a $K\times K$ transition matrix $Z=[z_{ij}],1\leq i,j\leq K$ where

[TABLE]

denotes the probability of transition from state $i$ at time $t-1$ to state $j$ at $t$. Only one latent process (and hence only one VAR process) is “active” (or turned on) at each time point $t$. The remaining latent processes are turned off. This allows recurring changes in the temporal interdependence structure of the factors as characterized by ${\bf\Phi}_{\bf f}^{[S_{t}]}$, which switches over time between the finite number of regimes, according to the regime indicator $S_{t}$ at time $t$. Compared to using a switching VAR model directly on ${\bf y}_{t}$, the specification of Equation (4) allows us to detect the change-points of the high-dimensional dependence structure based on a small number of factor series. We denote the model parameters by $\boldsymbol{\theta}=\{\boldsymbol{\theta}_{j}=({\bf\Phi}_{\bf f}^{[j]},\boldsymbol{\Sigma}_{\boldsymbol{\eta}}^{[j]}):j\in\{1,\ldots,K\}\}$ which are assumed unknown and to be estimated.

Factor Switching-VAR Model: We shall derive a high-dimensional switching VAR model from the non-stationary factor model with a regime-switching autoregressive factor process as defined by Equation (1) and (4). The regime-switching in the high-dimensional interdependence structure of observations $\{{\bf y}_{t}\}$ can be driven by that of the lower-dimensional SVAR factor process in (4). Substituting (4) into (1) yields

[TABLE]

Finally, we have a factor-based Markov-switching VAR (f-SVAR) for ${\bf y}_{t}$

[TABLE]

where ${\bf\Phi}_{\bf y}^{[S_{t}]}(\ell)={\bf Q}{\bf\Phi}_{\bf f}^{[S_{t}]}(\ell){\bf Q}^{\prime}$ and ${\boldsymbol{\nu}}_{t}={\bf Q}{\boldsymbol{\eta}}_{t}$. The model is a nonstationary generalization of the factor VAR model in (3), by allowing a regime-switching in the coefficient parameters. It provides a tool to capture the regime-switching in the large $N$-dimensional serial inter-dependence structure via a low-dimensional space. The model can quantify dynamics of a large-scale directed network with state-dependent changes in the network structure, i.e. switching according to distinct states. It enable the detection of the temporal change points in the network dependency structure, as well as estimation of the directed dependencies between massive number of nodes associated with each state.

State-Space Formulation: We propose a state-space representation for the factor model with regime-switching factors, to enable sequential estimation in time of the latent factors and the switching states. The latent switching VAR factor process (4) forms the state-equation which is projected to the high-dimensional space using the factor model (1) with an error as the observation equation. Defining the dynamic factor structure ${\bf F}_{t}=[{\bf f}^{\prime}_{t},{\bf f}^{\prime}_{t-1},\ldots,{\bf f}^{\prime}_{t-P+1}]^{\prime}$ as state vector, the model (1) and (4) are formulated in a switching linear Gaussian state-space form [Kim, 1994]

[TABLE]

The SVAR($P$) factor process (4) is re-written in a SVAR($1$) form of (10) in the state equation, where $\mathbf{w}_{t}=[{\boldsymbol{\eta}}^{\prime}_{t},{\bf 0}^{\prime},\ldots,{\bf 0}^{\prime}]^{\prime}$ is $rP\times 1$ state noise, and ${\bf A}_{\bf F}^{[S_{t}]}$ is a $rP\times rP$ state transition matrix switching with the state variables $S_{t}$, and of the form

[TABLE]

The matrix ${\bf A}_{\bf F}^{[S_{t}]}$ describes the directed connectivity that varies across states. The unobserved SVAR($P$) dynamic factors ${\bf F}_{t}$ now follows a (higher dimensional) latent SVAR(1) process. A noisy version of factor model is re-formulated from (1) as in the observation equation (11) by introducing an idiosyncratic noise $\boldsymbol{\epsilon}_{t}$, and with a $r\times rP$ mapping matrix $\mathbf{H}=[{\bf Q},{\bf 0},\ldots,{\bf 0}]$. We assume both $\{{\boldsymbol{\epsilon}}_{t}\}$ and $\{\mathbf{w}_{t}\}$ are white Gaussian noise, $\boldsymbol{\epsilon}_{t}\sim N(\mathbf{0},\boldsymbol{\Sigma}_{\boldsymbol{\epsilon}})$ and $\mathbf{w}_{t}\sim N(\mathbf{0},\boldsymbol{\Sigma}_{\bf w}^{[S_{t}]})$, with a time-constant state noise covariance matrices $\boldsymbol{\Sigma}_{\boldsymbol{\epsilon}}$ and the $\boldsymbol{\Sigma}_{\bf w}^{[S_{t}]}$ switching with $S_{t}$. Both the factor loadings ${\bf Q}$ and the noise covariance matrix in the observation equation are assumed to be regime-invariant and shared across regimes. The processes $\{{\bf f}_{t}\}$ and $\{{\boldsymbol{\epsilon}}_{t}\}$ are uncorrelated. Instead of hard state assignment for each time-point $t$, we can evaluate the probability of activation for each state, $P(S_{t}=j|{\bf y}_{1:T})$, which is termed “soft-alignment”. We denote all model parameters from each of the states as $\boldsymbol{\Theta}=\{\boldsymbol{\Theta}_{j}=({\bf A}_{\bf F}^{[j]},\boldsymbol{\Sigma}_{\bf w}^{[j]}):j\in\{1,\ldots,K\}\}$.

3 Estimation

We develop a three-step procedure for efficiently estimating the dynamic connectivity states in the high-dimensional fMRI data based on the proposed non-stationary factor model with regime-switching. In the first step, we explore the connectivity subspace assumed as common and shared across regimes, by fitting a stationary factor model (1) to the entire fMRI time series. We apply the method of PCA to estimate the factor loadings ${\bf Q}$ and the latent factors ${\bf f}_{t}$, and the Bayesian information criterion (BIC) of [Bai and Ng, 2002] to select the optimal number of factors. In the second step, we perform connectivity regime segmentation in the low-dimensional subspace relying on a factor model with Markov-switching VAR factor process (4). Based on the state-space representation (10)-(11), the latent factor process can be jointly estimated conditioned on the observational factor model. The temporal change-points of the regimes can be detected via the estimated state sequence $\{\widehat{S}_{t}\}$ by the SKF and SKS, and the factor VAR coefficient matrix ${\bf\Phi}_{\bf f}^{[j]}$ for each regime is updated iteratively using the EM algorithm. In the third step, we estimate the regime-dependent high-dimensional connectivity matrix ${\bf\Phi}_{\bf y}^{[j]}$ for the observation space using the estimated subspace parameters from the first two steps.

3.1 Step 1: Estimation of a Common Factor Model

PCA is a common approach to estimating approximate factor model based on the eigen-decomposition of sample covariance matrix [Bai, 2003, Stock and Watson, 2002]. Let ${\bf V}_{1}$, $\ldots$, ${\bf V}_{N}$ be $N$ orthonormal eigenvectors corresponding to the eignevalues of the $N\times N$ sample covariance matrix ${\bf S}_{\bf y}=\sum_{t=1}^{T}{\bf y}_{t}{\bf y}^{\prime}_{t}$, in a decreasing order such that $\widehat{\lambda}_{1}\geq\ldots\geq\widehat{\lambda}_{N}>0$. The PCA estimator of the loadings $\widehat{\bf Q}=\left[{\bf V}_{1},\ldots,{\bf V}_{r}\right]$ is defined by a matrix whose columns are the $r$ orthonormal eigenvectors corresponding to the largest $r$ eignevalues, and the factors can be estimated by $\widehat{\bf f}_{t}=\widehat{\bf Q}^{\prime}{\bf y}_{t}$. [Bai, 2003] has showed that the PCA estimators are consistent and asymptotically normal, under settings of large $N$ and large $T$. Besides, the estimates can be computed efficiently even under situations when $N<T$ (on the small $T\times T$ temporal covariance matrix instead of the huge $N\times N$ spatial sample covariance ${\bf S}_{\bf y}$). We can compute the noise covariance estimator as $\widehat{\boldsymbol{\Sigma}}_{\boldsymbol{\epsilon}}=\sum_{t=1}^{T}\widehat{\boldsymbol{\epsilon}}_{t}\widehat{\boldsymbol{\epsilon}}^{\prime}_{t}$ based on the residuals $\widehat{\boldsymbol{\epsilon}}_{t}={\bf y}_{t}-\widehat{\bf Q}\widehat{\bf f}_{t}$ from the fitted factor model. We fit an VAR model (2) to the estimated factors $\{\widehat{\bf f}_{t}\}$ by the least-squares (LS) method, and obtain the AR coefficient estimates $\widehat{\bf\Phi}_{\bf f}(\ell)$. For PCA estimation, the number of factors can be determined by model selection using BIC

[TABLE]

where $\|{\bf x}\|$ denotes the Euclidean norm of a vector ${\bf x}$ and $L_{r}$ is a bounded integer such that $r\leq L_{r}$.

3.2 Step 2: Estimation of Regime-switching Factor Model

Based on the state-space formulation (10)-(11), the objective is to extract the underlying states $\{\widehat{S}_{t}\}$, and to estimate the unknown coefficient matrix ${\bf\Phi}_{\bf f}^{[j]}$ and factor signals in ${\bf F}_{t}=[{\bf f}^{\prime}_{t},{\bf f}^{\prime}_{t-1},\ldots,{\bf f}^{\prime}_{t-P+1}]^{\prime}$ of the latent Markov-SVAR factor process given observations $\widehat{\bf y}_{t},t=1,\ldots T$.

Filtering and Smoothing: The inference of $S_{t}$ and ${\bf F}_{t}$ involve computing, sequentially in time, the filtered probabilities $Pr(S_{t}|{\bf y}_{1:t})$ and the filtered densities $p({\bf F}_{t}|{\bf y}_{1:t})$, given the available signal observations up to time $t$, ${\bf y}_{1:t}=\{{\bf y}_{1},\ldots,{\bf y}_{t}\}$, and the more accurate smoothed probabilities $P(S_{t}|{\bf y}_{1:T})$ and densities $p({\bf F}_{t}|{\bf y}_{1:T})$ given the available entire set of observations ${\bf y}_{1:T}=\{{\bf y}_{1},\ldots,{\bf y}_{T}\}$. We estimate the filtered and smoothed densities of ${\bf F}_{t}$ given state $j$ at time $t$, by the KF and the KS, respectively

[TABLE]

where ${\bf F}_{t|t}^{j}$ and $V_{t|t}^{j}$ are mean and covariance of the filtered density $p({\bf F}_{t}|{\bf y}_{1:t},S_{t}=j)$, ${\bf F}_{t|T}^{j}$ and $V_{t|T}^{j}$ are mean and covariance of the smoothed density $p({\bf F}_{t}|{\bf y}_{1:T},S_{t}=j)$ given state $j$ at time $t$, and $V_{t,t-1|T}^{j}$ is the cross-variance of joint density $p({\bf F}_{t},{\bf F}_{t-1}|{\bf y}_{1:T},S_{t}=j)$. The estimates of filtered and smoothed state occupancy probability of being state $j$ at time $t$ are also computed as

[TABLE]

EM Estimation: The estimates of the factor-subspace dynamic parameters in ${\bf A}_{\bf F}^{[j]}$ and $\boldsymbol{\Sigma}_{\bf w}^{[j]}$ can be obtained by the maximum likelihood (ML) method by maximizing the log-likelihood ${L}=\log p({\bf y}_{1:T}|\boldsymbol{\Theta})$ with respect to each parameter. Here, we use the EM algorithm for the switching state-space model suggested by [Murphy, 1998]. In the expectation step (E-step), the sufficient statistics are obtained from the smoothed estimates

[TABLE]

where ${\bf F}_{t|T}$, $V_{t|T}$ and $V_{t,t-1|T}$ are quantities of the smoothed densities $p({\bf F}_{t}|{\bf y}_{1:T})$ and $p({\bf F}_{t},{\bf F}_{t-1}|{\bf y}_{1:T})$, corresponding to (15) to (17) by marginalizing out the state variable $j$ of the $p({\bf F}_{t}|{\bf y}_{1:T},S_{t}=j)$ and $p({\bf F}_{t},{\bf F}_{t-1}|{\bf y}_{1:T},S_{t}=j)$ using Gaussian approximation. We retain the terms switching KF (SKF) and switching KS (SKS) to refer to KF/KS approach to estimating state parameters of the SVAR model, as in [Murphy, 1998].

In the maximization step (M-step), the estimates of the model parameters for regime $j$ are updated as follows

[TABLE]

where the weights $W_{t}^{j}=\mathbf{M}_{t|T}^{j}$ are computed from the smoothing step. The model parameters are iteratively until some convergence criteria are satisfied, to produce the ML estimates $\boldsymbol{\Theta}^{*}$. We used randomized initial estimates for entries of $\widehat{\bf A}_{\bf F}^{[j]}$. The factor loading matrix ${\bf Q}$ in $\mathbf{H}$ and the noise covariance $\boldsymbol{\Sigma}_{\boldsymbol{\epsilon}}^{[j]}$ which are assumed common to all regimes, remain fixed with the PCA estimates from Step 1, and not updated by the EM algorithm. Note that here the regime estimation is done based on the state equation of low-dimensional factors. This will lead to substantial computational reduction, and improve the identifiability of the individual subspace parameter estimators.

Regime Segmentation: Given the EM-estimated model parameters $\boldsymbol{\Theta}^{*}$, the reliminary temporal regime segmentation $\widehat{\bf\Phi}_{\bf f}^{[\widehat{S}_{t}]}$ in the subspace is defined by the latent state sequence estimated using the SKF, $\widehat{S}^{\text{SKF}}_{t}=\operatorname*{arg\,max}_{j}P(S_{t}=j|{\bf y}_{1:t})$ in (18) which indicates the most likely active state for each time point. This is then further refined by the SKS, $\widehat{S}^{\text{SKS}}_{t}=\operatorname*{arg\,max}_{j}P(S_{t}=j|{\bf y}_{1:T})$ in (19) based on both the past and future observations. We can also utilize this state-time alignment provided in $\widehat{S}^{\text{SKF}}_{t}$ and $\widehat{S}^{\text{SKS}}_{t}$ to partition the observed fMRI signals into their corresponding states, and the time-segments of each regime is then fitted with a separate stationary factor model to derive state-dependent estimators, as described in the next step.

3.3 Step 3: Estimation of Regime-dependent Connectivity Matrices

We investigate two different schemes for constructing the estimators for the high-dimensional VAR-based connectivity matrix or graph for each state $\widehat{\bf\Phi}_{\bf y}^{[j]}$, by plugging in the subspace parameter estimators obtained in the first two steps: (1) Coupled SVAR estimator (with common factor loadings) $\widehat{\bf\Phi}_{\bf y}^{[j]}(\ell)=\widehat{\bf Q}\widehat{\bf\Phi}_{\bf f}^{[j]*}(\ell)\widehat{\bf Q}^{\prime}$, by substituting in the f-SVAR model in (9) with the EM estimate $\widehat{\bf\Phi}_{\bf f}^{[j]*}(\ell)$ from Step 2 and the PCA estimate $\widehat{\bf Q}$ from Step 1. Note that conditioned on a common factor loading $\widehat{\bf Q}$, the factor coefficient matrices $\widehat{\bf\Phi}_{\bf f}^{[j]*}(\ell)$ of all regimes are jointly estimated by the EM, weighted at each state by the smoothed state occupancy probability $P(S_{t}=j|{\bf y}_{1:T})$ in (19). (2) Decoupled SVAR estimator (with state-dependent factor loadings) $\widetilde{\bf\Phi}_{\bf y}^{[j]}(\ell)=\widetilde{\bf Q}^{[j]}\widetilde{\bf\Phi}_{\bf f}^{[j]}(\ell)\widetilde{\bf Q}^{\prime[j]}$ by substituting in a separate f-VAR model in (3) for each state. ($\widetilde{\bf\Phi}_{\bf f}^{[j]}$, $\widetilde{\bf Q}$) are PCA estimates by fitting distinct stationary factor models (1) separately to each of the regime time-courses, derived from the SKS segmentation in Step 2. The limiting distribution of the factor-VAR estimator has been derived in [Ting et al., 2014] (Theorem 2). For ease of exposition, we drop the state index $j$ and focus on the VAR(1). The subspace estimator $\widetilde{\bf b}=vec(\widetilde{\bf\Phi}_{\bf y})=vec(\widetilde{\bf Q}\widetilde{\bf\Phi}_{\bf f}\widetilde{\bf Q}^{\prime})$ has an asymptotic normal distribution as $T\rightarrow\infty$

[TABLE]

where ${\bf G}=({\bf Q}\boldsymbol{\Sigma}_{\boldsymbol{\eta}}{\bf Q}^{\prime})\otimes({\bf Q}\boldsymbol{\Gamma}_{\bf f}{\bf Q}^{\prime})$ with ${\Gamma}_{\bf f}=cov({\bf f}_{t})$ and $\otimes$ denotes the Kronecker product. By replacing with the PCA estimates, the covariance matrix of the estimator can be estimated, defined by $\widehat{\bf G}$. Based on this, we test the significance of each subspace VAR coefficient in $\widetilde{\bf b}$ as being different from zero, with $H_{0}:b_{k}=0$ against $H_{1}:b_{k}\neq 0$, where $b_{k}$ is $k$-th element of ${\bf b}$. The test statistic is approximately distributed as $t_{k}=\widehat{b}_{k}/{\sqrt{\widehat{\bf G}_{kk}/T}}\sim N(0,1)$ when $T$ is sufficiently large, where $\widehat{\bf G}_{kk}$ is $k$-th diagonal entry of $\widehat{\bf G}$. A coefficient is significant if the $p-value<\alpha/D$ with $\alpha$ the significance level and $D=N^{2}$ the number of tested coefficients, implying corrections for multiple testing by Bonferroni method.

4 Simulations

In this section, we evaluate the numerical performance of the proposed factor-SVAR model-based estimators in identifying state-dependent changes in large-scale directed connectivity networks through simulations. The objective is to measure the ability of our estimation procedures in (1.) detecting the change-points of connectivity regimes via the estimated state sequence, and (2) estimating the high-dimensional directed connectivity matrix or graph between nodes for each regime.

Data Generation: We generated data from a regime-switching VAR(1) model with with $K=2$ states, with different coefficient matrix of the independent VAR for each state to characterize distinct connectivity patterns. To emulate the modular connectivity network structure, we assume a block-diagonal VAR coefficient matrix, formed by $10\times 10$ dimensional non-zero sub-blocks along the main diagonal. Each sub-block represents the directed connectivity in a sub-network of 10 nodes. The entries of the sub-blocks were randomly drawn from a uniform distribution. Here, we set the two state-dependent coefficient matrices with distinct structure as ${\bf\Phi}_{\bf y}^{[1]}:a_{ij}\sim\mbox{U}[-0.4~{}0.4]$ and ${\bf\Phi}_{\bf y}^{[2]}:a_{ij}\sim\mbox{U}[-0.2~{}0.2]$, for $i$ and $j$ in the same block. The entries of the off-diagonal blocks are zero. We set the same noise covariance matrix for both state $\boldsymbol{\Sigma}_{\boldsymbol{\eta}}=0.5{\bf I}$.

Locally-stationary time-series data with piece-wise stable connectivity structure over time, were obtained by concatenating the two VAR processes simulated independently. The simulated data consists of four time-blocks each from a VAR and of fixed length $T_{B}=50$ (total length of $T=200$), with 3 change points at times $t=50$, $t=100$ and $t=150$. The sample size available for the VAR model of each state is only $T=100$. To emulate the state-dependent recurring changes in the VAR connectivity structure, the successive time-blocks were generated according to the distinct coefficient matrices in a cyclic manner, alternating between the two connectivity states, following procedure in [Monti et al., 2014] for functional connectivity. Thus, the state labels and the corresponding state-dependent VAR coefficient matrices for each time points are considered known and used as ground-truth for evaluation, i.e. ($S_{t}=1$, ${\bf\Phi}_{\bf y}^{[1]}$) for $t=1,\ldots,50$ and $t=101,\ldots,150$; ($S_{t}=2$, ${\bf\Phi}_{\bf y}^{[2]}$) for $t=51,\ldots,100$ and $t=151,\ldots,200$.

We investigate the impact of increasing network dimensions on the estimation performance in terms of accuracy and consistency, by varying $N$ from 10 to 100 with an increment of 10 or one sub-block. The sample size $T$ is fixed to create the scenarios of dimensionality $N<T$ and $N\approx T$. The simulations were repeated 100 times. We computed factor-SVAR model-based estimates for the state sequence $\widehat{S}_{t},t=1,\ldots,200$ and the coefficient matrix for each state $\widehat{\bf\Phi}_{\bf y}^{[j]}$, using the estimation steps in Section 2. The number of factors was selected adaptively for each simulated data using BIC in 12.

Benchmark with K-means Clustering: We compare the performance of our factor-SVAR estimator with an recent approach based on K-means clustering of time-variant VAR coefficients proposed by [Samdin et al., 2016]. Here, a sliding window is first employed to estimate the time-evolving directed connectivity, by fitting stationary VAR model to shifted short-time windows of fixed length to obtain time-dependent estimates of VAR coefficients matrices. We used a rectangular window with a bandwidth of 30 samples and shift of 1 sample. The relatively short-segments may render the traditional ordinary least-squares (LS) fits of large-dimensional VAR matrices inaccurate, due to insufficient information to estimate the huge number of parameters. Therefore, we used the $L_{2}$-regularized or ridge estimator which imposes a $L_{2}$ norm penalty on the AR coefficients in the LS regression, to obtain a better-conditioned estimate particularly in high-dimensional settings. The regularization parameter was set $\lambda=0.1$, as suggested by [Korobilis, 2013] for VAR model estimation. Then, the K-means clustering algorithm is applied to the estimated time-variant VAR (TV-VAR) coefficients to partition the dynamic connectivity structure into the distinct states or regimes. As in [Allen et al., 2012], we used the L1 (Manhattan) distance which may be more effective for clustering high-dimensional data, compared to the L2 (Euclidean) distance.

Our proposed SVAR approach has more advantages than the K-means clustering of time-variant VAR coefficients, as discussed in [Samdin et al., 2016]. First, the sliding-window approach is limited by the choice of window size which is crucial: a large window leads to low statistical power for detecting abrupt and highly localized changes; a small window produces noisy estimates for smooth changes. In contrast, the SVAR model is capable of detecting changes at different time scales, both smooth and abrupt, avoiding the problems associated with fixed time windowing. Second, the K-means algorithm provides a ‘hard’ assignment of time points into states and does not account for the temporal correlation structure. In contrast, the SVAR estimator generates ‘soft’ state-time alignment by estimating sequentially, for each time point, the probability of the occupying states based on the entire observation time course.

Performance Measure: To measure the performance of the estimated VAR connectivity graphs within each regime, we computed for each simulation the total squared errors over all entries between the ground-truth and the estimators of the VAR coefficient matrix for each state $j=1,2$, $\|\widehat{\bf\Phi}_{\bf y}^{[j]}-{\bf\Phi}_{\bf y}^{[j]}\|_{F}^{2}$, where $\|{\bf H}\|_{F}={\mbox{tr}\,({\bf H}^{\prime}{\bf H})}^{1/2}$ denotes the Frobenius norm of matrix ${\bf H}$. To evaluate the connectivity regime change-point detection, we measure the percentage of correctly classified time points into the true states for each simulated time course.

Results: Figure 4.1 plots the averages and standard deviations of the state classification accuracies for different estimators over all replications, as a function of dimension $N$. Both the factor-SVAR model-based estimates $\widehat{S}^{\text{SKF}}_{t}$ and $\widehat{S}^{\text{SKS}}_{t}$ perform better in regime segmentation than the K-means clustering, with substantially higher accuracy consistently for all $N$, albeit with higher standard deviations. The refined smoothed estimates $\widehat{S}^{\text{SKS}}_{t}$ based on the entire observations are more accurate than the filtered estimates $\widehat{S}^{\text{SKF}}_{t}$. Moreover, it can be seen that the accuracy of K-means clustering drops as $N$ increases, while for both the switching Kalman estimates, it tends to stabilize for high dimensions when $N\geq 30$. This may be because the regime partitioning was done based on the noisy estimates of high-dimensional TV-VAR coefficients fitted on short-windowed samples, compared to the lower-dimensional, reliably estimated subspace of factors in our approach. Another reason is the inherent limitation of the K-means algorithm itself neglecting temporal evolution of the connectivity states, which instead can be captured by the Markov chain of the switching model.

Figure 4.2 and Figure 4.3 plots the estimation errors of the directed connectivity matrix for the two states by the K-means clustering-based and the factor-SVAR model-based procedures, for increasing network dimensions $N$. The results are averages and standard deviations over the 100 replications, which respectively indicate the accuracy (unbiasedness) and consistency of the estimator. It is shown that the f-SVAR subspace estimators clearly outperform the $L_{2}$-regularized VAR estimator based on K-means clustered regimes, for both states and particularly for large $N$, in terms of significantly lower estimation mean squared errors and standard errors, and only slightly underperformed when $N$ is small. We can also see a rapidly growing trend of estimation errors in the K-means-based $L_{2}$ estimator as $N$ increases and approaches the regime sample size. In contrast, the robustness of the proposed f-SVAR estimators in high-dimensional settings is evident from the slower error rates (Figure 4.2-4.2(a)) and the constancy of standard errors over the increased dimensions (Figure 4.2-4.2(b)). These results can be explained by the more accurate regime segmentation by the SKS conditioned on the EM-estimated parameters as shown in 4.1, and improved consistency of the factor-based estimator over the ridge estimator for high-dimensional VAR coefficient matrix in each regime. The asymptotic theory of our proposed estimator such as convergence rates will be further studied in future work. Among the f-SVAR methods, both the coupled (common ${\bf Q}$) and decoupled (regime-dependent ${\bf Q}$) estimators perform comparably, despite slight superiority of the later. This suggests that the difference in directed connectivity structure based on a block-diagonal VAR model is mostly explained by inter-dependence in the factors, and less so in the projection of the underlying subspace. Hence, it can be sufficiently approximated by regime-dependent factor process, with a constant factor loading matrix across regimes.

5 Application to Estimating Dynamic Brain Connectivity

In this section, we shall apply the proposed f-SVAR approach to estimating time-evolving effective connectivity in high-dimensional resting-state fMRI data, characterized by abrupt transition of underlying quasi-stable brain states.

5.1 Resting-state fMRI Data

*1) Data acquisition:*We studied the resting-state fMRI data of 10 subjects from the first scan of a dataset publicly available at NITRC (http://www.nitrc.org/projects/trt). A Siemens Allegra 3.0-Tesla scanner was used to obtain three resting-state scans for each subject. During scans, the subjects were asked to relax and keep their eyes open. BOLD functional images were acquired using a T2-weighted gradient-echo planar imaging (EPI) sequence (TR = 2000 ms; time echo (TE) = 25 ms; flip angle (FA) = $90^{\circ}$; field of view (FOV) = 192 mm; voxel size = $3\times 3\times 3$ mm3; matrix $64\times 64$; number of slices = 39). A time-series of $T=197$ EPI volumes was collected for each scan.

2) Preprocessing: The data were preprocessed using the AFNI and FSL software packages as in [Fiecas et al., 2013]. The steps included (1) Motion correction using six-parameter rigid body transformation, normalized correlation as cost-function and referencing to the middle volume; (2) Spatial normalization to the Montreal Neurological Institute (MNI) template; (3) Probabilistic segmentation of the brain to obtain white matter and cerebrospinal fluid (CSF) probabilistic maps, thresholded at 0.99. (4) Removal of the nuisance signals, namely the six motion parameters, white matter and CSF signals, and the global signal. (5) Spatial smoothing with a 6 mm full-width half-maximum (FWHM) Gaussian kernel.

3) Parcellation: We used the automated anatomical labeling (AAL) atlas to obtain an anatomical parcellation of the whole brain into 90 ROIs with 45 regions in each hemisphere . In this study, the ROIs were grouped into six pre-defined resting-state system networks (RSNs) of similar anatomical and functional properties, based on the templates in [Allen et al., 2012, Li et al., 2011]. The considered RSNs include sub-cortical (SCN), AN: auditory (AN), sensorimotor (SMN), visual (VN), attentional (ATN) and default mode network (DMN). We followed the ROI abbreviations in [Salvador et al., 2005].

5.2 Results

We analyzed the dynamic states of large-scale effective brain connectivity in the resting state. We fitted a three-state factor-SVAR(1) model using the EM algorithm to the resting-state fMRI time series concatenated for all subjects, to identify the state transitions and the high-dimensional directed dependencies within each state which are assumed to be shared across subjects (as measured respectively by the SKS-estimated state-time sequence and state-dependent VAR coefficient matrices). Here, the decoupled SVAR subspace estimator was used, and the number of factors selected for this data by using BIC was $\hat{r}=14$. We used the VAR model order of one, as typically assumed for fMRI data [Valdés-Sosa et al., 2005].

Figure 5.1 shows the estimated whole-brain directed connectivity matrices between ROIs for three distinct states, and the corresponding within-network connectivity graphs for three selected RSNs. Only significant connections are shown, tested based on the asymptotic normality of the factor-VAR coefficient estimator in (25), at $\alpha=0.05$ with Bonferroni correction. Our method identifies the modular organization of the resting-state networks over all states, where ROIs within a functionally relevant network tend to be densely connected, but sparsely connected between different networks, particularly pronounced in VN, DMN and SMN. This characteristic has been reported in previous studies of static fMRI functional connectivity, e.g. [Ferrarini et al., 2009]. In consistency with findings in dynamic functional connectivity states [Allen et al., 2012, Hutchison et al., 2013a], our results also show the distinct large-scale connectivity patterns across different brain states in terms of variability in both the strength and sign of the connectivity and the network modularity. More interestingly, our method further reveals state-related difference in the directionality of the connections not reported previously, as evident from the asymmetry of the estimated VAR coefficient matrices. We discuss few apparent patterns that differ between the effective connectivity states. It is shown that the states are differentiated by the ROI-wise connectivity for both between-networks and within-networks. For the within-network connectivity, we observe the strongest connections between ROIs in state 1 for all the three RSNs, generally. For the sensorimotor networks, the directed interactions between central regions in the primary motor cortex is the strongest in state 1, which however shows disrupted connections with the parietal regions. In state 2, we found strong uni-directional influences from both the superior parietal nodes (SPG.L and SPG.R) to the supplementary motor area (SMA) with negative correlation (as indicated by blue edges), which are not present in states 1 and 3. For the attentional networks, we identified the lateral frontal-parietal network (similar to ventral attention network [Vincent et al., 2008]) between regions e.g. middle frontal gyrus and inferior parietal gyrus in all states, with the strongest connections occurred in state 1. However, we found denser directed information flows across both left and right hemispheres in states 2 and 3, compared to state 1. Particularly, it is interesting that the cross-hemisphere connections between parietal regions detected in states 2 and 3 were completely absent in state 1. For the default mode networks, state 1 also reveal the strongest and densest connections between ROIs related to posterior cingulate cortex (PCC)/precuneus, medial prefrontal cortex and the left and right inferior parietal lobule, with PCC correctly identified as a major hub of the DMN, strongly connected with other regions, as reported in numerous studies [Fransson and Marrelec, 2008].

To examine the transitions of the connectivity states in Figure 5.1 as a function of time, the estimated state-time alignment for the 10 subjects is shown in 5.2. The results suggest that the effective connectivity states changes over time and the pattern of changes varied across subjects. However, the connectivity states reoccur over time and shared across subjects. It also exhibits slow dynamics, where the connectivity tends to be assigned to single discrete states for long periods, with occasional fast switching between states. Moreover, the degree of the non-stationarity differs between subjects, from the rapid transitions between states (subjects 3, 5 and 6) to almost time-constant connectivity remained in particular states, i.e state 3 (subjects 2, 4, 8), state 2 (subject 2) and state 1 (subject 9). Note that state 1 (yellow) with enhanced connectivity for all RSNs exhibits the lowest occurrence in the time-courses over all subjects. Figure 5.3 shows the estimates for subject 6. The connectivity regime changes in the observed fMRI signals (Figure 5.3(a)) can be reflected in the lower-dimensional factor time series (Figure 5.3(b)). Besides, the SKS refines the state estimates by SKF, smoothing the spurious spikes and producing more stable regimes, as shown in (Figure 5.3(c)).

6 Conclusion

We developed a novel approach to identifying dynamic effective connectivity states with a large number of brain regions from fMRI data, based on a regime-switching factor model. The proposed approach first characterizes the high-dimensional fMRI data via a small number of factors for dimension reduction using a factor model in the observation space, and then performs connectivity regime segmentation in this low-dimensional latent factor subspace. By specifying the factor dynamics to follow a Markov-switching VAR process with a state-space formulation, it enables a reliable and efficient detection of change-points of the connectivity states using the Kalman smoothing and EM algorithm, and estimation of high-dimensional connectivity matrix for each state by projection of the estimated subspace parameters. The use of a regime-switching VAR specification allows us to examine state-driven changes in another important feature of connectivity, i.e. the directionality of connections, which are not addressed in earlier studies of dynamic functional (un-directional) connectivity. Hence, our approach provides a unified parametric framework for estimating both the time-varying connectivity structure and its quasi-stable state partitions, as distinct to using the separate steps of sliding-window connectivity analysis followed by K-means clustering. Moreover, the shortcomings of K-means clustering producing spurious fluctuations of states due to fixed-time windowing of time-varying connectivities and its failure to account for the temporal structure, can be overcome by the modeling with Markov chain which can capture both stable periods and the abrupt alternations of states via the transition probabilities.

Simulation results demonstrate the superiority of our approach over the K-means clustering of TV-VAR coefficients, giving more accurate estimation of the dynamic states, and the within-state connectivity graph, particularly in the high-dimensional settings. In analyzing the resting-state fMRI data, the proposed estimator confirms previous findings of non-stationary, re-occurring brain states during rest, and state-dependent modulation of large-scale connectivity patterns and modular structure. Furthermore, we produced new evidence for across-state difference in both the strength and directionality of directed information flows within resting-state networks. Future works will investigate different variants of the proposed framework, e.g. by allowing regime-switching in the factor loadings, instead of the factor dynamics itself. The method can also be extended to analyze time-varying directed coherence which measures connectivity at specific frequency of brain activity.