Deep Markov Factor Analysis: Towards Concurrent Temporal and Spatial Analysis of fMRI Data
Authors: Amirreza Farnoosh
Sarah Ostadabbas
Venue: NeurIPS 2021 Poster (2021)
Keywords: Deep generative models
Bayesian dynamical matrix factorization
Variational inference
Cognitive neuroscience
Functional MRI (fMRI)

Abstract: Factor analysis methods have been widely used in neuroimaging to transfer high dimensional imaging data into low dimensional, ideally interpretable representations. However, most of these methods overlook the highly nonlinear and complex temporal dynamics of neural processes when factorizing their imaging data. In this paper, we present deep Markov factor analysis (DMFA), a generative model that employs Markov property in a chain of low dimensional temporal embeddings together with spatial inductive assumptions, all related through neural networks, to capture temporal dynamics in functional magnetic resonance imaging (fMRI) data, and tackle their high spatial dimensionality, respectively. Augmented with a discrete latent, DMFA is able to cluster fMRI data in its low dimensional temporal embedding with regard to subject and cognitive state variability, therefore, enables validation of a variety of fMRI-driven neuroscientific hypotheses. Experimental results on both synthetic and real fMRI data demonstrate the capacity of DMFA in revealing interpretable clusters and capturing nonlinear temporal dependencies in these high dimensional imaging data.

Deep Markov Factor Analysis: Towards Concurrent Temporal and Spatial Analysis of fMRI Data
Amirreza Farnoosh
Augmented Cognition Lab (ACLab)
Department of Electrical & Computer Engineering
Northeastern University
Boston, MA 02115
farnoosh. www.ooostee.ott.teteeeeeeeee.ee edu
Sarah Ostadabbas
Augmented Cognition Lab (ACLab)
Department of Electrical Computer Engineering
Northeastern University
Boston, MA 02115 ostadabbas@ece. neu. edu
Abstract
Factor analysis methods have been widely used in neuroimaging to transfer high dimensional imaging data into low dimensional, ideally interpretable representa- tions. However, most of these methods overlook the highly nonlinear and complex temporal dynamics of neural processes when factorizing their imaging data. this paper, we present deep Markov factor analysis (DMFA), generative model chain of low dimensional temporal embeddings
that employs Markov property together with spatial inductive assumptions, all related through neural networks, to capture temporal dynamics in functional magnetic resonance imaging (IMRI) data, and tackle their high spatial dimensionality, respectively. Augmented with a discrete latent, DMFA is able to cluster IMRI data in its low dimensional temporal embedding with regard subject and cognitive state variability, therefore, enables validation of a variety of fMRI-driven neuroscientific hypotheses. Experimental results on both synthetic and real fMRI data demonstrate the capacity of DMFA in revealing interpretable clusters and capturing nonlinear temporal dependencies these high dimensional imaging data.
Introduction
Functional magnetic resonance imaging (fMRI) has been extensively used in cognitive neuroscience to study brain structures and their interactions at rest or during a cognitive task [Glover] 2011]. fMRI also provides insights into how brain's functional connectivity changes during different experimental conditions [Preti et al. 2017 Azari et al. 2020). However, due to high dimensional nature of fMRI data (tens of thousands of voxels in few seconds-long sessions), analyzing functional connectivity of brain could get very challenging interest (ROI)-based approaches reduce size of data for processing [Poldrack] 2007 [tzel et al.) 2009 Farnoosh and Soltanian-Zadeh) 2017] Gadgil et al. 2020]. However, averaging across many voxels within each region could wash out signals from small number of task-relevant voxels with noise from contaicivteed voxels, therefore results in a loss of information [Podrack 2007)
2013). Most classical methods employ region of
35th Conference on Neural Information Processing Systems (NcurIPS 2021)

---

Recently, few approaches have been proposed, based on probabilistic generative models, for topo- graphic factorization of IMRI data into a weighted summation of few localized activation sources (i.e., temporal weights and topographic spatial factors), among which topographic factor analysis (TFA) MManning et al., 2014b) and its multi-subject extension, hierarchical TFA (HTFA) Maaning et al. 2018], and neural TTA (NTFA) (Sennesh et al., 2020) are the most noted ones. This factor- ization serves as a necessary preparation step for subsequent statistical analysis that can effectively characterize subject- and stimulus-level variations and reveal task- or cognitive state-related networks in brain. However, TFA approaches assume a prior in which temporal weights are conditionally independent function of time, which means they do not encode temporal dynamics. Given the non-linearity and complex dimeneepenssscies inherent in fMRI, a model is required that can capture and represent these dependencies.
In this paper, we propose deep Markov factor analysis (DMFA)] a Bayesian model for factorization of fMRI data that learns a deep generative Markovian prior to reason about nonlinear temporal dynamics. This is realized by a chain of low dimensional temporal embeddings related through neural networks. This prior is further augmented by a discrete latent for multimodal dynamical estimation, and clustering subject- and task-level variations directly in its low dimensional temporal embedding. To accommodating high spatial dimensionality, DMFA generatively parameterize spatial factors from a low dimensional spatial latent through neural networks.
We evaluate the performance of DMFA on synthetic and two real large-scale fMRI datasets. Our experiments demonstrate that DMFA uncovers meaningful clusters in these data and achieves better predictive performance for unseen data relative to the state-of-the-art
2
Related Work
Factor Analysis in fMRI: Factor analysis in neuroimaging includes a wide range of approaches for reducing data dimensionality to facilitate their interpretability and computational tractability. Principal component analysis (PCA) [Pearson 1901] and independent component analysis (ICA) [Comon et al. 1991] are among the most well-known classical factor analysis methods. To accommodate tensor data and mitigate scalability issues, multilinear versions of PCA and ICA have been proposed
in Vasilescu and Terzopoulos [2005], Beckmann and Smith (2005). Cichocki] [2011), Richard and Montanari [0114). Hopkins et al. [2015]. Specifically, for multi- suujies fMRI study, Lee et al. [2008] proposed independent vector analysis (IVA) and Richard et al. [2020) developed Viewi ICA to model shared responses. Likewise, Chen et al. 2015] developed a shared response model (SRM) for aggregating multi-subject fMRI data and highlighting group differences, Van Kesteren and Kievit (202T) incorporated structured residuals into the exploratory factor analysis (EFA) framework. Karahanoglu and Van De Ville (2015) deconvolved hemodynamic response from rest fMRI time series and then performed temporal clustering on the resulting whole-brain innovation signals to recover the corresponding spatial patterns. However, spatial factors obtained by these methods unstructured, often have the same size as the images in the original dataset, and may include many small and large voxel clusters across the brain, therefore are not directly interpretable [Manning et al.) 2014b). More important, these methods are permutation invariant along temporal dimension (i.c., do not assume any relationships between temporal dimensions), therefore do not model temporal dynamics Yu et al. 2016).
To enhance spatial interpretability, Manning al. 2014b] introduced topographic factor analysis (TFA), a probabilistic technique that casts each brain image as a weighted sum of Gaussian spatial factors. In this way, each Gaussian blob (a.k.a topographic factor) can be thought of a brain region for which activation levels are estimated over time, and is easily interpreted through its set of parameters. Later on, authors proposed hierarchical version of TFA (HTFA) [Manning et al.) 2018) to incorporate data from multiple subjects and enable hypothesis testing across subjects by applying hierarchical Gaussian priors for global and subject-specific parameters of temporal and topographic factors. These linearly-depende of unimodal distributions and limit the model's 's expressivity. Recently, [Sennesh et al. 2020 proposed neural TFA (NTFA) for task fMRI, which extends TFA by incorporating neural networks onto its framework. NTFA assumes separate latent embeddings for participants and stimuli and map them into the temporal and spatial latents with neural networks. However, methods in this line of work
hicrarchical Gaussian priors, though, favor estimation
TThe source code is available at btttss//ggiubbccommossaaaaaa/btp

---

essentially assume prior in which temporal weights are conditionally independent as function of time, which means that they do not encode temporal dynamics.
Dynamical Factorization: A number of matrix/tensor factorization approaches have been proposed for modeling temporal dynamics in sequential data. These methods employ GauesiGn state-
space models [Sun et al. 2014]. autoregressive temporal regularizer [Bahadori et al.
2014 Yu et al.
2016 Takeuchi et al. 2017]. and multilinear dynamical systems [Rogers et al.
2013 Cai al.
2015 Jing et al. 2018]. In contrast to these methods which provide point estimates,
Bayesian
dynamical factorization has been proposed in Xiong et al. (2010) Charlin et al. [2015), Sun and
Chen [2019), in which linear temporal dynamics are applied to factor latents. While some of these methods have been successful in dynamical modeling of sequential data, they are less effective for
high dimensional spatial data like fMRI as they do not explicitly adopt any structural constraints for spatial factors. Moreover, their linear dynamical assumptions lack the capacity to characterize nonlinear dependencies.
Motivated by recent advances in deep learning, several studies have implemented neural networks
into Gaussian state space models for nonlinear dynamical modelling Krishnan et al., 2015
Watter
al. 2015
Chung al. 2015 Karl et al. 2017 Krishnan et al. 2017 Fraccaro et al. 2017 Becker
ct al.
2019 Farnoosh et al. 2021]. A common practice is to learn temporal latent model followed
by mapping to data space, a.k.a. decoder, and an encoder for amortized variational inference.
However, this encoding/decoding framework is not tractable in very high dimensional spatial data like fMRI. To be more specific, is computationally intensive feed high dimensional data for amortized estimation, discussion). Also, several studies have employed recurrent neural networks (RNNs) for temporal analysis of fMRT data (Hjelm et al.] 2018 Dvornek et al.) [018 Yan et al.] [019 Wang et al., 2019]. RNN-based methods are not probabilistic and can only process reduced dimensional ROI data. Although convolutional RNN models can capture both temporal and spatial correlations, by construction they do not provide the familiar spatial correlation maps that we expect in neuroimaging analysis and that are necessary for neuroscientific interpretabilityy In addition, RNN model are harder to train on long sequences due to the vanishing/explodig gradient problem as they feed parameters sequentially. While this problem is much alleviated in LSTM. GRU, and Transformer models, these architectures still result in a huge computational graph and intense GPU memory consumption on long sequences, which leads to very long training times. In contrast, DMFA, will be explained below, samples from a distinct posterior distribution at each time point and can fit arbitrarily long sequences without any issues (by allowing parallel computations over time), and at the same time infers smooth spatial factors automatically.
map to them directly from latent space (see Section 4.1 for detailed
3 Deep Markov Factor Analysis (DMFA)
Let's assume corpus of N IMRI data {Y, n E [1:N]] for Y, E TV, where T and V are the number of time points and voxels, respectively. DMFA defines hierarchical dynamical deep generative model over this corpus in factorization framework:
c
W
W2
WT
Yn ~ Norm(W, Fn, T) W, Norm (1² (2,") ow (Z"))
F, [ RBF Pk: 91)
P.7 Norm(16(=7), (z))),
Figure 1: Graphical model representation for deep Markov factor analysis (DMFA). DMFA incorporates a deep generative Markovian prior po po(2 ( to represent temporal variations in W, which is conditioned on a discrete latent, C, for data clustering. The K spatial factors conditioned
2 Norm(0, I),
2 ~ po 12. CC,).
C, Cat(T)).
where F, E RKVV are K < V spatial factors and Wn E KXT are their associated temporal weights, consistent with factorization framework. po(Z. C.) a deep generative Markovian prior over a set low dimensional temporal latents Z {> a it E [0:T]]
shared latent, z". Latent nodes and observations are represented by solid and gray-shaded circles, respectively. The solid black squares denote nonlinear mappings pa- rameterized by neural networks.
3

