Multiscale Bayesian State Space Model for Granger Causality Analysis of Brain Signal
Sezen Cekic, Didier Grandjean, Olivier Renaud
TL;DR
This paper introduces a multiscale Bayesian state space model utilizing wavelet decomposition to analyze time-varying, frequency-specific Granger causality in brain signals, capturing complex neural interactions during tasks.
Contribution
It extends Bayesian state space models with wavelet-based multiscale analysis for dynamic, frequency-specific Granger causality estimation in neuroscience data.
Findings
Revealed complex frequency-based cross-talk between amygdala and orbitofrontal cortex.
Demonstrated effectiveness of multiscale Bayesian approach on intracranial data.
Captured short- and long-range dependencies in brain signals.
Abstract
Modelling time-varying and frequency-specific relationships between two brain signals is becoming an essential methodological tool to answer heoretical questions in experimental neuroscience. In this article, we propose to estimate a frequency Granger causality statistic that may vary in time in order to evaluate the functional connections between two brain regions during a task. We use for that purpose an adaptive Kalman filter type of estimator of a linear Gaussian vector autoregressive model with coefficients evolving over time. The estimation procedure is achieved through variational Bayesian approximation and is extended for multiple trials. This Bayesian State Space (BSS) model provides a dynamical Granger-causality statistic that is quite natural. We propose to extend the BSS model to include the \`{a} trous Haar decomposition. This wavelet-based forecasting method is based on a…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Multiscale Bayesian State Space Model for Granger Causality Analysis of Brain Signal
Sezen Cekic
Methodology and Data Analysis, Department of Psychology,
University of Geneva,
Didier Grandjean
Neuroscience of Emotion and Affective Dynamics Lab,
Department of Psychology,
University of Geneva,
and
Olivier Renaud
Methodology and Data Analysis, Department of Psychology,
University of Geneva
Abstract
Modelling time-varying and frequency-specific relationships between two brain signals is becoming an essential methodological tool to answer theoretical questions in experimental neuroscience. In this article, we propose to estimate a frequency Granger causality statistic that may vary in time in order to evaluate the functional connections between two brain regions during a task. We use for that purpose an adaptive Kalman filter type of estimator of a linear Gaussian vector autoregressive model with coefficients evolving over time. The estimation procedure is achieved through variational Bayesian approximation and is extended for multiple trials. This Bayesian State Space (BSS) model provides a dynamical Granger-causality statistic that is quite natural. We propose to extend the BSS model to include the à trous Haar decomposition. This wavelet-based forecasting method is based on a multiscale resolution decomposition of the signal using the redundant à trous wavelet transform and allows us to capture short- and long-range dependencies between signals. Equally importantly it allows us to derive the desired dynamical and frequency-specific Granger-causality statistic. The application of these models to intracranial local field potential data recorded during a psychological experimental task shows the complex frequency based cross-talk between amygdala and medial orbito-frontal cortex.
keywords: À trous Haar wavelets; Multiple trials; Neuroscience data; Nonstationarity; Time-frequency; Variational methods
The published version of this article is
Cekic, S., Grandjean, D., Renaud, O. (2018). Multiscale Bayesian state-space model for Granger causality analysis of brain signal. Journal of Applied Statistics. https://doi.org/10.1080/02664763.2018.1455814
1 Introduction
In many neuroscientific experiments, data are recorded in an experimental situation where stimuli are presented at fixed times and are expected to induce a reaction. For psychologists and neuroscientists, being able to model and explain the dynamics of the functional and effective links between neural and behavioural signals recorded during the experiment is of primary interest. They often have strong prior hypotheses about these causal links and therefore need reliable statistical tools to draw valid conclusions.
1.1 Granger causality
The question of how to operationally formalize and test for causality is a fundamental and philosophical problem. A mathematical solution, which relies on the causal nature of predictability, was provided in the 60’s by the economist Clive Granger and was latter coined “Granger causality”. According to Granger [25], if a signal “Granger-causes” a signal , then the history of should contains information that helps to predict above and beyond the information contained in the history of alone. The axiomatic imposition of a temporal ordering is the crucial element that enables us to interpret such dependence as causal. The presence of this relation between and will be referred to “Granger causality” throughout the text. In the 1960s, Granger [25] adapts the definition of causality proposed by Wiener [53] into a practical form and since that time Granger causality has been widely used in economics and econometrics. It is however only since last few years that it became popular in neuroscience, see Cekic et al. [12] for a recent review.
1.2 Existing methods and limits
In the context of linear Gaussian autoregressive models, for which the restrictions are Gaussianity and linearity, which imply stationarity in most cases, a common way to test for Granger causality between two series is to estimate a vector autoregressive model (VAR) and then test the significance of the off diagonal coefficients of interest [26, 34].
However, in neuroscience the data are usually nonstationary and this characteristic is moreover of interest. In the simplest case, the data are stationnary up to a particular point where the Gaussian process has been perturbed away from its stationary distribution (perhaps by some external intervention). We therefore want to derive a causality-statistic that allows us to capture the causal structure differentially for each time.
Basic causality statistic has therefore to be extended to the nonstationary case in order to be suitably applied in a neuroscience context, which can be achieved by letting the VAR model evolve in time.
Time-varying VAR model estimate implies three challenges: over-parametrization, model order selection and multiple trials.
The two widely used approaches allowing us to deal with the nonstationarity are the windowing approach, based on the locally-stationary assumption and the adaptive estimation approach, based on the slowly-varying assumption of the parameters (see Ding et al. [17], Schlögl [44] and [13]).
The windowing approach consists in estimating VAR models in short temporal sliding windows where the underlying process is assumed to be (locally) stationary (see Ding et al. [17] for a methodological tutorial of windowing estimation approach in neuroscience).
The windows size is a trade-off between the accuracy of the parameter estimates and the resolution in time. The choice of the model order is a delicate issue, and depends on the choice of the segment length. Some criteria have been proposed in order to optimized simultaneously the windows length and the model order [32].
The usual approach with multiple trials is to average the estimation or do a global optimization to get an overall estimation.
In Cekic [10] we found that this windowing methodology presents several limits. First, the improvement of the time resolution implies short time-windows and so few residuals for assessing the quality of the fit. In addition the size of the temporal windows is subjective (even if it depends on a criterion) as is the overlap between the time-windows. The order of the model in turns depends on the size of the windows and so the quality of the estimation strongly relies on several subjective parameters.
The adaptive estimation approach consists on estimating a different model at each time, where the observations at time are expressed as a linear combination of the past with coefficients evolving slowly over time. The differences between the methods consist on the way the transition and the update from coefficients at time to those at time are processed [see 44].
All these adaptive estimation methodologies depend on a quantity that acts as a tuning parameter and defines the relative influence of past values and innovation noise on the recursive estimation. Generally this free tuning parameter determines the speed of adaptation as well as the smoothness of the time-varying VAR parameter estimates. The algorithms are very sensitive to this tuning parameter [44] and therefore the estimation quality strongly depends on it. The “ad-hoc” nature of this tuning parameter is obviously a problem in term of statistical inference and this issue was not raised in the development of these algorithms coming from the engineering field. The model order and the tuning parameters are usually optimized together by a Mean Square Error criterion [45].
Kalman filtering algorithm [29] can be used in order to estimate time-varying VAR models when we express it in a state space form [1, 4].
If the transition matrix and variance-covariance matrices of the observed and state equation are known, the Kalman smoother algorithm gives the best linear unbiased estimator for the state vector [29], which in this specific case contains the time-varying VAR coefficients. In the engineering and neuroscience literature, these matrices are systematically set to fixed values or estimated through some “ad-hoc” estimation procedure (see Schlögl [44] and Hesse et al. [27], Arnold et al. [1] for applications in neuroscience). There is moreover the very important issue of model order selection which becomes very tricky with model complexity and the plurality of the trials.
1.3 Neuroscience data specificities
Our model has to be applied to experimental neuroscience data, whose intrinsic specificities must be taken into account in its development. Therefore, in order to derive a suitable dynamical causal statistic, we need a model that allows us to get a reliable estimate of the dynamical VAR coefficients based on several trials and, last but not least, that also allows us to capture short- and long-range causal dependencies between signals due to specific frequency characteristics of the data.
1.4 Proposal
Faced with data with a time-varying structure (like neuroscience data), none of the above methods relies on a tailored statistical model that provides satisfying estimation and inference procedures and proposes a solution to deal with short and long range causal dependencies potentially present in the data.
We propose a new methodology for suitably modelling multivariate nonstationary time series in order to get a reliable Granger-type dynamic causal statistic. It is based on a linear Gaussian vector autoregressive (VAR) model with coefficients evolving over time according to a linear dynamical system. Given that this model is strongly over-parametrized, we propose to place it in a Bayesian framework and to use the variational method to estimate all the densities. This variational Bayesian methodology [6] estimates all the necessary quantities. The Bayesian nature of the model moreover offers a natural criterion for model order selection. In Section 2, we describe our Bayesian state space (BSS) model and discuss its technical specificities and the estimation procedure. We extend it to deal with multiple trials (or epochs), in a proper manner for the estimation and the inference procedure (section 3.3). In Section 4, we propose an additionnal extension of the BSS model, called the multiscale Bayesian state space (MSBSS) model, which is based on the à trous multiscale wavelet transform. The latter approach has never been used in this context of time-varying VAR coefficient estimate and we will show that it allows better estimate of short- and long-range specific dependencies between recorded signals and offers a very simple way to deal with time-frequency uncertainty bounds. In Section 5, we will present a Bayesian dynamical Granger-causality statistic based on the time-varying estimated VAR coefficients and in Section 6, we present simulation studies to assess our proposed methodology. Finally, we present in Section 7 the results of the application of the method to intracranial local field potential recorded during a psychological experiment in specific brain areas, namely in the amygdala and orbitofrontal cortex.
2 The Bayesian State Space Model
We propose to write the dynamic VAR model in a state space form with an observation equation in which the dynamic VAR coefficients are driven by the state equation [this modelling proposal was previously made by 9]. This leads to the following system of equations:
[TABLE]
where is the value of the signals at time , are the time-varying VAR coefficients (up to order ) and the vector of size contains all the time-varying VAR coefficients that have to be estimated for the time . The matrix is the transition matrix of the state vector , is the variance-covariance matrix of the state equation, and is the variance-covariance matrix of the observation equation. We will use the notation and to denote the entire set of values from to . Although the state equation seems to be only a first-order autoregression, note that the vectors and contain all the coefficients up to order , and therefore the direct dependency of on past values is actually unlimited (and driven by the choice of the order ). This formulation is actually similar to the state-space representation of AR(p) or ARMA(p,q) models (see e.g. [4] examples 12.1.4-5)
With a slight abuse of notation, we can write and , which are respectively the observation, the state and the initial state densities.
Given the Bayesian framework, we are interested in obtaining the posterior distribution of the unknowns of this model. In simple cases, these distributions can be obtained analytically but here, the curse of dimensionality makes the computation intractable and forces us to rely on approximation techniques.
2.1 Variational Bayes
Variational approximation was applied to the linear Gaussian state space model in Ghahramani and Beal [21] and Cassidy [9], and more recently in Luessi et al. [33]. This approximation methodology has been extensively explored and used during the past years; see for example Titterington [49], Beal [6], Ormerod and Wand [38] and Fox and Roberts [18]. It is not widely known within the statistical community dominated by Monte Carlo and Laplace approximation methods but is however much faster than Monte Carlo, especially for large models, and allows us to deal with models containing a very large amount of parameters (see Friston et al. [19] for a comparison of variational and Laplace approximations).
For sake of clarity, we will define the set of unknown parameters as , so the full set of unknowns for model (1) is .
The target quantity is the posterior distribution of the parameters which is an intractable integral of very high dimension. The variational approach allows us to approximate this posterior density by a variational posterior density that will be selected to be optimal according to the Kullback–Leibler distance dissimilarity criterion [6] (see Appendix A for further details on variational Bayesian approximation).
2.2 Mean-field approximation
Variational Bayesian methodology allow the approximating density to factorize over groups of parameters. The researcher has the choice here, but the explicit link between the ’s in equation 1 urge not to factorize . The less intricate links between , , and allow for the following factorization (see Figure A.1 from the supplementary material)
[TABLE]
Note that it may lead to serious degradation in the resulting inference if unsatisfied [38, 49, 6].
2.3 The variational evidence lower bound
As presented in details in Appendix B, the variational Bayesian methodology provides a quantity F that is a lower bound for the evidence of the model and that can be computed efficiently. This leads to a natural criterion for model order selection which is crucial to estimate our strongly over-parametrized time-varying VAR model. We can therefore perform model selection by comparing the F quantities computed for each model order and select the model that exhibits the highest . The specific analytic form of for our model is derived in Appendix B.
3 Model specification
We will now describe the fully hierarchical Bayesian model that we propose.
The optimal form for and of course depends on the choice of the prior distributions and . However, the analytical form of will be of the same distributional form as the prior distributions if the complete-data likelihood is part of the exponential family and if the hidden and parameter prior distributions and are conjugate to this complete-data likelihood (this condition is known as “conjugate exponential”, see Beal [6]).
Recalling the model (1), the complete data likelihood is Gaussian. We will now describe the model by defining conjugate priors over the model parameters.
3.1 Prior distributions
3.1.1 Prior for
The state space structure of the model (1) leads to the conditional prior distributions for the hidden state defined in Section 2 and the prior mean and variance for are set to and respectively.
3.1.2 Prior for
We propose a diagonal structure for the matrix, meaning that a specific causal coefficient at time will only depend on its own past value plus a white noise and not on the past values of other VAR coefficients. This assumption seems reasonable in terms of brain connectivity, where the dynamic of a particular causal input can be assumed to be driven by its own trajectory and not by that of other causal or auto-causal inputs. To satisfy the conjugacy condition explained in details in Appendix A, We place a Gaussian prior on each diagonal entry of the matrix:
[TABLE]
where is the hyperparameter variance of the specific diagonal element whose distribution will be discussed below. We choose to impose a conservative prior for by setting each location hyperparameter to .
As we have no subjective input for the variances , it is desirable to have priors that exhibit very little information. The main recommendation in Gelman [20] is to use half-t priors on standard deviation parameters to achieve arbitrarily high noninformativeness. The definition of the half-t distribution can be found e.g. in Wand et al. [50]. They showed that the half-t distribution can be written as a scaled mixture of inverse-gamma distributions. So the half-t prior is obtained for each element through the auxiliary variable construction
[TABLE]
which ensures that . The hierarchical representation in equation (4) respects the conditional conjugacy condition due to the conjugacy properties of the inverse-gamma distribution.
By its noninformativeness property, this half-t prior choice will let each variance element get an appropriately high or low value, thereby allowing them to play the role of shrinkage parameters for the distribution of . In our specific case, it will allow us to be conservative and to tend to avoid an erroneous causality assessment.
3.1.3 Prior for
The accuracy of the Markovian conditional distribution followed by the VAR coefficients is established through the variance-covariance matrix . Throughout this work we assume that this matrix is diagonal (for the same theoretical reason discussed for the matrix) and that all its elements are equal. In fact there is no reason to assume that the variability of one VAR coefficient is different from that of another. This choice is also motivated by the concern to limit the number of parameters to estimate.
In order to satisfy the conjugacy conditions and for the same theoretical reason mentioned earlier, we set the same weakly-informative half-t prior distribution for the single standard deviation parameter of the diagonal matrix:
[TABLE]
which again ensures that .
3.1.4 Prior for
The variance-covariance matrix of the observed equation is not supposed to be diagonal due to the theoretical interdependence between brain signals modelled in the system.
We impose for a generalisation of the multivariate case of the Half-t prior used for and . Huang et al. [28] derived this prior and explained that with a suitable hyperparameter choice, it induces an Half-t distribution for each standard deviation term corresponding to the diagonal of , as well as a marginal uniform distribution for all correlations. Explicitly we have
[TABLE]
where denotes a diagonal matrix. This hierarchical structure ensures that each diagonal element is distributed as , and the particular choice of leads to marginal uniform distributions over for all correlation terms [28].
Lastly, we note that the selected prior distributions impose the choice of hyperparameters , , and . We will choose for the uniform property to hold. For the remaining hyperparameters, the larger they are the less informative the priors are. We follow Menictas and Wand [36] and set them to .
3.2 Update equations
The variational Bayesian iterative algorithm and optimal posterior distributions are based on the results presented in Appendix A and factorization chosen in Section A.0.2. The optimal variational posterior distribution for the hidden state sequence (Variational E-Step) is multivariate Gaussian at each time :
[TABLE]
where the sufficient statistics , as well as the cross-moments
, are obtained by the Kalman–Rauch–Tung–Striebel (KRTS) smoother algorithm [29, 30], applied to an augmented system of equations, see Appendices C, D and E. All the derivations for the Variational M-steps can be found in Appendices F, G, H, I, J, K and L.
3.3 Multiple trials
One important contribution of the present article is to show that our model can be modified to deal with conditionally independent sequences which are supposed to have the same hidden state. This reflects the case that arises during an event-related experimental paradigm, where many trials on the same condition are measured.
In Beal [6] and Cassidy [9], this extension is treated by first estimating the necessary sufficient statistics for each sequence independently in the E-step, and then by averaging these statistics to get only one set of sufficient statistics, which is representative of the entire set of independent sequences before performing the M-step. However this approach does not take into account the complex dependence of the variational posterior on the whole dataset . In fact, all the computations done so far can be adapted to multiple trials. A detailed derivation of the model and related variational posteriors distributions in a multiple trials setting can be found in Appendix M.
4 The Multiscale Bayesian State Space Model
In a neuroscience context, an important limitation of many models, including the BSS model, is the inability to capture both the short- and long-range possible causal dependencies between signals. Causal interactions in a neuroscientific experiment context may not be instantaneous, but delayed over a certain time interval that must be subjectively chosen depending on the research hypothesis. Another important subjective parameter is the time-lag (), that determines the interval between two data points As shown in Barnett and Seth [5] and Solo [48] these choices strongly influence classical Granger statistics in their ability to detect causalities, whereas the result should be as invariant as possible to arbitrary choices like the chosen sampling frequency or added time-lag in the prediction model.
Based on these considerations, and because Granger causality is based on predictive ability, if the auto-causal information contained in the history of the predicted signal in model (1) is not well represented in , the Granger-causality evaluation, which is based on the predictability improvement of by adding the information contained in the history of the second variable , may be spuriously assessed as significant. On the other hand, if the predictive ability of the history of the causal signal is not informative enough, we can miss some crucial information about causal interdependencies between signals and .
We propose a new solution that has the ability to appropriately select the short- and long-term causal histories of and by combining the BSS model with the à trous multiscale decomposition methodology and that remains within the conjugate exponential framework [41].
4.1 The à trous Haar wavelets transform
For that purpose, we will perform the à trous multiscale decomposition of signals and contained in in model (1), in order to use these quantities as predictive histories in the matrices . The reader is referred to Renaud et al. [41] and references therein for a complete overview of the method. We define as the à trous wavelet coefficient of the signal at time for scale . Since we will use the à trous wavelet coefficients for prediction, they should not be based on future values and the only family that satisfies this constraint is the à trous Haar wavelet transform. Then , where the finest scale is the original series .
4.2 The multiscale Bayesian state space model
Based on the derivation in Section 4.1, we can modify the VAR model in equation (1). We keep the quantity to predict in the time domain, but the histories of the series and will equal the à trous Haar wavelet transforms of series and respectively. We can therefore define the set that contains the decompositions of the histories of series and as
[TABLE]
and therefore adapt equation (1) replacing with .
In particular, we underline that for any , for a suitable choice of , is an orthogonal transform of and therefore in this specific case BSS and MSBSS models are equivalent. For any choice of the ’s, all results in Section 2 and all the estimation procedures in Section 3 can be applied. The only difference is that the matrix is replaced by , and that the dimension of is changed accordingly. The model order for the different scales as well as the number of scales to be taken in the model must now be selected. As argued in Renaud et al. [41], the relative non-overlapping frequencies used in each scale motivate an independent selection of the model order for each scale. Defining a maximal scale decomposition and a maximum model order for each scale , the model order selection procedure set in Section A.1 is iteratively applied to select in a stepwise manner for each scale from 1 to . The free-energy quantities related to the models are then compared and the that exhibits the highest free energy is selected. The model order for the smooth is finally selected.
This à trous extension can thus be viewed as a generalisation of the BSS model that contains information relative to the frequencies. Each wavelet scale indeed is directly related to a specific frequency band and the resulting dynamic Granger causality statistic is therefore directly interpretable in terms of frequencies (see Section 7).
5 Bayesian Granger-Causality Statistic
Based on the model proposed in this article, a necessary and sufficient condition for not to be Granger-causal for at a given time , is that each element of the subset of , that contains all the causal coefficients of interest, equals zero.
The most appropriate approach to evaluate the compatibility of this type of hypothesis with the data would be to compute for each time a Bayes factor between the VAR model under the restriction for just one value of () and the VAR model without restrictions (). For each , this would requires the computation of the evidence of model , which seems untractable: one would need a Markovian process for that is conditionned (sort of bridge) on the restriction that for just one given . At the very least, one would need to estimate a different (conditional) model for each and compute its free energy. Additionnaly, we have no guarantee that the free energy approximation is of the same magniture for all these models.
We will use a simpler approach that rely only on the posterior density of the (unconditional) model. The use of highest posterior density (HPD) regions for Bayesian testing was introduced in Box and Tiao [7] and used in the Bayesian literature, as for example in Kim and Press [31] and West and Harrison [52, p. 280]. For a given time , consider the sub-vector and let be its dimension. A suitable partition of the posterior parameters and gives the marginal posterior density . The contribution of each element of to the prediction of may be assessed by considering the compatibility in this marginal posterior with the value , where is a zero vector of dimension . The question is therefore to know whether the parameter point is included in the highest posterior density region (HPD) of size . This happens if and only if , where is the quantile of the standard distribution with degrees of freedom. As stated in Box and Tiao [7, p. 125], it follows that the parameter point is covered by the HPD region of content if and only if . Equivalently, we can search the region of minimum coverage that contains .
The above approach gives only pointwise evaluations (i.e., for a given time ). When jointly testing a set of values for a complete time, frequency or time-frequency connectivity map, it is important to suitably correct the significance threshold for multiple comparisons. We do not correct it in Section 6 as we are interested in the separate evaluation for each time but we do for the application in Section 7.
6 Assessment of Accuracy
We now turn our attention to the accuracy of variational Bayesian inference for our model (1) under the mean-field assumption described in Section A.0.2 and priors discussed in Section 3. We provide here a study of the quality of variational Bayesian estimate for model order selection and Granger causality detection. The first simulation study in Section 6.2 evaluates the ability of the BSS and MSBSS models to detect Granger causalities between two signals and in Section 6.2, we present a systematic comparison between the proposed BSS and MSBSS models and the windowing approach in term of Granger causality detection ability. Note that a Monte-Carlo type of simulation is not feasible due to the very large amount of unknown parameters of our model and so direct comparison of the obtained variational posterior with true posterior density was unfortunately not feasible.
6.1 Practical implementation
To initialize the BSS and the MSBSS model algorithms, we run 10 iterations of the simple EM algorithm [47]. We thereby obtain reliable starting values for under initial conditions defined in Section 3. We then iteratively update the parameters of the variational approximate posteriors through the Variational Bayes EM until the relative free-energy criterion described in Section A.1 between two consecutive iterations changes less than a tolerance value that we choose here to be equal for the model order selection study and to for the Granger causality detection study following Menictas and Wand [36].
6.2 Granger-causality detection
We will now assess the ability of the proposed models BSS and MSBSS to detect Granger causality. We simulated signals with parameters that vary slowly in time, which is a reasonable simulation of neuroscientific data. It will show how the method performs with data that are not generated according to the model. The simulation consists in replications of a bivariate BSS model but with hidden variables that evolve slowly and deterministically through time. Based on these deterministic , data are generated from the observed equation in model (1). The simulations were carried out for model orders , series lengths , number of trials and also for different values of the causal parameter in order to see the limits of the methods. The simulations for the causal parameter values were carried out for all model orders and number of trials, but only for a series length of time points (which is the case where the method will break first).
The matrix was fixed to the values and the slowly-varying parameters are all set to zero except for the entries related to the causal parameters and for order . The values of the simulated parameters can be seen on the two panels (which are identical) on the top of Figures in Appendix N. Signals were simulated with normal and non-normal errors. Additional simulation settings and additional results are available in [11]
Data generated with slowly-varying parameters and normal errors
Concerning the MSBSS model, for each of the simulation, the model order for the different scales as well as the number of scales to be taken in the model are selected as described in Section 4. The maximum number of scale was set at and the maximal model order per scale at . In multiple trial scenarios, the number of scales as well as the model order selection procedure was carried on one single trial only (for computational simplicity). Estimation procedure is always conditioned on the same number of time points, allowing us to compare the models. The Granger-causality detection capability of a method is quantified by the true negative rate (TNR) and the true positive rate (TPR). For each time, TNR is the percentage of the simulations for which the % HPD region defined in Section 5 contains the causal statistic when it must actually contain it, and TPR is the percentage of the simulations for which the % HPD region does not contain the causal statistic when it should not. Secondly, the BSS model is estimated with the true model order (oracle order) and the relative TNR and TPR are computed for that model. Finally, for the BSS model, the model order is selected with the free-energy criterion as described in Section A.1 and the model is estimated based on this selected order . The relative TNR and TPR are calculated.
The causality detection ability for the BSS model with estimated model order is not shown in the results, because the model order was selected correctly by the free-energy criterion for each simulation (the causality detection ability is thus exactly the same as this for the BSS model with oracle order).
Figures 1 and 2 show Granger-causality detection accuracy for a model order , series length of and causal parameter equal to and respectively. By construction, the causality arises between times , as shown in the two graphs on the top of the figures. For Figure 1, the BSS model performs slightly better than the MSBSS model with trial and both methods yield very similar results, correctly recovering the underlying directional influences with trials. One can observe that the causality detection evolves as the causal parameter changes from zero to one. The TPR is indeed gradually increasing as the parameter changes from zero to one, slightly faster for the BSS than for the MSBSS model for , and identically for . For Figure 2, with such a low causal parameter, the TPR related to the MSBSS model is worse than the one related to the BSS model for trial, but both method perform very well with trials although the causal parameter is low.
Figures for Granger-causality detection accuracy related to all other model orders, series lengths and causal parameter values for normal errors are given in Appendix N. Globally, for a causal parameter that equals , the MSBSS and the BSS models show very good Granger-causality detection accuracy in terms of TPR and TNR for and trials. The only cases which display poorer results are the causality detection for the MSBSS model with trial for model order and series length , and . When the causal parameter value decreases, the causality detection ability decreases as well. Poorer causality detection results are globally observed for a causal parameter from to for trial and for a causal parameter value of for trials, especially for the MSBSS model.
Data generated with slowly-varying parameters and non-normal errors
We also simulated data with the same settings as above, except for the observation equation errors which are here multivariate t-distributed with parameters and . All graphical results can be found in Appendix N. Globally, the detection accuracy remains satisfactory with data generated with non-normal errors. The Granger-causality detection ability of the proposed methods is therefore robust to this model assumption departure.
Comparison with the windowing approach.
We will now compare our proposed models (BSS and MSBSS) with the windowing approach. The comparison is performed on data generated with slowly-varying parameters with normal errors as explained in Section 6.2. Simulations are performed for model orders , series length , number of trials and causal parameter values . Data are fitted using the sliding window methodology proposed in Ding et al. [17]. The overlap parameter between the time-windows is chosen to be equal to time point and the subjective windows size is chosen to be equal to time points. The model estimation procedure is performed using the Viera–Morf algorithm implemented in BSMART [15] and GCCA toolboxes [46]. A first estimate was obtained with the true model order (oracle). For a second estimate, the model order selection is performed by considering the mode of the optimal model order calculated in each temporal window based on the Bayesian Information Criterion (BIC) as proposed in the BSMART toolbox [15], in the GCCA toolbox [46] and in the SIFT toolbox [37]. Based on the estimated models, a time-domain Granger-causal statistic is computed in each temporal window and its significance is assessed through the asymptotic distribution (see [15] and [46] for further details).
In Figures 1 and 2, the dash and the dot lines represent the results for the windowing approach. For Figure 1, the MSBSS and the BSS models perform much better than the windowing approach with the true model order and with the model order selected based on the BIC criterion for trial. MSBSS and the BSS models detect the causality faster than the windowing approach with trials. The TPR is worse for the windowing estimation than for the MSBSS and the BSS models for both and trials. For the case with a low causal parameter (Figure 2), the windowing estimation with oracle order performs almost identically as the MSBSS model and performs much worse than the BSS model for trial. MSBSS and the BSS models detect the causality faster than the windowing approach with oracle order with trials. The windowing estimation with the model order selected based on the BIC criterion display very poor detection accuracy for and for trials.
All other graphical results related to the windowing estimation procedure can be found in Appendix N. Globally, for the simulations with trial, the model order selection procedure based on the BIC performs well for order and is inaccurate for orders and and the related Granger-causality detection fails. One can observe that the causality detection evolves as the causal parameter changes from zero to the maximal value of the causal coefficient. For the simulations with trial, the windowing estimate detects the causality pattern globally slower than the MSBSS and the BSS models, whereas the TNR remains almost identical. When the value of the causal parameter decreases, the TPR for the MSBSS model decreases as well and becomes worse than the TPR for the windowing estimate in some cases (e.g., for the causal parameter and model order ). The TPR for the BSS model estimate decreases with the value of the causal parameter as well, but its TPR is always better than the TPR for the windowing estimation procedure. When the causal parameter equals , however, the MSBSS and BSS model estimates do not detect anything whereas the windowing approach displays a TPR around %.
For the simulations with trials, the model order selection procedure based on the Bayesian Information Criterion performs well overall. Windowing estimate (with model order selected by the BIC and oracle model order) globally detects the causality pattern slower than the MSBSS and the BSS models, whereas the TNR is almost identical everywhere. When the value of the causal parameter decreases, the model order selection procedure based on the BIC performs less well and its TPR degrades even more than the one for the MSBSS and the BSS models (e.g., for order and causal coefficient and and for order and causal coefficient , and ).
7 Application
We will now apply the proposed MSBSS model to real iEEG (intracranial electroencephalogram) data recorded during a psychological experimental situation. Brain recordings are localized within the amygdala (AMY) and medial orbito-frontal cortex (mOFC) regions in order to study the dynamics of neuronal processes between these regions in response to emotional prosody exposure. It is known in the literature that the emotional content of the stimulus induces the presence of causal links AMY mOFC and mOFC AMY (Grandjean et al. [23],Grandjean et al. [24],Grandjean and Scherer [22]).
7.1 Results
Here we present the results for one patient and the experimental condition anger. The patient was exposed to short pseudowords pronounced with an angry prosody (Grandjean et al. [23]). The data were acquired with a high resolution of Hz during sec. There are trials available. As high frequency behaviour of data was not of interest, we downsampled the data by a factor of four with an exponential smoother (it is of course vital to process with a smoother that is not based on future values). A MSBSS model was estimated following the procedure described in Section 4. As researchers were already aware about frequency bands where they expect a causal relationship, the number of scales was fixed to (on the downsampled data), corresponding to four frequency bands respectively around , , and Hz, plus a smooth that represents the frequency content below Hz with a maximal order per scale of . These particular choices for are motivated by the occurrence of the onset of the stimulus at ms. Indeed, due to the use of past values in the model, the estimation procedure starts at time and we do not want to miss the onset of the stimulus in the model estimation. The model order per scale is selected on the whole set of trials.
Based on the MSBSS model, a significant Granger causality from signal 1 to signal 2 at frequency and at time means that the energy of signal 1 at frequency significantly improves the prediction of the value of signal 2 at time .
The testing procedure was performed as follows: we first tested the overall set of causal VAR parameters of interest and then the scale-specific causal VAR parameters () by considering the compatibility in their marginal posterior with the value , where is a zero vector of dimension . We computed the coverage of the smallest highest posterior density (HPD) region that contains . We will call the complementary probability of this coverage the significance level.
The results are reported for the two directional causalities AMY mOFC and mOFC AMY. The first line of each graph in Figure 3 represents the results related to the overall statistic and each following line represents the results related to each scale corresponding to the frequency bands respectively around , , and Hz.
The number of tests provided by the proposed method is basically proportional to the number of time points. To circumvent the multiple testing problem, a solution that seems suitable is the cluster mass test which consists in defining clusters of neighbouring time regions using a permutation scheme to assess its significance. We applied it for the overall testing for each frequency with a threshold corresponding to a level of and this procedure controls for the family wise error rate (FWER) [35].
7.1.1 Testing the scale specific causality
Figure 3 shows the results of the estimated Granger causality for the two causal links of interest AMY mOFC and mOFC AMY. Vertical lines represent the onset and the offset of the stimulus which occur respectively at ms and ms. They are displayed for the interpretation of the results but were not provided to the model. The estimated model order for each scale plus the smooth is [5 5 3 1 1]. Only the clusters assessed significant by the cluster mass test are displayed and the colors represent the value of the individual statistic. As mentioned in Section 4, the smooth represents the frequency content of the series from the largest scale to the lowest frequency in the signal.
These results give a partial answer to the question of how AMY and mOFC regions are functionally and causally related during the exposure of auditory emotional stimuli. Actually, an initial Granger causality event from mOFC AMY is observed in gamma range (Hz) just before the stimuli onset and during stimulus exposure followed by an AMY mOFC in beta range (Hz) during this same period. At the offset of the stimuli the Granger causality is again from mOFC AMY in gamma range (Hz), this directional Granger causality is sustained during the period after the stimuli presentation. These results are compatible with a known complex cross-talk between these two brain regions during emotional stimulus exposure and the related meaning of such event for the organism. Of course such results should be extended to several patients with similar brain recordings in targeted brain regions and compared to the processing of other emotional stimuli.
8 Conclusion
We derived a time-varying Granger-causality statistic through a Bayesian nonstationary multivariate time series model with dynamic coefficients. While similar models have been proposed [9], one of the main contributions of this article is to provide an assessment of accuracy of the method, and especially an extension to the à trous Haar wavelets transform [41]. This very flexible à trous Haar procedure enables us to capture short- and long-range dependencies between signals with only few parameters to be estimated and allows us to be specific for the frequency in the assessment of causality, which is a main point of interest in the neuroscience community. The central finding of this article is that variational Bayesian estimate of time-varying VAR models enables us to achieve good (if not excellent) accuracy in terms of Granger-causality detection for normal, but also for non-normal errors and that the method performs much better than the commonly used windowing methodology in terms of Granger-causality detection accuracy. This model thus provides a very powerful tool for dynamical spectral causality analysis in a neuroscientific context.
Further points to highlight are the suitability of the model to deal with multiple trial in a fully correct way and the potential extensibility of the model to deal with more than two signals, or two sites of interest.
The choice of the sampling frequency (e.g. for EEG recordings) and the preprocessing steps commonly perform by researchers in neuroscience are two delicate issues for subsequent Granger-causality analysis, because different choices may lead to different Granger-causality results. Due to the multiscale decomposition, the MSBSS model is probably much more robust in this regard, and this would be an interesting topic for further research.
Toolbox and Appendix
An open matlab toolbox called MSGranger is available at the following url: https://www.unige.ch/fapse/mad/services/matlab/msgranger/. Single and multiple trials are implemented and both the Bayesian State Space (BSS) and the multiscale Bayesian state space (MSBSS) models can be used to obtain dynamical and frequency-specific Granger-causalities.
Appendices with all the derivations, Figures referenced in the text, as well as the cluster mass test and data analysed in Section 7 are available in the appendix.
Acknowledgements
The authors gratefully acknowledge Center for advanced modelling science, Swiss National Science Foundation under Grant 100014_156493 and Lifebrain H2020-SC1-2016-201 under Grant 732592.
Appendix A Elements of variational Bayes
The most common type of variational Bayesian methodology, known as mean-field approximation uses the Kullback–Leibler distance (KL distance) between and as a dissimilarity function.
A.0.1 Learning rules
By simple algebra, we will decompose the evidence of the model by inserting the variational density . Using the fact that , we have
[TABLE]
and therefore
[TABLE]
where denotes expectation and its subscript denotes the density used for this expectation. Equation (LABEL:VB_freeNRJ) is the fundamental equation of variational Bayesian methodology. By necessary positiveness of the KL distance, we have obtained a lower bound for the logarithm of the evidence
[TABLE]
where \text{F}\big{(}q(\varphi_{1}^{T},\Omega_{1}^{b}|Z_{1}^{T})\big{)} is called the negative free energy. By minimization of the KL distance between and , we maximize \text{F}\big{(}q(\varphi_{1}^{T},\Omega_{1}^{b}|Z_{1}^{T})\big{)}. Furthermore, since the KL distance is equal to zero if and only if the two densities and are identical, the functional quantity \text{F}\Big{(}q(\varphi_{1}^{T},\Omega_{1}^{b}|Z_{1}^{T})\Big{)} equals the model evidence if and only if equals the true target posterior . The aim is thus to find a density for which the integrals in \text{F}\Big{(}q(\varphi_{1}^{T},\Omega_{1}^{b}|Z_{1}^{T})\Big{)} are tractable and which is close to .
A.0.2 Mean-field approximation
The choice underlying the variational Bayesian methodology, known as mean-field approximation in physics, is to allow the approximating density to factorize over groups of parameters [see 19, for a comparison between Laplace and variational Bayesian assumptions]. We will suppose here that the approximating density factorizes as
[TABLE]
and the lower bound for the model evidence can be rewritten based on this factorization as
[TABLE]
Depending on the model at hand, mean field approximation may have minor to major impacts on the resulting inference. For example, if is such that and have a high degree of dependence, then the restriction = will lead to serious degradation in the resulting inference [38, 49, 6]. The factorization in equation (12) is obviously not unique. For instance, some authors factorize also into [51].
A.0.3 Variational EM algorithm
With the use of the calculus of variations (hence the name variational Bayes), it can be shown that under assumption (12), the variational distributions and , that maximize the functional \text{F}\Big{(}q(\varphi_{1}^{T}|Z_{1}^{T}), q(\Omega_{1}|Z_{1}^{T}),\dots,q(\Omega_{b}|Z_{1}^{T})\Big{)}, can be expressed and therefore maximized in an iterative way [6, 18]. First,
[TABLE]
where superscript denotes the iteration number. The other steps for are
[TABLE]
where is the expectation over all the distributions at iteration except . See Beal [6] and Ostwald et al. [39] for all proofs.
The distributions and are known as full conditionals in the MCMC literature. The mutual dependence of the optimal variational posterior densities in equations (14) and (15) suggest a similarity with Gibbs sampling [8] which involves successive draws from the full conditionals. Mean field approximation indeed leads to tractable solutions in situations where Gibbs sampling is also applicable.
Thus, if we set equal to and equal to , we have maximized the lower bound for the model evidence F under the (12) constraint.
The form of equations (14) and (15) define a circular dependence which explains the use of an iterative algorithm whose convergence can be assessed by monitoring the relative increase of F. This iterative algorithm is known as the variational Bayesian expectation-maximisation algorithm [6]. The result is that the formulas for the sufficient statistics of each unknown distribution can be expressed as a series of equations with mutual dependencies. As discussed in Beal [6] and Cassidy [9], this variational Bayesian EM algorithm reduces to the ordinary frequentist EM algorithm for maximum likelihood estimate [16] if the parameter priors are flat.
A.0.4 Directed acyclic graphs and Markov blanket theory
The model proposed in Section 2 can be viewed as a hierarchical Bayesian model and hence can be represented as a probabilistic directed acyclic graph (DAG). This DAG representation is very useful for visualising the relationships between hidden variables (), parameters () and observations (), each of them being represented as nodes. Typically, square nodes indicate observed variables, round nodes indicate latent random variables and arrows act for conditional dependence. Figure 4 contains the DAG for the full model considered in the present article. Moreover for models having such a structure, the variational Bayesian algorithm benefits from a graphical-related concept hailing from machine learning theory, known as variational message passing [54]. More specifically, the benefits for variational Bayesian models are directly related to the concept of Markov blanket that we will first define.
Definition A.1
The Markov blanket of a node in a DAG () is defined as the set of its parents, , children, , and co-parents, . Two nodes are defined co-parents if they have at least one child node in common [54].
The point of particular interest for variational Bayesian theory is that the variational sequential update equation for a node only depends on expectations over variables in its Markov blanket. It directly follows that equation (15) can be rewritten as
[TABLE]
Fox and Roberts [18] show that equation (16) can be rewritten in an even simpler form as
[TABLE]
Equation (16) is much simpler than equation (15). Similar simplifications can be obtained for equation (14).
A.0.5 Conjugate exponential
The optimal form for and of course depends on the choice of the prior distributions and . The analytical form of can be assessed via the following theorem [6].
Theorem A.2
For models with observed variables , hidden variables and parameters , the mean field variational Bayesian approximation has the following characteristic: if the complete-data likelihood is part of the exponential family (in its “natural form”, meaning parametrized by its natural parameter), and if the hidden and parameter prior distributions and are conjugate to this complete-data likelihood, then the corresponding variational approximate posterior distributions that maximize F, and , are of the same distributional form as the prior distributions and respectively.
As stated in Beal [6, p. 160], with suitable priors, the state space model with unknown parameters is in the conjugate-exponential family, but the “natural form” parametrization required by Theorem A.2 presents a parameter-to-natural parameter mapping that is non-invertible. We can however use the fact that all the nodes of the model defined in Section 2 are conditionally conjugate, implying that the optimal variational approximate posterior distributions and will be of the same distributional form as, respectively, the prior distributions and (a node is said to be conditionally conjugate when its conditional distribution given its Markov blanket (see Definition A.1) is in the same family as its conditional distribution given its parents).
Theorem A.2 moreover ensures that the analytical form of the variational distributions and does not change during iterations. Since this property does not hold for general equations (14) and (15), variational posterior distributions become quickly unmanageable outside the conjugate exponential framework. All proofs of Theorem A.2 and related properties can be found in Beal [6].
A.1 The variational evidence lower bound
In Bayesian analysis, the evidence (or marginal likelihood) provides a natural criterion for model selection by comparing the evidences obtained for the models to be compared. Well-known criteria based on evidence comparison are the Bayesian Information Criterion (BIC) and the Bayes Factor.
As observed in equation (11), the variational Bayesian methodology has the advantage of providing a quantity, the free energy, that is a lower bound for the marginal likelihood of the model and that can be computed efficiently. This leads to a natural criterion for model order selection, which is crucial to estimating our time-varying VAR model, where we have many more variables to estimate than available observations. We can therefore perform model selection by comparing the free-energy quantities computed for each model order and select the model that exhibits the highest . This latter comparison supposes that we have placed uniform priors over each model structure , thereby considering them as equiprobable.
Recalling equations (11), (12) and letting be the model estimated for a specific order , the free-energy quantity can be re-expressed as
[TABLE]
where the first right-hand side term is the average log-likelihood over the entire set of parameters and hidden states that therefore acts as an accuracy term, and the second right-hand side term is the Kullback–Leibler distance between the prior and the variational posterior distributions for the entire set of parameters . Since the Kullback–Leibler distance increases with the number of parameters, this second term acts as a penalty.
The choice of (instead of ) as a criterion for model order selection implicitly assumes that the free energy lies at the same distance to the evidence whatever the model order . This seems reasonable, given that the dataset is the same for all models. A similar procedure for model order selection can be found in Corduneanu and Bishop [14] and Roberts and Penny [42].
It has been shown that in the large sample limit, the free energy becomes equivalent to the Bayesian information criterion (BIC) [2]. This popular model order selection criterion can therefore be seen as a limiting case of the variational Bayesian framework [6, 40].
The specific analytic form of for our model is derived in Appendix B.
Appendix B Computation of the Free Energy
It is straightforward to re-express the free-energy quantity as
[TABLE]
where we omit the conditional dependence to the model for notational simplicity. The first term on the r.h.s. of equation (19) represents the Kullback–Leibler divergence between the prior and the variational posterior for the entire set of parameters , the second r.h.s. term represents the entropy of the variational posterior distribution of the hidden variables and the third r.h.s. term is the average log-likelihood of the data and hidden state parameters taken over the entire set of parameters and hidden state .
Recalling the mean-field factorization assumed in our model as well as the priors distributions set in Section 4, the free-energy quantity takes the explicit form:
[TABLE]
where the conditional dependance of the variational posterior distributions to the data is now omitted for notational simplicity.
The KL divergences of term 1 have closed form for conjugate exponential distributions and are therefore straightforward to obtain despite their intensive computation. Let us focus on the entropy term that appears in term 2. Following [6], this can be rewriten as
[TABLE]
Term 3 therefore disappears in equation (20) and the quantity becomes
[TABLE]
Assuming now that the parameters have a point mass density rather than their variational posterior distribution , the quantity becomes
[TABLE]
The quantity can be obtained just after the forward recursion step and therefore amounts to
[TABLE]
where the quantities and are defined in Appendix D [43].
Appendix C Mean and Fluctuation Theorem
The Mean and Fluctuation Theorem is a decomposition proved by [3]. Using our notation, based on the model and the Conditional distributions defined in Section 2 and the set of unknown parameters , the following decomposition of the density holds
[TABLE]
Note that the formulation in [3] considers the matrix as time-invariant.
Appendix D Unified Inference Theorem
For the BSS model, recall equation (14), which is a key quantity we want to evaluate
[TABLE]
In the situation where are random parameters rather than fixed values, Barber and Chiappa [3] propose an elegant solution based on a suitably augmented system of equations that allows to infer through classical state space model inference algorithms. The theorem states that the above density can be written as
[TABLE]
where the augmented elements are
[TABLE]
and where and are defined as the Cholesky decompositions of
[TABLE]
In our specific case, we have that , due to the non-randomness of the matrix , and so the unique quantity to define is . In the case where and are diagonal as defined in Section 3, we have that
[TABLE]
where is the variational posterior variance of the -th entry of the matrix defined in Section H and is the variational posterior mean of the -th entry of the inverse variance-covariance matrix defined in Section J that can be straightforwardly obtained due to the properties of the inverse-gamma and gamma distributions.
For our model, the augmented system of equations then gives
[TABLE]
where is defined in equation (27). Complete proofs can be found in [3] and [39].
Appendix E E-step: Computation of the Distribution of
As discussed in Section A, variational Bayesian algorithms lead to EM-like iterative equations for the optimal variational densities . The variational update equation form for the hidden variables was state in equation (14). The derivation of equation (14) can be found in Appendix E. Due to the conjugacy condition discussed in Section A and equation the Conditional distributions stated in Section 2, the variational posterior distribution is multivariate Gaussian of dimension . As explained in Beal [6] and Cassidy [9], if the parameters were, as they called, point estimated, equation (14) would be straightforwardly resolved with classical state space model tools.
However in our variational Bayesian scenario, parameters are random variables governed by a specific variational distribution , and so has to be computed for each variable with respect to its variational posterior distribution , sequentially for all parameters in the set . As described in Beal [6] and Barber and Chiappa [3] and fully explained in Ostwald et al. [39], the simplification in equation (14) no longer holds in this random parameter scenario because
[TABLE]
The difference between the two terms in equation (29) is evaluated in Barber and Chiappa [3] and its decomposition is known as the mean and fluctuation theorem, reproduced in Appendix C. To nevertheless use standard state space model algorithms to solve in the variational Bayesian framework, and therefore to capitalize on the vast literature and results that exist on the topic, Barber and Chiappa [3] prove the so-called unified inference theorem recalled in Appendix C. The idea is to reformulate as a standard density with known parameters , by suitably augmenting the first equation in the system defined in Section 2. The classical KRTS smoother algorithms may then be suitably applied to the new density , and allows us to derive the target quantity . This finally corresponds to , and by extension to by equation (14).
The optimal variational posterior distribution for the hidden state sequence is therefore multivariate Gaussian at each time :
[TABLE]
where the sufficient statistics , as well as the cross-moments , are obtained through the KRTS smoother recursive equations applied to the suitable augmented system of equations discussed above.
Forward recursions.
This step implements the recursive equation for . Let , , , and denotes the expectation with respect to the suitable variational distribution . These quantities are obtained recursively as
[TABLE]
where the Kalman gain is given by
[TABLE]
Backward recursions.
The backward recursions implement the recursive equations for . Let now and . These quantities are obtained recursively as
[TABLE]
We also define the cross-time quantities:
[TABLE]
This last element represents .
Appendix F M-step: Computation of the Distribution of
We will derive here the update equation for in the case where is different for each diagonal entry of . Extension to the situation where is a unique parameter related to a unique variance parameter is straightforward.
Recalling the prior form for defined in Section 3 as well as equation (15), the optimal form for the variational posterior becomes:
[TABLE]
where throughout this Supplementary Material, ct will denote the normalization constant for the given density or log-density. Recalling equation , we have that
[TABLE]
and
[TABLE]
The variational posterior for each element is therefore an inverse-gamma distribution with shape and scale parameters defined as:
[TABLE]
Appendix G M-step: Computation of the Distribution of
We will derive here the update equation for in the situation where is diagonal. Extension to the situation where is full is straightforward.
Recalling the prior form for each defined in Section 3 as well as equation (15), the optimal form for the variational posterior becomes:
[TABLE]
where all terms not depending on are put in the constant term ct. Let us consider the two terms separately. Under the prior specified for , term 1 simply equals
[TABLE]
Recalling equation for , term 2 can be rewritten as
[TABLE]
The whole expression for then becomes
[TABLE]
and so the variational posterior for each element is an inverse-gamma distribution with respectively shape and scale parameters defined as:
[TABLE]
Appendix H M-step: Computation of the Distribution of A
The optimal form for the variational posterior is:
[TABLE]
where all terms not depending on are stacked in the constant term ct and the conditional dependence on the data is dropped for sake of brevity. Let us consider the two terms separately. The development of term 2 is the same whatever the form for (full, diagonal or proportional to identity). It becomes
[TABLE]
We can now rewrite (45) in term 2 of equation (44):
[TABLE]
If we suppose the matrix diagonal, as implied by the prior form as well as equation (15), term 1 in equation (44) is straightforward:
[TABLE]
The full expression then gives:
[TABLE]
Therefore
[TABLE]
Appendix I M-step: Computation of the Distribution of
We will derive here the update equation for in the situation where the parameter is unique. Extension to the situation where is related to different is straightforward. Recalling the prior form for in Section 2 as well as equation (15), the optimal form for the variational posterior becomes:
[TABLE]
Recalling the equation for , we have that
[TABLE]
and
[TABLE]
The variational posterior is therefore an inverse-gamma distribution with shape and scale parameters
[TABLE]
Appendix J M-step: Computation of the Distribution of
If the matrix is proportional to identity, i.e. with only one element , the optimal variational form for the posterior becomes:
[TABLE]
where all terms not depending on are included in the constant term. Recalling the equation of , we have that
[TABLE]
and
[TABLE]
The whole expression for can the be rewritten as
[TABLE]
The variational posterior distribution is therefore gamma with shape and scale parameters
[TABLE]
where \Gamma=\sum\limits_{i=1}^{k}\big{[}M_{1_{i,i}}-2M_{3_{i,i}}\psi_{i}\theta_{i}+M_{2_{i,i}}[(\psi_{i}\theta_{i})^{2}+\theta_{i}^{2}]\big{]} and are defined in equation (49).
Appendix K M-step: Computation of the Distribution of
We will derive here the update equation for the parameters . Recalling the prior form for defined in Section 3 and equation (15), the optimal form for the variational posterior becomes:
[TABLE]
Recalling equation for , we have that
[TABLE]
and
[TABLE]
We therefore have
[TABLE]
And then regrouping term 1 and term 2 yield to the inverse-gamma distribution for the variational posterior with shape and scale parameters defined as
[TABLE]
Appendix L M-step: Computation of the Distribution of
Recalling the prior form for , , and equation (15), the optimal form for the variational posterior is
[TABLE]
where all terms not depending on are included in the constant term ct. Let us consider first the term 1. As the prior for defined in Section 3 is inverse-Wishart, term 1 becomes
[TABLE]
Let us now consider term 2:
[TABLE]
Assembling terms 1 and 2 in equations (65) and (66) yields to
[TABLE]
and then
[TABLE]
Appendix M Multiple trials
The model can be modified to deal with conditionally independent sequences which are supposed to have the same hidden state. This reflects the case that arises during an event-related experimental paradigm, where many trials on the same condition are measured.
In Beal [6] and Cassidy [9], this extension is solved by first estimating the necessary sufficient statistics for each sequence independently in the E-step, and then by averaging these statistics to get only one set of sufficient statistics, which is representative of the entire set of independent sequences before performing the M-step. However this approach does not take into account the complex dependence of the variational posterior on the whole dataset . We can however write and solve the full evidence of all the data if we rewrite
[TABLE]
where is block diagonal of dimensions with each diagonal element being identically distributed. The state equation remains the same as in model defined in Section 2. We can observe in model (69) that the time-varying parameter vectors are unique, taking into account the whole dataset . The evidence of the complete model can be rewritten as
[TABLE]
All the computations done so far can easily be adapted to this new setting. The variational posterior distributions relative to model (69) are now conditional on the whole dataset . A detailed derivation of the model and related variational posterior distributions with multiple trials can be founded in Section M.
In the DAG depicted in Figure 4, we can see that conditional dependence on the observed variables is present only for nodes . For the hidden variables, the variational E-step KRTS algorithm can be run on the whole state space system (69) and the resulting variational posterior density will therefore be conditional on the whole dataset . In Appendix E, the update equations remain the same, except for the Kalman gain in equation (32) which becomes
[TABLE]
Looking at the update equations for the parameter in Appendix K, we see that for multiple sequences equation (59) becomes
[TABLE]
where ct denotes the normalization constant for the given density. Therefore equation (62) can be rewritten as
[TABLE]
The resulting shape and scale parameters for therefore become
[TABLE]
Looking at the update equations for the matrix in Appendix L, we see that for multiple sequences, equation (64) becomes
[TABLE]
where each block-diagonal entry is conditional on a unique dataset , and will therefore present the same properties as the unique matrix in Appendix L. The resulting parameters for the variational posterior inverse-Wishart distribution
can therefore be rewritten as
[TABLE]
Appendix N Granger-Causality Detection Results
N.1 Data generated with slowly-varying parameters and normal errors
N.2 Data generated with slowly-varying parameters and non-normal errors
N.3 Windowing estimation procedure
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Arnold et al. [1998] M. Arnold, X. H. R. Milner, H. Witte, R. Bauer, and C. Braun. Adaptive AR modeling of nonstationary time series by means of Kalman filtering. IEEE Transactions on Biomedical Engineering , 45(5):553–562, 1998.
- 2Attias [1999] H. Attias. Inferring Parameters and Structure of Latent Variable Models by Variational Bayes. In Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence , UAI’99, pages 21–30, San Francisco, 1999. Morgan Kaufmann. ISBN 1-55860-614-9. URL http://dl.acm.org/citation.cfm?id=2073796.2073799 .
- 3Barber and Chiappa [2007] D. Barber and S. Chiappa. Unified Inference for Variational Bayesian Linear Gaussian State-Space Models. In B.Schölkopf, J.C.Platt, and T.Hofmann, editors, Advances in Neural Information Processing Systems 19 (NIPS 2006) , pages 81–88. MIT Press, 2007.
- 4Barnett and Seth [2015] L. Barnett and A. K. Seth. Granger causality for state-space models. Physical Review E , 91(4):040101, 2015. doi: 10.1103/Phys Rev E.91.040101 .
- 5Barnett and Seth [2017] L. Barnett and A. K. Seth. Detectability of Granger causality for subsampled continuous-time neurophysiological processes. Journal of Neuroscience Methods , 275:93–121, 2017. ISSN 01650270. doi: 10.1016/j.jneumeth.2016.10.016 .
- 6Beal [2003] M. Beal. Variational Algorithms for Approximate Bayesian Inference . Ph D thesis, Gatsby Computational Neuroscience Unit, University College London, London, 2003.
- 7Box and Tiao [1973] G. E. P. Box and G. C. Tiao. Bayesian Inference in Statistical Analysis . Addison-Wesley, Reading, Massachusetts, 1973.
- 8Casella and George [1992] G. Casella and E. I. George. Explaining the Gibbs sampler. The American Statistician , 46(3):167–174, 1992.