---

conditioned on a discrete latent Cv for clustering. The temporal weights W,, are sampled from Gaussian distribution whose mean and covariance are a function of 2 and are parameterized by neural networks 14 and og respectively. The spatial factors F2 are determined by set of K radial basis functions (RBFs) whose parameters {p. 7 are sampled from a Gaussian distribution parameterized by neural networks 1i and og from a low dimensional spatial latent ef All networks have parameters, which are collectively denoted by 0. Finally, ox denotes observation noise. The graphical model for DMFA is depicted in Fig. |
We assume that each fMRI data Y, belongs to a specific cluster s [1:S] which is declared by its associated discrete latent Cn ~ Cat(TT) sampled from a categorical prior with assignment probability vector и. This discrete latent conditions the first temporal latent 2n,0 resulting in Gaussian mixture distribution over this temporal latent:
Clustering Latent:
po(=w,olC., 8) Norm(us, E.);
where μla and Eg are determined by cluster assignment.
Deep Generative Markovian prior:
We assume that the temporal latents Z are related through Markov chain, for which the transition probability pe(= ".,,,l=,,,-)) = Norm (1z( (= oz (2W,e_))) a Gaussian distribution whose mean and covariance are parameterized by neural networks from
2n.t fn. We further blend the estimated mean from neural network with linear transformation of
w
n to support both nonlinear and linear transitions (we dropped n for brevity):
( (1 g) 3 M4 (=x-) + + @ plg N. (z("))
where (-) represents a linear mapping, fig (1, (.) is the nonlinear mapping of neural network, and 9 E [0, 1] a weighting vector estimated from = with another neural network.
Temporal Weights & Spatial Factors: As with the transition model, we assumed Gaussian distri- butions for temporal weights, such that their means and covariances are function of their associated temporal latents:
Wn,t Norm (=w,)), (=w,)))
where M6 and ow are neural network functions. The high spatial dimensionality in fMRI data encourages the need for a hierarchical analysis that summarizes spatial factors with fewer parameters. Therefore, consistent with Manning al. [2014bla 2018). we represent each spatial factor as radial basis function (i.e., a Gaussian blob):
fku == RBF (u; 7k) = exp (- Ilp, -0--PA px.l?
exp (1x)
where Pe R23 denotes the position of voxel with index v, and Pk: Th R denote the center and extent of the Gaussian blob, respectively. The set of all RBF parameters {pk, 7k k E [I:K]] ~ Norm (u (= (=x)) are sampled from a Gaussian distribution whose mean and covariance are function of the shared spatial latent z,. Introducing " as a low dimensional spatial embedding
encourages estimation of multi-modal distribution among spatial factors.
3.1 Learning DMFA with Variational Inference
We train DMFA with stochastic variational methods (Hoffman et al. 2013 Ranganath et al. 2013 These methods approximate the posterior po (W, Z 2", C, >\Y) using variational distribution (W, z" 2F, C, 7), where á denotes parameters of the variational model, by maximizing a lower bound on the log- likeliihldood of the data:
[(0, Φ) log logpo (Y) (L (W, 22, ll 2 2
C.p.)ly)))
By maximizing this bound with respect to the parameters 0 and we learn a deep generative model and perform Bayesian inference, respectively.
Parameterizing Variational Distribution: We assume fully factorized variational distribution and introduce trainable variational parameters for each data as mean and covariance of diagonal Gaussian:
K (W, (W, 2" "", ,,,,,,,C,, C, = II ((Cw)a(=W,)))(t) II APa,k )ink) II al ((>,,))(ww..)))
km

---

Table I: Network architectures in DMFA. The fully connected (FC) layers parameterize the three Gaussian distributions in the generative model: po(2 1=) po(w))=()) and po (p. 2/2"). second row denotes inputs to neural networks.
Model
14. 0o R R2
H4, o8 : IR22 TRKK
1 o R - R
Layer
Input
2t-1 IR2
2/" E R²
2" 122 FC PReLU 4C4 PReLU
FC2 De PReLU
FC2 D. PReLU
FC D, 2 Sigmoid
FC D. x 2D PReLU
3
9 [0, 1]²
FC 2D, x 2K
FC x 6K
4
FC x D, PReLU
k² log o" E
jf E R
5
FC D: x 2 p. E R
FC x PReLU
FC x PReLU
7
PReLU
FC x 6K
logoz o² R²
log of goK
9
FC2 x 2 2 IR?
where q(-) for continuous latents are parametric Gaussian distributions. Note that due the high dimensionality of fMRI data (V >> N, i.e., number of voxels are much larger than number of data), we do not utilize an amortized inference (where a single parametric function maps each data to set of variational parameters) as it could be computationally very expensive and result in over-fitting (see Section 1.11. Although we can define variational parameters for categorical distributions ((C.). we approximate it with posterior p(C,, [o,l== to compensate information loss induced by the mean-field approximation:
P(C, = =)p(>,,00f CC =
(C,, =) pCC, = s|2,,0)
Σ3 E. P(C,
s)
w IC,
This approximation has two-fold advantage: spares the model additional trainable parameters for variational distribution, and further links variational parameters (=(,=))) to generative parameters of po("...) and ps (c), hence results a more robust learning and inference algorithm.
Derivation of Evidence Lower Bound (ELBO): Once the generative model and variational distribu- tion are determined, we can derive the ELBO which contains a reconstruction term and regularization terms for each latent to bring their variational posterior as close as possible to their generative prior:
Cn (0, Φ)) Ey (Wa))a Jogpo(Y, Wn, Fw)
(Likelihood)
Σ (Cn)
(Clusters)
Σ Byg 1)
KL KLL(ss (1o0(==)))))))) (x,
(Markov Transitions)
Σ Byg (effee)
( (wn, (wo
(Weights)
- (+t)) KL Thokk Tpox)lp m =KL(a(<<)llb(z²))
(Factors)
We compute Monte Carlo estimate gradient of the ELBO using a reparameterized sample from the variational distribution of continuous latents. For the discrete latent, C,,, we compute the expectations over (C,) by summing over all the possibilities, hence no sampling performed. We analytically calculate the KL terms of ELBO for both multivariate Gaussian and categorical distributions, which leads to lower variance gradient estimates and faster training as compared to e.g., noisy Monte Carlo estimates often used in the literature.
5

---

Zg
Neurrsi ActiVity
group?
Real
BOLD Petpo non: neural rectuity convakved with HAT
0.5-
Recovared Wilght: in- ma-ale
Presiitee weigte ii co-plaggi piaio.
0.4-
trajectory
Reconstructes
0.0-
-0.2
100
time
(a)
(b)
(c)
Figure 2: (a) DMFA recovered the three clusters of activation in our synthetic IMRI dataset, unsupervised. The mean dynamical trajectory (composed of three consecutive rotational dynamics) shows the inferred trajectory of each cluster mean over time in the temporal latent and is consistent with the periodic activation of sources in data clusters. (b) Real and reconstructed brain images. (c) The learned generative model's 's predictions for selected activation source show that DMFA encoded the nonlinear hemodynamic response function its deep temporal generative model.
4 Training Details
We described the network architectures of DMFA in Table where D, and De are the dimensions of hidden layers for (Hg, og) and ""), respectively. We did all the programming in PyTorch v1.3 and used Adam optimizer with learning rate of 0.01. We initialized all the parameters randomly except for locations of Gaussian blobs for which set initial values to local extrema in their averaged fMRI data. We clipped Gaussian blob parameters to the confines brain if needed. We used a linear KI annealing schedule, [Bowman al., 2016), which increased from 0.01 to 1 over 100 epochs. We learned/tested all models on an Intel Core 17 CPU @3.7 GHz with Gigabytes of RAM, which proves tractability of the learning process. Per-epoch training time varied from 0.1 minutes in small datasets to 6.0 minutes in larger experiments and 200 epochs sufficed for most experiments.
4.1 Discussion Parameter Count for DMFA
The number of learnable parameters for variational distribution in DMFA is O(NTK). DMFA has O(KD++D,) parameters for its temporal generative model and O(K) parameters for its spatial generative model. Note that the clustering latent, C. does not impose additional parameters to the variational distribution, while only adds O(S) parameters to the temporal generative model.
While DMFA introduces extra features and more complex modeling assumptions for fMRI experi- ments compared to TFA methods Manning et al. [2014b]] a 2018), (i.e., infers nonlinear temporal dynamics and performs clustering), we emphasize that it has the same order of parameters these methods. TFA methods similarly have O(NTK) parameters as they employ fully factorized varia- tional distribution. We also want to highlight that DMFA is tractable in very high dimensional spatial data like IMRI in contrast to nonlinear state-space models of Watter et al. [2015), Krishnan et al. 2017 Fraccaro al. 2017]. Karl et al. [2077) Becker et al. 2019]. The encoder/decoder structure in these works, i.e., go W kn)/po Yn, W.). immediately scales both generative and variational parameters to least OKV), where V 10 NT in fMRT data, hence causes extensive computa- tional burden and more importantly overfitting. Furthermore, these methods do not learn a generative (((p7) 7), and as result are not able reason about subject-level variabilities in this respect. We overcome these challenges in DMFA by carefully designing our non-amortized variational inference and imposing functional form assumptions on spatial factors in factorization framework. The proposed learning and inference algorithm keep generative parameters O (KD, + D,) & O(KV), and variational parameters in ONNKK), where NT 102 -10S, yield an observation to parameter ratio of O NTV NIK = 0 k for all the experiments, therefore permit an
model for spatial factors, i.e., Po(p,
efficient learning process.

---

Jubbin Oeton
Real
Callach
MoxMMn
EBL
Stanford
Figure 3: (Left)) DMFA's clustering results show that the ASD and control groups can be partially separated. (Middle) Real and stronetd brain images, showing the smoothing given by sparse factorization. (Right) downstream classification task showed that DMFA and SVF outperformed regions of interest +SVM) and ICA-SVK in the Caltech, SBL, Stanford, and Yalc subsets of the dataset at p-value < 0.005, corrected for multiple comparisons. ROI+SVM and ICA+SVM performed better in the NYU subset.
5 Experimental Results
We analysed the performance of DMFA on a synthetic and two large-scale fMRI datasets. In the first experiment, we verified DMFA's capability in recovering the true clusters and capturing the underlying temporal dynamic in a synthesized fMRI dataset (visualized in Fig.2). Next, we evaluated the performance of DMFA on a large scale resting-state 1MRI data, Autism dataset [Craddock et al.) 2013], and a task IMRI data, Depression dataset [Lepping et al. 2016]. We assessed the clustering Teature of DMFA on both datasets in terms of disease and cognitive state separation tasks, visualized in Fig. and Fig. We further provided a quantitative comparison with two state-of-the-art Bayesian generative models for fMRI data, HTFA [Manning et al. 2018) and NTFA [Sennesh al., [220]. and deep state-space model, RKN [Becker et al. [2019]. in terms of synthesis quality of the generative models on both datasets by computing held-out likeliood and prediction accuracy, respectively. in Table 2
5.1 Synthetic fMRI
We generated synthetic fMRI data using MATLAB package provided by Manning et al. (014bb) which is known to be useful for analysing fMRI models. The synthesized brain image for each trial (time point) is weighted summation of number of radial basis functions (spatial factors) randomly located in brain. This synthesized fMRI data is then convolved with a hemodynamic response function (HRF) and zero-mean Gaussian noise with a medium-level signal-to-noise ratio is added. Here, we considered 30 activation sources (spatial factors) randomly located in a standard MNI- -3mm brain template with roughly V 270, 000 voxels and 150 trials. We randomly split these 30 activation sources into 3 groups, each having 10 of the Gaussian blobs. These three groups of sources are periodically activated according to some random weights, one after another, for trials. We generated non- non-overlapping sequences of T = 5 time points from this synthetic fMRI data. This resulted in 10 data points for each activation group (i.e., N = 30). To train DMFA, we set T = 5, K = 30, D. = 2, D. 8, S 3, and o = 0.01. As depicted in Fig. (a), DMFA was able to successfully recover the 3 clusters of activation that were present in this dataset totally unsupervised. The dynamical trajectory of each cluster mean in the temporal latent (i.e., May les; C) is predicted sequentially from the learned generative model and is visualized in the bottom-right of Fig. 2 (a). and appears to be partitioned into three consecutive rotational dynamics. This is consistent with the periodic activations of sources in data clusters (which come in tandem). Predictions of the learned generative model for a selected activation source are visualized in Fig. (d) for the next 50 time points, estimated follows: w; ~ p( p( (z),, where 2 ~ p(zz 21-1) for t = [151..... 200}. These predicted samples perfectly follow hemodynamic response function, confirming DMFA's capacity in capturing the underlying nonlinear HRF by using neural networks in its deep temporal generative model.
5.2 Autism Dataset
We used the publicly available preprocessed resting state fMRI (rs-fMRI) data from the Autism Brain Imaging Data Exchange (ABIDE) collected at 16 international imaging sites (Craddock et al., 2013]. This dataset includes rs-fMRI imaging from 408 individuals suffering from Autism Spectrum
7

---

JUCE
Real
Reconstructed
W
MDD+Control
MDD
Control
(a)
Figure 4: (a), Left: Training DMFA clustered together temporal latent variables associated with each subject without supervision, while partially separating clusters of points associated with the MDD group from those associated with the control group. The MDD group appears more concentrated into the center of the temporal latent space, while the control group have their temporal latent variables dispersed more broadly across the latent space, (a), Middle, Right: DMFA enabled us to partially separate "positive" and "negative" stimuli per-subject with Gaussian clusters. (b) Real and reconstructed brain images.
Disorder (ASD), and 476 typical controls. Each scan has T 145 ~ 315 time points at TR 2. and V = 271, 633 voxels. We split the signals into sequences of 75 time points. We take two approaches to evaluate the performance of DMFA in separating ASD from control: (i) Cluster data directly in the low dimensional latent, Z", using the clustering feature of DMFA (called DMFA), A), (ii) Extract functional connectivity matrices, (Hull et al. 2017). from learned temporal weights, W, followed by 10-fold SVM for classification (called SVM)+ As baselines, we performed 10-fold SVM classification on extracted connectivity matrices from (i) averaged signals of 116 ROIS in automatic anatomical labeling (AAL) atlas [Kazeminejad and Sotero 2019 (called R and (ii) hundred time courses obtained with nonlinear spatial ICA of Hyvarinen and Morioka 2017 (called ICA+SVM).
Several studies have been done on this dataset to differentiate ASD group from control [Abraham
et al. [017]] Parisot et al. 2017 Singhh et al. 2017 Kazeminejad and Sotero 2019 Wang et al. 2021 Sun al. [222]). all of which using supervised methods, and could achieve accuracies up 72% using the signals extracted from anatomically labeled regions in the brain by carefully splitting data to be as homogeneous as possible and redueing site-related variability.
We set T = 75, K = 100, D. = 15, D, = 5, S = 2, = 0.01, and trained DMFA for 200 epochs on the entire dataset (Full), and also datasets from 9 imaging sites (with more balanced datasets) separately: Caltech, Leuven, MaxMun, NYU, SBL, Stanford, UM, USM, Yale. As shown in Fig. 3 (Right), DMFA and DMFA+SVV outperformed baselines in Caltech. SBL, Stanford, and Yale (at p-value < 0.005, corrected for multiple comparisons), while ROI+SVM and ICA+SVM only performed better in NYU dataset. DMFA+SVM performed slightly better than baselines on the entire dataset (Please note that DMFA is clustering approach, hence, no error bars are provided for it in Fig. 33.. Clustering results for Caltech. Maxmun, SBL, and Stanford are shown in Fig. (Left) in which ASD and control seems to be partially separable (see Fig. S1 in supplementary for more visualization results).
5.3 Depression Dataset
In this dataset (Lepping et al. 2016). 19 individuals with major depressive disorder (MDD) and 20 never-depressed control participants listened to standardized positive and negative emotional musical and nonmusical stimuli during fMRI scanning. Each participant underwent 3 musical, and 2 nonmusical runs each for 105 time points at TR=3 with V 353, 600 voxels. During each run, each stimulus type (positive, and negative) was presented for 33 seconds (~ 11 time points) interleaved with instances of neutral tone of the same length. We discarded instances of neutral tone, and split each run into non-overlapping sequences of T = 6 time points in agreement with stimuli design (each stimuli block is split into two sequences). other words, each run has 4 sequences associated with "positive stimuli", and 4 with "negative stimuli" resulting in total of 8 data points for each run.

---

Table 2: Comparison of held-out log-likelihood and prediction accuracy. DMFA results in models with higher held-out likelihood and better prediction accuracy, therefore, is better fit compared to the competitors.
Negative Held-out Log Likelihood (Nats)
Prediction (NRMSE%)
Model
Ablation
DMFA
NTFA
HTFA
DMFA
RKN
Dataset
DMFA
DMFA
Autism (Caltech)
2.52
2.61
2.82
2.75
2.63
5.44
6.19
Depression
5.82
5.97
6.64
6.26
6.02
6.35
7.49
Synthetic
3.47
3.64
3.71
3.83
3.61
2.07
2.84
The best results are highlightod in bold fonts.
Likelthood values are scaled by 10 9.10 .000 and 10 for Autism, Depression, and Synthetic datasets, respectively.
In the first experiment, we trained DMFA the entire musical runs (N = 39 x 3 x 105-330 - by setting T = 6, K = 100, D. = 15, D, = 5, a = 0.001 for 200 epochs. The results are shown in Fig. 4 (a. Left). We observed that DMFA fully separated the data points associated with each subject into distinct clusters across the low dimensional temporal latent space. In other words, DMFA was able to re-unite pieces of signals associated with each subject without any supervision. More importantly, DMFA was able to partially separate the data points associated with MDD group from control. As seen in Fig. (a, Left), the MDD group data points are fairly populated in the center of temporal latent while control group are dispersed across the latent space. However, DMFA was not able to meaningfully separate "negative" and "positive" music pieces in its low-dimensional latent space from a subject-level perspective, since the variation between runs of subject dominates the stimulus-level variation. For this reason, in a second experiment, we focused on subjects, and their first musical run from both MDD and control group and trained DMFA respectively. Again, as expected, data points from each subject were distinctly clustered in the latent space (see middle and right columns in Fig. (a)). Additionally. DMFA was able to fit two partially separating Gaussians to "positive", and "negative" stimuli per subject. However, since the number of data points for each subject and run is limited, the significance of these clusters are not conclusive. A dataset with longer runs could possibly answer that.
5.4 Comparison with State-of-the-Art
We further evaluated DMFA against HTFA Manninng et al., 2018). an established probabilistic generative model for multi-subject fMRI analysis, which uses unimodal Gaussian priors for both temporal weights and spatial factor parameters, and its neural network-based extension NTFA [Sennesh al. 2020 in terms of held-out log-likelitee likelihood, and a state-of-the-art deep state-space model, recurrent Kalman networks (RKN) [Becker et al. [2019], in termss of next- nexet-time-eee time-point prediction accuracy of test set (see T Held-out likelihood is an established metric for evaluating generative models (in lack of a ground -tuthh and measures how probable an unseen data is under the generative model. For next-time-poin prediction, we adopt rolling prediction scheme as in Chen et al. [2019). ) and predict next time point on the test set sequentially from historical data using the generative model and spatial factors learned on the train set.
We split our data to train and test as follows. We used the Caltech site subset from autism dataset and split each subject's fMRI time series into two half (each with T 70). We trained the models on the first half, and tested on the second half. For the depression dataset, we considered 4 sequences from each subject's 's run for training and tested on the remaining sequences. We set K 100 spatial factors for all the models and trained each with their default hyperparameters. For the RKN model, we fit AAL atlas ROI signals as it is intractable to feed raw high-dimensional data. We reported average normalized prediction eтror (NRMSE%) ) for both RKN and DMFA on AAL atlas regions to make fair comparison. For the synthetic fMRI, we set K 30 and picked half of the data for test, and used the 30 Gaussian blob activations as ROI signals for reporting prediction accuracy in RKN and DMFA. The results are shown in Table 2 which proves that DMFA resulted in models with higher likelihood on the test sets, hence it is a better fit when compared to NTFA and HTFA models. DMFA also resulted in better predictions on the test sets of autism, depression, and synthetic datasets when compared with RKN.
9

---

5.5 Ablation Study
We performed two ablation studies further evaluate (and quantify) the impact of deep neural networks and Markov temporal apriori. First, we trained a version of DMFA in which we removed nonlinear activation functions (i.e., PReI.U), denoted as DMFA. As reported in Table 2 our datasets, this resulted in decrease in held-out likelihood. Also, the inferred temporal latents lacked the interpretable patterns seen in Fig. and Fig. F in terms of clinical/cognitive separability as latents collapsed into the conditional Gaussian priors. Second, we trained version of DMFA in which we removed the temporal connections (i.e., temporal transition model), denoted as DMFA. Again, this resulted in drop in held-out likelihood as reported in Table 2 however, the temporal latents preserved their interpretable patterns.
Conclusion
We presented deep Markov factor analysis, DMFA, new probabilistic model for robust factor analysis of high dimensional fMRI data. We employed chain of low dimensional Markovian latents connected by deep neural networks and conditioned on a discrete latent as stte-space embedding for temporal weights to account for nonlinear dynamics, enable data clustering in a low dimensional space, and provide informative visualizations about data. To tackle high spatial dimensionality in fMRI, we employed low dimensional spatial embedding to generatively parameterize spatial factors. DMFA proves fast and capable on large-scale fMRI data.
Acknowledgment
This work was supported by award numbers 1944964 and 1835309 from the national science foundation (NSF). The authors would also like to thank the anonymous reviewers for their insightful feedback.
References
Alexandre Abraham, Michael Milham, Adriana Martino, R Cameron Craddock, Dimitris Samaras. Bertrand Thirion, and Gael Varoquaux. Deriving reproducible biomarkers from multi-site resting-state data: An autism-based example. Neurolmage, 147:736-7455 745, 2017.
Bahar Azari, Christiana Westlin, Ajay B Satpute. J Benjamin Hutchinson, Philip A Kragel, Katie Hoemann, Zulqamnai Khan, Jolie B Wormwood, Karen S Quigley. Deniz Erdogmus, et al. Comparing supervised and unsupervised approaches emotion categorization in the human brain, body, and subjective experience. Scientific reports, 10(1):1-17, 2020.
Mohammad Taha Bahadori, Qi Rose Yu, and Yan Liu. Fast multivariate spatio-temporal analysis via low rank tensor learning. In Advances in neural information processing systems, pages 3491-3499, 2014.
Philipp Becker, Harit Pandya, Gregor Gebhardt, Cheng Zhao, C James Taylor, and Gerhard Neumann. Recurrent kalman networks: Factorized inference in high-dimensional deep feature spaces. In International Conference on Machine Learning, pages 544-552. PMLR, 2019.
Christian F Beckmann and Stephen M Smith. Tensorial extensions of independent component analysis for multisubject finri analysis. Neuroimage, 25(1):294-33 2005.
Samuel R Bowman, Luke Vilnis, Oriol Vinyals, Andrew M Dai, Rafal Jozefowicz, and Samy Bengio. Generating sentences from continuous space. In 20th SIGNLL Conference on Computational Natural Language Learning, CoNLL 2016, pages 10-21. Association for Computational Linguistics (ACL), 2016.
Yongjie Cai, Hanghang Tong, Wei Fan. Ping Ji, and Qing Facets: Fast comprehensive mining of coevolving high-order time series. In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 79-88, 2015.
Laurent Charlin, Rajesh Ranganath, James McInerney, and David M Blei. Dynamic poisson factorization. In Proceedings of the 9th ACM Conference Recommender Systems, pages 155-162, 2015.
Po-Hsuan Chen, Janice Chen, Yaara Yeshurun, Uri Hasson, James V Haxby, and Peter J Ramadge. A reduced- dimension firr shared response model. NIPS, volume 28, pages 460-468, 2015.
10

---

Xinyu Chen, Zhaocheng He. Yixian Chen, Yuhuan Lu, and Jiawei Wang. Missing traffic data imputation and pattern discovery with bayesian augmented tensor factorization model. Transporiation Rexearch Part C: Emerging Technologies. 104:66-77, 2019.
Junyoung Chung, Kyle Kastner, Laurent Dinh, Kratarth Goel, Aaron Courville, and Yoshua Bengio. recurrent latent variable model for sequential data, Advances in neural information processing systems, 28: 2980-2988, 2015.
Andrzej Cichocki. Tensor decompositions: new concepts in brain data analysis? Journal of the Society of Instrument and Control Engineers, (007)507-5166 2011.
Pierre Comon, Christian Jutten, and Jeanny Herault. Blind separation of sources, part ii: Problems statement. Signal processing, 24(1):1 1-20, 1991.
Cameron Craddock, Yassine Benhajali, Carlton Chu, Francois Chouinard, Alan Evans, András Jakab, Bud- hachandra Singh Khundrakpam, John David Lewis, Qingyang Li, Michael Milham, al. The neuro burcau preprocessing initiative: open sharing preprocessed neuroimaging data and derivatives. Frontiers in Neuroinformatiis, 7, 2013.
Nicha C Dvornek, Daniel Yang. Pamela Ventola, and James Duncan. Learning generalizable recurrent neural networks from small task-fmri datasets. In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 329-337. Springer, 2018.
Joset A Etzel, Valeria Gazzola, and Christian Keysers. An introduction to anatomical roi- fmri classification analysis. Brain research, 1282:1 14-125, 2009.
Amirreza Farnoosh and Hamid Soltanian-Zadeh. Information theoretic evaluation of brain connectivity compared cross correlation using simulated resting state fmri data. In 2017 Iranian Conference Electrical Engineering (ICEE), pages 107-112 JEEE, 2017.
Amirreza Farnonsh, Bahar Azari, and Sarah Ostadabbas. Deep switching auto-regressive factorization: Applica- tion to time series forecasting. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 7394-7403, 2021.
Marco Fraccaro, Simon Kamronn, Ulrich Paquet, and Ole Winther. A disentangled recognition and nonlinear dynamics model for unsupervised learning. In Advances in Neural Information Processing Systems, pages 2017.
Soham Gadgil, Qingyu Zhao, Adolf Pfefferbaum, Edith V Sullivan, Ehsan Adeli, and Kilian M Pohl. Spatio- temporal graph convolution for resting-state fmri analysis. International Conference Medical Image Computing and Computer-Assisted Intervention, pages 528-538. Springer, 2020.
Gary H Glover. Overview of functional magnetic resonance imaging. Neurosurgery Climics, 22(2):133-139, 2011.
Devon Hjelm, Eswar Damaraju, Kyunghyun Cho, Helmut Laufs, Sergey Plis, and Vince D Calhoun. Spatio-temporal dynamics of intrinsic networks in functional magnetic imaging data using recurrent neural networks. Frontiers in neuroscience, 12:600, 2018
Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley. Stochastic variational inference. The Journal of Machine Learning Research, 14(1):1003-144) 2013.
Samuel B Hopkins, Jonathan Shi, and David Steurer. Tensor principal component analysis via sum-of-square proofs. In Conference Learning Theory, pages 956-1006, 2015.
Jocelyn V Hull, Lisa B Dokovna, Zachary J Jacokes, Carinna M Torgerson, Andrei Irimia, and John Darrell Van Horn. Resting-state functional connectivity in autism spectrum disorders: review. Frontiers in psychiatry, 7:205, 2017.
Aapo Hyvarinen and Hiroshi Morioka. Nonlinear temporally dependent stationary sources. In Artificial Intelligence and Statistics, pages 460 469. PMLR, 2017.
Peiguang Jing, Yuting Su, Xiao Jin, and Chengqian Zhang. High-order temporal correlation model learning for time-series prediction. IEEE transactions on cybernetics, 2366))2385 2018.
Fikret 1şik Karahanoglu and Dimitri Van De Ville. Transient brain activity disentangles fmri resting-state dynamics in terms of spatially and temporally overlapping networks. Nature communications, 101)11100, 2015.

---

Maximilian Karl, Maximilian Soelch, Justin Bayer, and Patrick van der Smagt. Deep variational bayes filters: Unsupervised learning of state space models from raw data. In 5th International Conference on Learning Representations. ICLR 2017. Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net, 2017. BtLLhttss//pppeeeviieee net/
Amirali Kazeminejad and Roberto Sotero. Topological properties of resting-state fmni functional networks improve machine learning-hasee based autism classification. Frontiers in neuroscience, 12:1018, 2019.
Rahul Krishnan, Uri Shalit, and David Sontag, Structured inference networks for nonlinear state space models. In Proceedings of the AAAI Conference Artificial Intelligence, volume 31, 2017.
Rahul G. Krishnan, Uri Shalit, and David A. Sontag, Deep kalman filters. CoRR, abs/1511.05121, 2015, URL http://arxiv.ooe www.ii..rrr/aaa.. www.//g//gggr/.orae 05121
Jong-Hwan Lee, Te-Won Lee, Ferenc A Jolesz, and Seung-Schik Yoo. Independent vector analysis (iva):
multivariate approach for fmri group study. Neuroimage, 400 10:6-10,, 2008
Rebecca J Lepping, Ruth Ann Atchley, Evangelia Chrysikou, Laura E Martin, Alicia A Clair, Rick E Ingram, W Kyle Simmons, and Cary Savage. Neural processing of emotional musical and nonmusical stimuli in depression. PloS one, (66) 56859, 2016.
Jeremy R Manning, Rajesh Ranganath, Waitsang Keung, Nicholas Turk-Browne, Jonathan D Cohen, Kenneth Norman, and David M Blei. Hierarchical topographic factor analysis. In 2014 International Workshop on Pattern Recognition in Neuroimaging, pages 1-4. IEEE, 2014a.
Jeremy R Manning, Rajesh Ranganath, Kenneth A Norman, and David M Blei. Topographic factor analysis: bayesian model for inferring brain networks from neural data. PloS one, 9(5)e99914, 2014b.
Jeremy R Manning, Xia Zhu, Theodore L Willke, Rajesh Ranganath, Kimberly Stachenfeld, Uri Hasson, David Blei, and Kenneth A Norman. probabilistic approach to discovering dynamic full-brain functional connectivity patterns. Neurolmage, 2018.
Sarah Parisot, Sofia Ktena, Enzo Ferrante, Matthew Lee, Ricardo Guerrerro Moreno, Ben Glocker, and Daniel Rueckert. Spectral graph convolutions for population-based disease prediction. In International conference on medical image computing and computer assiised intervention, pages 177- Springer, 2017.
Karl Pearson. Liii, on lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(1 1):559-572, 1901.
Russell A Poldrack. Region of interest analysis for fmrr. Social cognitive and affective neuroscience, 2(1):67-70. 2007.
Maria Giulia Preti, Thomas AW Bolton, and Dimitri Van De Ville. The dynamic functional connectome: Statc-of-the-art and perspectives. Neuroimage, 160:41-54, 2017.
Rajesh Ranganath, Chong Wang. Blci David, and Eric Xing. An adaptive learning rate for stochastic variational inference. International Conference on Machine Learning, pages 298-306, 2013.
Emile Richard and Audrea Montanari. statistical model for tensor pea. Advances Neural Information Processing Systems, pages 2897-2905, 2014.
Hugo Richard, Luigi Gresele, Aapo Hyvärinen, Bertrand Thirion, Alexandre Gramfort, and Pierre Ablin. Modeling shared responses in neuroimaging studies through multiview ica. In Advances in Neural Information Processing (NeurIPS2020), 2020.
Mark Rogers, Lei Li, and Stuart J Russell. Multilinear dynamical systems for tensor time series. In Advances in Neural Information Processing Systems, pages 2634-2642, 2013.
Eli Sennesh, Zulgarnain Khan, Yiyu Wang. Benjamin Hutchinson, Ajay Satpute, Jennifer Dy, and Jan Willem van de Meent. Neural topographic factor analysis for fimr data. Advances in Neural Information Processing Systems, 33, 2020.
Chandan Singh, Beilun Wang, and Yanjun Qi. A constrained, weighted-11 minimization approach for joint discovery of heterogeneous neural connectivity graphs. CoRR, abs/ 1bb/1709900490, 2017. URL http://arxx org/ wwwpab///........ 04090
John Z Sun. Dhruv Parthasarathy, and Kush R Varshney. Collaborative kalman filtering for dynamic matrix factorization. IEEE Transactions Signal Processing, (2(14):399-3300 2014.
12

---

Lei Sun, Yanfang Xue, Yining Zhang, Lishan Qiao, Limei Zhang. and Mingxia Liu. Estimating sparse functional connectivity networks via hyperparameter-Inee free learning model. Artificial Intelligence in Medicine, 111: 102004, 2021.
Lijun Sun and Xinyu Chen. Bayesian temporal factorization for multidimensional time series prediction. CoRR, abs/1910.06366, 2019. http:tttt///xiiiorrr//// wtt.:///rggb///.....oo.. 06366
Koh Takeuchi, Hisashi Kashima, and Naonori Ueda. Autoregressive tensor factorization for spatio-temporal predictions. In 2017 IEEE International Conference on Data Mining (ICDM), pages 1105 1110. IEEE, 2017.
Nicholas B 2013.
Functional interactions big data in the human brain. Science, 342(6158):580-5884)
Erik-Jan Van Kesteren and Rogier A Kievit. Exploratory factor analysis with structured residuals for brain network data, Network Neuroscience, 5(1):1-27, 2021.
M. Alex O Vasilescu and Demetri Terzopoulos. Multilinear independent components analysis. In 2005 IEEE Computer Society Conference on Computer Vision and Puttern Recognition (CVPR '05), volume 1, pages 547-553. IEEE, 2005.
Lebo Wang. Kaiming Li, Xu Chen, and Xiaoping P Hu. Application convolutional recurrent neural network for individual recognition based on resting state fmri data. Frontiers neuroxcience, 13:434, 2019.
Mingliang Wang, Jiashuang Huang, Mingxia Liu, and Daoqiang Zhang. Modeling dynamic characteristics of brain functional connectivity networks using resting-state functional mri. Medical image analysis, 71:10200 71:102063, 2021.
Manuel Watter, Jost Springenberg, Joschka Boedecker, and Martin Riedmiller. Embed to control: locally linear latent dynamics model for control from raw images. In Advances in neural information processing 2015.
systems, pages
Liang Xiong. Xi Chen, Tzu-Kur Kuo Huang, Jeff Schneider, and Jaime G Carbonell. Temporal collaborative filtering with bayesian probabilistic tensor factorization. In Proceedings the 2010 SIAM international conference on data mining, pages 211 211-221 1-222. SIAM, 2010.
Weizheng Yan, Vince Calhoun, Ming Song. Yue Cui, Hao Yan, Shengfeng Liu, Lingzhong Fan, Nianming Zuo, Discriminating schizophrenia using recurrent neural network applied on time
Zhengyi Yang. Kaibin Xu, courses of multi--ite fmri data. EBioMedicine, 47:543-555 2019.
Hsiang-Fu Yu, Nikhil Rao, and Inderjit Dhillon. Temporal regularized matrix factorization for high-dimensional time series prediction. In Advances neural information processing systems, pages 847-855, 2016.
13


Reviews:
[{"review_id": "zaf1UNBklnD", "reviewer": "Reviewer_tynp", "summary": "The paper presents a deep Markov factor analysis (DMFA), which is a generative model that employs Markov property in a chain of low dimensional temporal embeddings as well as spatial inductive assumptions in order to capture temporal dynamics in functional magnetic resonance imaging (fMRI). The paper shows that DMFA has capability to cluster fMRI data in its low temporal embedding in regards to subject and cognitive state variability. The advantages of DMFA are demonstrated both through synthetic and application data.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper meets and exceeds the criteria of originality, quality, clarity and significance. The paper clearly indicates previous methodologies in the direction of their method (i.e., factor analysis and dynamic factorization), while providing ample citations for each. The development of DMFA is clearly presented with explicit and easy-to-understand diagrams. The estimation method (variational inference) is also clearly presented and makes it easy to reproduce for any interested practitioner. The simulation experiments show the clear utility of the  method for the generative model to learn from complex spatiotemporal dynamics presented in fMRI data. And finally, several application datasets are presented (e.g., in Autism, Depression, Ablation, etc.). The reviewer sees no major issues with the current work, and sees the paper as well-polish and an excellent example of a high-quality submission.\n\nBut due to some comments from other reviewers, I do see some significant issues with methodology/innovation, as well as a lack of proper sensitivity analysis. Therefore I will reduce my score.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "0.5"}, "extra_text": {}}, {"review_id": "h2FSH79yd6L", "reviewer": "Reviewer_FHJp", "summary": "This paper introduces Deep Markov Factor Analysis (DMFA). DMFA can model spatiotemporal BOLD fMRI time series by using time-varying latent factors that govern time-varying spatial maps consisting of spatial mixtures of RBFs. Latent time series are associated to a categorical clustering variable in order to enable grouping. Certain links in the Bayesian diagram are modeled using deep neural networks. Optimization is done by a variational approach.\nExperiments are performed on a simulated data set with known ground truth, the ABIDE autism data set, and a Depression study data set. Some clustering of relevant clinical outcomes is observed.\n\n\n", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "DMFA lies at the intersection of the Topographical Factor Analysis approach (ref 14), unsupervised TV-regularized resting-state analysis as performed in (https://www.nature.com/articles/ncomms8751) and variational autoencoders for time series (https://arxiv.org/abs/1506.02216). The latter two papers are not mentioned.\n\nSeveral concerns arise around this work:\n\nWhat is the point in using Gaussian blobs to model spatial maps? The argument that the idea is inherited from TFA is insufficient. Why not just use voxel maps at a certain resolution? Not all brain regions are shaped like spherical Gaussians or sparse sums thereof. The reconstructed maps show this very clearly -- pretty much all detail is lost.\n\nWhat is the point of *temporally varying* spatial Gaussian blobs? Surely brain regions remain in place over time. One could make the argument that “activity moves around”, but it remains true that certain spatially fixed regions receive more or less blood over time. It seems that this feature adds unnecessary complexity and hinders interpretability\n\nThe synthetic model is not very helpful since it basically generates data according to the model used. Why not generate some fake brain data based on a different model, like ICA components?\n\nThe experimental evaluation does not seem to contain any clearly stated numbers enabling comparisons to other methods. It is great that the method finds some factors in an unsupervised manner, but please state the accuracy over the full datasets for the task. It may be helpful to look at Abraham et al 2017 (ref from the paper) for the data curation procedure as this is an important thing to get right.\n\nWhy is there no comparison to the method presented in (https://www.nature.com/articles/ncomms8751), adapted to this task? This provides a simpler baseline.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "There is no reason to believe there would be a negative societal impact to this work. That there would be a strongly positive one is also unlikely.\n\nLimitations were not really discussed.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "YbTEtr1vzYy", "reviewer": "Reviewer_RgVr", "summary": "The paper proposes deep Markov factor analysis (DMFA) that uses the Markov process to capture the temporal dynamics in the fMRI dataset and maps the high spatial dimensions to a low-dimensional feature space. DMFA can cluster fMRI responses into “low dimensional temporal embedding” based on subject or cognitive state. The empirical studies on synthetic and real fMRI datasets illustrate that DMFA generates better performance in comparison with state-of-the-art techniques --- such as NTFA and HTFA.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The following are the major concerns and minor comments:\n\n1) My first question is about the prior and trail distribution(s). You can replace the combination of Normal distribution and RBF kernel with the Gaussian Process (Wu, 2019, J. Machine Learning Research) or even a Deep Gaussian Process (Roger Frigola-Alcalde, 2015, Thesis at University of Cambridge). The author(s) should explain why they chose these assumptions and how they can affect the performance of the analysis.\n\n2) The proposed method is compared mostly with the same type of factor analysis. I would like to see comparisons with other paradigms of analyzing using similar datasets. For instance, comparing the performance of the proposed method with the graph-based approach (e.g., Mingliang Wang et al., 2021, Medical Image Analysis) or functional connectivity-based ones --- e.g., Lei Sun et al. 2021 Artificial Intelligence in Medicine.\n\n3) The proposed method can be summarized in the form of an algorithm --- perhaps at the end of Section 3.\n\n4) In this paper, the notations are confusing. In the regular papers, scalers are denoted by small letters, vectors are defined with small letters (highlighted by bold), matrices are denoted by capital letters using bold. In this paper, they are a lot of conflicts. It is so hard to trace what is a set, a matrix, or even a distribution.\n\n5) The equations should have numbers. It is hard to trace them in the current format.\n\n6) In this paper, \"embedding\" is used to describe the temporal low dimensional space. Mathematically, embedding spaces have greater dimensions than the original spaces (not lower). For clarity, I suggest the phrase \"low-dimension feature space\" be used here instead.\n\n7) The brain images in Figures 2--4 and S1 should have value bars to make it clear what each color represents. Additionally, some subfigures do not have legends. The font size in these figures should also be increased.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Please see the main comment section.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "30"}, "extra_text": {}}, {"review_id": "REQiQBUeV8V", "reviewer": "Reviewer_mRP4", "summary": "The authors present deep Markov factor analysis (DMFA) to capture temporal dynamics in functional magnetic resonance imaging data. The method chosen by the authors is relevant since existing methods overlook the highly nonlinear and complex temporal dynamics of neural processes when factorizing their imaging data. The authors claim that their method (DMFA) is able to cluster fMRI data in its low dimensional temporal embedding, which enables validation of fMRI related hypotheses.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Reasons to accept:\n\n- If the authors claims are correct, the study is ground breaking in fMRI data analysis since the current models do not encode temporal dynamics. The voxel activations have non-linearity and inherent temporal dependencies, which are explored using the methods mentioned by authors.\n\n- The paper is well structured and explained. The theory behind Deep Markov Factor Analysis (DMFA) is explained perfectly.\n\n- The authors thoroughly explored the current state-of-the art related to the paper. The comparisons are clear in the text, table as well as the figures.\n\n- Limitations an ablation study are explained clearly.\n\n- Ablation study is clearly explained by the authors.\n\n \n\nReasons to decline:\n\n- Are the claims mentioned in following lines correct? The authors mentioned in line number 95-98 that RNN-based methods do not estimate any spatial factors, and are harder to\n\ntrain on long sequences. Is this verified on Convolutional RNN models (Wang et al., 2019) and LSTM/GRU models (Yan et al., 2019)? Are there no recent papers using transformer based models?\n\n- Better visualizations. Hard to interpret Figure 2 (b) and Figure 3 (Middle) images.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Summary:\n\nThe authors explored a relevant problem and can be widely used in neuroimaging to transfer high dimensional imaging data into low dimensional and interpretable representations.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "72"}, "extra_text": {}}, {"review_id": "9YL0YCKo53", "reviewer": "Reviewer_1fHu", "summary": "The paper Deep Markov Factor Analysis (DMFA): Towards Concurrent Temporal and Spatial Analysis of fMRI Data proposes a scalable approach to factorize the data of each subject into a product of a small number of spatial and temporal components.\nFollowing on the work of (Manning et al., Topographic Factor Analysis, Plos ONE 2014) on topographical factor analysis, each spatial component is given by a radial basis function. The temporal components are modeled from a set of latent satisfying the Markov property. The first latent is conditioned on the data class giving to the model clustering capabilities. Variational inference is used to derive a lower bound of the log-likelihood which is then optimized via stochastic gradient descent.\nIn a first experiment on resting state data (Autism dataset), DMFA is used to obtain subject specific atlases that allow to seggregate ASD and control subjects with an accutacy of 60 % (it is definitely hard to read).\nOn a second experiment (Depression dataset), MDD and controls are exposed to different musical stimuli and DMFA is able to seggregate well the data of different subjects without supervision. In some extend it gives some hint about the kind of stimuli and the conditions of the subjects but results are not so convincing (the authors are honest about this). \nIn a third experiment they show that DMFA gives better results than other methods in terms of held-out log-likelihood and prediction of the next time point given the history.\nThe code for DMFA is available. Pre-trained model are given and the code to generate the figures from pre-trained model works.\n", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "# 1. Originality\nThe model proposed by the authors uses the same spatial prior as in (Manning et al., Topographic Factor Analysis, PloS one 2014) and a temporal prior similar to (Krishnan et al. Deep kalman filters, 2015.).\nThe combination of these two techniques is new.\nI find that the related work section is well done overall and I find that the difference with previous work is clearly stated. In the section \"Factor Analysis in fMRI\" I would have added techniques related to IVA (Lee JH et al. Independent vector analysis (IVA): multivariate approach for fMRI group study., Neuroimage 2008), MultiViewICA (Richard Hugo et al., Modeling Shared responses in Neuroimaging Studies through MultiView ICA, NeurIPS 2020) or TensorICA (Beckmann et al., Neuroimage, 2005).\n\n# 2. Quality\n## 2.1 Sensitivity analysis needed\nThe proposed method has many hyper-parameters to tune (number and size of hidden layers, size of latent factors,  number of spatial maps to cite just a few). Some of these parameters vary from one experiment to the next such as  the size of hidden layers between synthetic and real experiments or the value of $\\sigma_0$ between the autism and depression dataset.\nIdeally, it would be nice to have sensitivity curves for all hyper-parameters but the fact that there are no sensitivity curves for the hyper-parameters that vary across experiments is problematic I think.\n\n## 2.2 Autism experiment\n1. Unless I missed something, the model is fit on the full dataset (train and test set) so that the time-courses of interest can be recovered in all subjects. I find this worrying because there is leakage of information from the test set. Ideally, the parameters of the model would be learnt on the train data and fixed when the test data are used. However with DMFA this is not possible because there is one set of parameters per subjects. I suggest here another experiment using naturalistic stimuli that I believe yields a cleaner pipeline: the authors could use the forrest dataset where subjects listen to an audio version of forrest gump (Hanke et al., Scientific data, 2014) to learn the spatial representations. Then they could use the second part of the dataset where the same subjects are listening to music of different genre and attempt to recover the music genre from the projected data.\n2. I find it quite surprising to use temporal ICA sources to build a functional connectivity matrix. These sources are typically independent and in most ICA algorithms, they are normalized so the functional connectivity matrix should always be identity unless I am wrong. Maybe it makes more sense to use spatial ICA sources to get the atlases and then project the data onto these atlases ?\n3. A standard baseline would be to use the 444 regions of the BASC Atlas (Bellec, Multi-level bootstrap analysis of stable clusters in resting-state fmri., NeuroImage 2010) to get connectivity matrices before SVM is applied on them (with python the atlas is easy to get using nilearn https://nilearn.github.io/modules/generated/nilearn.datasets.fetch_atlas_basc_multiscale_2015.html ). I fill it is a more natural baseline than the ones tried by the authors.\n4. The result for DMFA (without SVM) are missing when the full dataset is used. Why did the authors make this choice ?\n5. I would have expected to see the results of competitive unsupervised methods (NTFA, HTFA) as well as the ablation study done here so that we can compare the methods based on the decoding accuracy they yield.\n6. The authors state that competitive methods for the same task reach 69 % accuracy (l256). I believe it would be fair to display this information in Figure 3 so that we can actually see where DMFA stands (I believe DMFA is below 69 % isn't it ?)\n\n## 2.3 Depression experiment\nThere is clear evidence that the different subjects are in different clusters. However, what is interesting here in my opinion is whether the model can retrieve the condition of the different subjects (or some characteristics of the displayed stimuli) better than alternatives. About the condition of subjects authors say: (l284) \"MDD group data points are fairly populated in the center of temporal latent while control group are dispersed across latent space\". I find this statement a bit fuzzy. As it is, this experiment does not allow me to formulate a clear conclusion about whether the model works better than alternatives or not. \n\n## 2.4 Comparison with state of the art\nI am not sure I fully understood this part. \n1. This comment asks for more pedagogy. The likelihood is intractable. How can the authors evaluate it ?\n2. Leaving aside the problem of information leakage from the test set, I think the previous experiments on task fMRI or rest fMRI data provide more classical benchmark in neuroscience and more importantly yield results that are relevant in clinical applications (if a method can classify better than alternatives whether a patient is ASD or not this would really yield practical benefits).  Although the held-out likelihood / prediction of next sample suggests that DMFA is a better generative model, I find it hard to know whether these metrics are really reliable and if it is of any use in neuroscience.\n\n# 3. Clarity\n\nI find this submission extremely difficult to follow.\n\n## 3.1 Notation and organisation issues\n1. The notations become clearer in Table 1 but when reading section 3, I don't know the dimensions of the latent variables $z_n^F$, the domain and co-domain of function $\\mu_{\\theta}^W$ etc. The space in which any parameter leaves should be defined straight away at least in the text if it is not done directly in the formulas. This is really an important point as it makes the submission difficult to understand.\n2. The notations $\\mathcal{N}$ and $Norm$ are redundant unless I missed something. All of these notations should be properly defined in a notation section.\n3.  Are $\\sigma_{\\theta}^F$ matrices ? I believe so but I am not sure. Matrices, vectors and reals should have a different notation (using uppercase, bold and lowercase for example). Another example  l132 $g$ lives in $[0, 1]$ but in Table 1 it lives in $[0, 1]^2$. I understand these are typos but with clearer notations, it becomes obvious to know what is correct.\n4. I believe $\\sigma_0$ is never defined.\n\n## 3.2 Motives for modelization choices\n1. Having a common name for groups of parameters such as $\\theta$ or $\\phi$ may be practical but I find it confusing in the sense that it hides the fact that some parameters are shared and some are not. In general it is not always clear which parameters are shared or not. The notations give some hints but it would be better if this was explicitely stated in the text. In addition, I believe it is important to understand why these choices are made. For example, you assume that the spatial components are subject specific as well as the temporal components. These assumptions are not made in all models (for example in (Calhoun, A method for making group inferences from functional MRI data using independent component analysis.,  Human Brain Mapping 2001) the temporal part is assumed to be shared across views). Other parameters are shared (such as the neural network functions) and we might want to know why.\n2. Table 1 describes a particular architecture. It would be nice to understand how the authors end up chosing this particular one and what is the performance of other choices.\n\n## 3.3 Readability\n1. The figures are too small especially Figure 3 (right) or the legend in Figure 4 (a). In Figure 3 (right) this is particularly problematic as it is difficult to read the accuracy values. A grid might also help.\n2. Table 1 is difficult to understand. There is this layer column which goes down to 10 layers whereas if I understand correctly, none of the networks have more than 1 hidden layer.\n3. In order to facilitate the reading of the derivation of the ELBO it would have been nice to see a more detailed version in appendix.\n4. (l309) \"trained each with their default hyperparameter\": what are the defaults hyper-parameters for DMFA ? For other models ? According to which paper or which package ?\n\n## 3.4 Reproducibility\n1. l226 \"These three groups of sources are periodically activated according to some random weights\". The type of randomness should be precised.\n2. In the depression dataset there is no value for $S$. \n3. I find it difficult to reproduce the figures from scratch using the code as it involves downloading the full dataset and pre-processing it myself (in the README there is no script to get the dataset although there seems to be one to do the pre-processing).  What I would have hoped for is a small synthetic experiment that I can run in a short amount of time and is just there to show me that the code works.\n\n# 4 Significance\n1. Ignoring the fact that some information may leak from the test set, results displayed in Figure 3 are only marginally better than SVM after applying an atlas. For a method that takes 12 000 minutes to run, I think we need a marked improvement if we want other researchers to use it. \n2.  However, I would like to highlight that with fMRI data it is quite difficult to study the temporal structure when the TR is as high as 2s. I would suggest the authors to take datasets where the TR is much lower such as (Feinberg, et al. \"Multiplexed echo planar imaging for sub-second whole brain FMRI and fast diffusion imaging.\" PloS one 2010) so that they have a higher chance of seeing strong effects.\n\n\n# 5 Minor comments\nl157 The comment that amortized variational inference could lead to over-fitting is wrong in my opinion. Amortized variational inference is strictly less powerful than variational inference as in amortized variational inference the parameters are restricted to leave in the image of the parametrized function while in variational inference they are completely free.   ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Authors have a broader impact section in appendix.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "15"}, "extra_text": {}}, {"review_id": "6aiY9powGp", "reviewer": "Reviewer_m9hn", "summary": "Extending topographical factor analysis (TFA), this paper proposes a model, deep Markov factor analysis (DMFA), that models temporal dynamics in embedding space. The model is applied to both simulated and real fMRI data.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality\n\nThe proposed model is a novel extension of TFA with implementation details carefully worked out.\n\n\nQuality\n\nThe method seems solid, but the results are not as impressive. A few comments below.\n1. Is TFA always better than parcellating and averaging? Parcels are often generated from external datasets whereas latent representation are often built from datasets that are being analyzed. Would the former provide cleaner prediction evaluation?\n2. How is using Gaussian blob more interpretable than spatial ICA, which tends to generate biologically interpretable networks?\n3. The current simulation seems to be based on assumptions of the proposed model, and model hyperparameters are set to the ground truth. Also, quantitative comparison against TFA, HTFA, NTFA, and RKN should be provided.\n4. How were hyperparameters selected for real data?\n5. DMFA only beats ROI/ICA + SVM for some ABIDE datasets. Using p<0.02 as threshold would not survive multiple comparison correction for the number of datasets, i.e. 0.05/10 = 0.005. \n6. What model is used for estimating likelihood in Table 2? Is likelihood estimated per subject? What is the standard deviation? If consider the errorbar, is DFMA’s Nats significantly lower than other methods?\n\n\nClarity\n\nThe paper is clear and concise. A few minor comments below.\n1. Would ST-GCN fall under the scope of this paper? If so, please compare.\n2. More description on TFA would be useful for readers to get a better sense of the mentioned limitations.\n\n\nSignificance\n\nThe model would be useful for modeling any spatiotemporal data, but the results do not convincingly show better performance over existing methods.\n\n\nPost Author Responses\n\nI have read through the other reviews and responses, and will keep my original score mainly since classification improvements are only seen for half of the tested datasets.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "An elaborate discussion societal impact is provided in the paper. Limitation on implementation is also provided.", "review_text": "", "extra_scores": {"ethical_concerns": "There are no ethical issues in this paper.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Dear authors,

the overall consensus is that your work has merit in terms of method and experiments
although there are a number of concerns about some limited comparison with alternative
baseline approaches. Despite this concern I am inclined to flip the decision on the
positive side as code is available and you can still insert some numbers in the camera
ready as proposed.

Best regards,
The AC

Decision: Accept (Poster)