FIAR: An R Package for Analyzing Functional Integration in the Brain

Functional integration in the brain refers to distributed interactions among functionally segregated regions. Investigation of effective connectivity in brain networks, i.e, the directed causal influence that one brain region exerts over another region, is being increasingly recognized as an important tool for understanding brain function in neuroimaging studies. Methods for identifying intrinsic relationships among elements in a network are increasingly in demand. Over the last few decades several techniques such as Bayesian networks, Granger causality, and dynamic causal models have been developed to identify causal relations in dynamic systems. At the same time, established techniques such as structural equation modeling (SEM) are being modified and extended in order to reveal underlying interactions in imaging data. In the R package FIAR, which stands for Functional Integration Analysis in R, we have implemented many of the latest techniques for analyzing brain networks based on functional magnetic resonance imaging (fMRI) data. The package can be used to analyze experimental data, but also to simulate data under certain models.


Introduction
One of the most intruiging problems in simultaneous recordings of two or more signals from the nervous system is the detection of information flow such as causal relations and the timing between them. This principle of brain organization is known as functional integration (Friston et al. 2005). Functional integration refers to distributed interactions among functionally segregated regions. Studies of functional integration try to understand how regional responses are mediated by connections between brain areas and how these connections change with experimental manipulations or disease.
In the case of functional magnetic resonance imaging (fMRI) data this is especially challenging, since the neural data are blurred with the hemodynamic signal. The result is that the signal of interest is only measured in an indirect way. Furthermore, the relation between the neural signal and the measured hemodynamic signal differs from region to region and even from subject to subject (Aguirre et al. 1998;Handwerker et al. 2004), so the hemodynamic response function can only be estimated. Nevertheless, distinguishing between efferent and afferent connections in brain networks is crucial to construct formal theories of brain function. As a consequence, many statistical methods to study brain connectivity based on hemodynamic measurements have been developed.
The three most widely used methods to study functional integration are (1) dynamic causal modeling (DCM), (2) structural equation modeling (SEM), and (3) Granger causality (GC). In the past years, DCM (which is part of the Statistical Parametric Mapping software, SPM, Ashburner et al. 2008) has become a "gold standard" for studying effective connectivity between brain regions, i.e., the direct influence one brain region has over another (Friston 2009). DCM employs an explicit forward model for explaining which (neural) states caused the (hemodynamic) data. DCM assumes that hemodynamic signals are caused by changes in local neural activity, mediated by experimental inputs (e.g., the presentation of a visual stimulus or the instruction to attend to motion) and the distributed neural interactions among brain regions. DCM is based on a model of this distributed processing and is parameterized by the strength of coupling among the neural regions. This neural model is then supplemented with a hemodynamic model that converts the neural activity into predicted hemodynamic signals. The convolution or impulse response function, mapping from underlying neural activity to observed fMRI responses, is called a hemodynamic response function (HRF, Buxton et al. 1998). DCM runs in the MATLAB (The MathWorks, Inc. 2010) environment and to our knowledge, there is no package for R (R Development Core Team 2011) that can perform DCM.
Perhaps the most widely used method to study effective connectivity in the brain is SEM. Although developed in the fields of econometry and the social sciences, it has been succesfully introduced into the neurosciences to study causal pathways in the brain (McIntosh and Gonzales-Lima 1994). Since its introduction, extensions of the classical SEM framework have been developed to specifically analyze time series data. For example, Kim and colleagues proposed a unified structural equation model (USEM) approach for modeling brain connectivity (Kim et al. 2007). USEM unifies a vector autoregressive (VAR) model , represented by longitudinal pathways, and a conventional SEM, represented by contemporaneous pathways. Consequently, USEM is able to model the autoregressive nature within each time series and the correlations between the d-dimensional time series simultaneously.
USEM may be performed with standard SEM software such as Mplus (Muthèn and Muthèn 2004), LISREL (Jöreskog and Sörbom 2005), EQS (Bentler 1995), or some packages in R like sem (Fox 2010) or lavaan (Rosseel 2011). However, the autoregressive connectivity model and data matrix need to be specified manually. For networks with a large number of brain regions and a high autoregressive order, manually specifying the model is difficult. An R function extending the model and data matrix automatically to the desired autoregressive order and consecutively computing the model fit is not available.
A third very popular tool to study brain connectivity is Granger causality (Granger 1969). The idea of GC is that the causal influence of one time series on another can be conceived by the notion that the prediction of one time series is improved by incorporating knowledge about the other. Only very recently the standard framework for GC has been extended to the multivariate case, where predictor and dependent variables are no longer constrained to be univariate (Barrett et al. 2010).
There are several R packages like lmtest (Zeileis and Hothorn 2002) and vars (Pfaff 2008) that provide a GC test. Unfortunately, only bivariate relations may be examined where one time series is caused by one single other time series. However, when studying brain networks usually more than two regions are considered at the same time and repeated bivariate analysis may be misleading. For example, one time series may falsely appear to cause another if they are both influenced by a third time series but with a different delay. Therefore, it would be useful to be able to perform multivariate GC on fMRI time series. These multivariate GC techniques have not been implemented into an R package yet.
These facts led us to believe that there is a need for an assembled package to perform some of the most popular and recent techniques for studying functional integration in brain networks. We call this package FIAR, which stands for functional integration analysis in R. The package is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/ package=FIAR and currently covers DCM, (unified) SEM, (multivariate) GC. We plan to keep it up to date and even extend it with other techniques in the future. In the next section the implemented techniques are briefly discussed and the corresponding R functions demonstrated. In the final section, all techniques are applied to the attention to visual motion dataset that may be downloaded from the SPM website (http://www.fil.ion.ucl.ac.uk/spm/data/).

Theoretical background
Whereas SEM and GC were developed in other areas of science, DCM has been specifically designed for the analysis of functional imaging data. In fMRI, the blood oxygen level-dependent (BOLD) signal is observed, while the signal of interest is the (hidden) neural activity in each region. DCM therefore estimates two models simultaneously: a causal model, indicating which neural activity in one region causes changes in activity in other regions and a hemodynamic model of how the fMRI signals were produced by complicated physiological events, initiated by changes in neural activity.
Overall, DCM models the temporal change in the neural activityż as a bilinear function of the current state z, the inputs u (usually the experimental design) and some neural coupling parameters A, B, and C:ż where A represents the anatomical connections between brain regions . These connections can be seen as average connections between regions, irrespective of the inputs. The B connections are modulated by the inputs j (modulatory or functional connections) and can be added to the A connections in order to obtain the total strength of a connection under input j. Finally, C connections are the direct inputs to the nodes of the system. These parameters can be expressed as partial derivatives: DCM combines this neural model with a plausible and experimentally validated hemodynamic model (the so called Balloon model, Buxton et al. 1998;Friston et al. 2000) with six more parameters that describes the transformation from neural activity to BOLD activity (Stephan et al. 2007). The hemodynamic model is described in full detail in Friston et al. (2000).
Combining the neural and hemodynamic states in a joint state vector x, and their respective parameters into a joint parameter vector θ, results in the state equatioṅ which can be integrated and passed through the output nonlinearity λ to predict the BOLD signal y y = λ(x).
This so called forward model is the basis for estimating the neural and hemodynamic parameters from the measured data. The term "causal" in DCM also refers to this forward model, because it models how hidden neural changes are causing the observed data. A fully Bayesian approach is used to estimate the neurodynamic and hemodynamic parameters . Experimentally validated priors are used for the hemodynamic parameters and conservative shrinking parameters for the neural coupling parameters, embedded in an expectation maximization (EM) algorithm.
Although DCM is probably the most sophisticated of the three models discussed here, we would like to point out two drawbacks. First, due to the large number of free parameters DCM is computationally demanding and this restricts analysis to a small number of regions. Also, the deterministic input output nature of the model only allows to model variations explained by the inputs. This forces all dynamic information to be represented in the design matrix. In other words, it assumes that all neural dynamics can be captured without error from the chosen inputs. This assumption makes DCM estimation very dependent on the exact number and form of the exogenous inputs. Recently, stochastic extensions of DCM have been presented (Daunizeau et al. 2011), but these are not yet widely used and are currently not implemented in FIAR.
Second, the specification of the prior neural and hemodynamic model will put restrictions on the information that can be captured by them. For instance, in DCM the hemodynamic model has much more affordance for delayed coherent variations than the neurodynamic model. Therefore, the delay will be put into the hemodynamics in the fitting of the model. Not because this is actually the case in the data at hand, but because the model has assumed this to be true (Roebroeck et al. 2011).

Model specification
A DCM analysis begins with the specification of the model. In FIAR the model may be specified in three ways. One way is to manually create a DCM list containing all necessary   or to manually enter all parameters when auto = FALSE (default). For a full overview how to specify a DCM, we refer to Appendix A where an example is presented.

Data generation
FIAR allows to generate data under a specified model with a desired signal to noise ratio (SNR) and autoregressive (AR) coefficient. The function dcmGenerate will create the simulated time series with length v (rows) of n regions (columns). When SNR is set to 0, the pure  signal is returned. When SNR > 0 the signal is mixed with Gaussian white noise to achieve the specified SNR when ar = 0. When ar > 0, the noise will be autocorrelated with a coefficient ar. When SNR = 0, the argument ar has no function. Finally, names is a string vector that allows to give names to the regions. As an example we take the connectivity model in Figure 1.
The model contains three brain regions V1, V2, and V3. There is an anatomical connection between region V1 and V2 with a strength of 0.7 Hertz and an anatomical connection between region V2 and V3 of 0.4 Hertz. The first experimental input, which has an onset every 60 scans for a duration of 30 scans, directly influences region V1 with a strength of 0.4 Hertz and it also influences (modulatory influence) the connection between region V2 and V3 with a strength of 0.2 Hertz. The second experimental input, with an onset every 30 scans and a duration of 15 scans, influences region V2 with a strength of 0.3 Hertz and creates a functional pathway from region V3 to region V2 with a strength of 0.2 Hertz. The DCM object in Appendix A contains all necessary scanner and model parameters to specify the model. We can generate data from it as follows: R> set.seed(11111112) R> ts <-dcmGenerate(DCM, SNR = 1, ar = 0.2, names = c("V1", "V2", "V3")) R> head ( This will produce three time series V1, V2, and V3 that are integrated as specified in DCM$a, DCM$b, and DCM$c, with a SNR = 1 and autocorrelation coefficient of 0.2.

Model estimation
As already mentioned, DCM is a Bayesian method, meaning that posterior model parameters are estimated based on prior information and the data. The function dcmEstimate takes as arguments the object DCM, which contains the model, and ts, which represents the time series. The estimation process results in posterior values of the model parameters and their probabilities. The object ts may be the simulated data or experimental data extracted from brain regions.

R> DCM <-dcmEstimate(DCM, ts = ts)
EM-step(-): The posterior parameters are denoted by capital letters in analogy to the small letters of the priors. For example DCM$A are the posterior anatomical connections. They are represented by an n × n matrix where the colums represent the "from" region and the rows the "to" region. The unit of the connections is Hertz. For example, the following output

V1
V2 V3 V1 -1.0000000 0.0000000 0 V2 0.4789123 -1.0000000 0 V3 0.0000000 0.3011896 -1 means that there is an inhibitory anatomical connection from region V1 to itself with a strength of 1 Hertz, an exhibitory anatomical connection from region V1 to region V2 with a strength of 0.48 Hertz, and so on. Only the off diagonal elements are important to interpret. The fields where we did not expect connections in our prior model are the fields that contain zero.
The probability that these connections differ from zero can be found in the DCM$pA field that is created during the estimation process. Figure 2: Posterior connections in the example model.

V1
V2 V3 V1 0.000000 0.0000000 0 V2 0.999976 0.0000000 0 V3 0.000000 0.9992914 0 We have very strong evidence that the connections from V1 to V2 and from V2 to V3 are different from zero. The posterior connections after estimation are relatively close to the true parameters, as can be seen in Figure 2.
Based on the posteriors, we can compute the model evidence with the function dcmEvidence. This will add the fields DCM$AIC and DCM$BIC to the DCM object which produce the Akaike information criterion (Akaike 1973) and Bayesian information criterion (Schwartz 1978), respectively. For example: We can do this for multiple DCMs and compare the fitvalues with the function R> dcmCompare(DCM1, DCM2) where DCM1 and DCM2 represent the different estimated DCMs that are being compared. The returned values of this latter function are Bayes factors (BF, Raftery 1995). This BF can be used to choose one model over another (see example below).

Theoretical background
Structural equation modeling (SEM) was developed in the field of econometrics and the social sciences and first applied to neuroimaging data by McIntosh and Gonzales-Lima (1994). They comprise a set of regions with causal relations between them. Causal relations are thus not computed by the data but are assumed a priori. The strengths of the a priori connections are tuned in such a way that they minimize the discrepancy between the observed and expected instantaneous correlations between regions.
There are two major drawbacks to SEM for fMRI analysis. The first is that the observed correlations in the data are treated as if they reflect causal relations between neural activity. SEM rests upon a phenomenological model of dependencies among the (hemodynamic) data, without reference to how the data were caused at the neural level. The second drawback is that SEMs do not make use of temporal information. If the time indices of the time series were randomly permuted, this would not change the correlations between time series. In other words, SEM does not take into account the dynamics of the time series, but rather computes static correlations.
In the last years, extensions of this classical SEM framework have been developed to better fit the specific nature of fMRI data. One approach that is implemented in FIAR is unified SEM (Kim et al. 2007). USEM unifies a vector autoregressive (VAR) model ) and a conventional SEM. Consequently, USEM is able to model the autoregressive nature within each time series and the correlations between the d-dimensional time series simultaneously.
Specifically, longitudinal temporal relations are defined as relationships between brain regions involving different time points, and are represented in the form of a multivariate autoregressive model of order p, MAR(p). Conversely, contemporaneous relations reflect relationships between brain regions at the same time point, and involve conventional SEMs. Let y j (t) be the jth variable (e.g., the average BOLD intensity for the jth ROI) measured at time t, j = 1, 2, . . . , m. The m-dimensional MAR(p) with an added component of contemporaneous relations can be written as: Here is an (m×1) vector of white noise with zero-mean and error covariance θ ; A is the parameter matrix of the contemporaneous relations, and Φ(i), i = 1, . . . , p is a series of (m×m) parameter matrices representing the longitudinal relations. The diagonal elements in the Φ(i) represent the coefficients of the autoregressive process for each variable, and the off-diagonal elements represent the coefficients of the lagged relation between the variables.
If we denote φ as the set of free parameters of Φ, A, and θ, than these parameters are estimated so that the implied covariance matrix Σ(φ) is close enough to the sample covariance matrix S by minimizing a discrepancy function F (S, Σ(φ)). The most widely used fitting function for general structural equation models is the maximum likelihood function defined by: with q the number of variables in the model (Kim et al. 2007).

Software implementation
The package FIAR contains the function ARsem which is a wrapper around the sem function from the R package lavaan (Rosseel 2011). The function takes as first argument the connectivity model of a classical SEM analysis. The model is specified as a vector that takes 1 when we assume a connection between 2 regions and 0 otherwise. The columns represent the "from" regions and the rows the "to" regions. The argument data should only contain the time series (rows) of the regions (colums) in the model. It is important that the time series are given a name. The third argument order represents the autoregressive (AR) order of the connectivity model. The function automatically transforms both model and data to the extended model of the specified AR order and the lagged dataset that it needs respectively. Notice that setting order = 0 (default) is equivalent to performing a classical SEM analysis. Consider for example a three node network with a connection from region X to Y and from region Y to Z: SEM can be used to investigate causal relations between regions X, Y , and Z, but causality has a different meaning than in DCM. Causality in (7) is assumed a priori and tested as the presence of instantaneous correlations between regions. Since correlations are bidirectional, this SEM will produce the same model fit as the symmetrical model shown in (8).
The causality tested in (7) is different than in (8), because it is assumed to be different a priori. Nevertheless, due to their symmetry the models will produce identically the same fit and can not be compared. Autoregressive SEM allows to compare symmetrical models and their different causal structures. In a first order USEM, Y at time t (Y t ) is no longer solely influenced by X at time t (X t ), but also by the previous state of X (X t−1 ) and its own previous state (Y t−1 ).
The arrow from X t to Y t is a contemporaneous arrow and reflects the covariance between the two ROIs. The contemporaneous arrows have the same meaning as in the classical approach. The extension is in the longitudinal arrows (e.g., X t−1 →X t ), which try to model the autocorrelations within each region in a direct way. Notice that this model is no longer symmetrical to the AR model of order 1 we would obtain if we extend the model from (8). Hence, we can compare the models and test which causality structure is more likely given the data. Fitting the classical SEM from (7) in FIAR is performed via R> model0 <-c(0, 0, 0, 1, 0, 0, 0, 1, 0) R> fit0 <-ARsem(model0, data = semdata) R> summary (fit0) Lavaan ( The example data semdata contains three variables. Based on the model, we test whether there is a connection from X to Y and from Y to Z. The function ARsem returns a fit object of class lavaan for which several methods are available, including a summary method. The standard output produces the model fit which is not good in this case (χ 2 1 = 1340, p < 0.001) and the estimates for the connections (X → Y = 0.12, Y → Z = −0.15).
All new variables are of the form originalname_order and represent the lagged variables up to the specified order, with variables of order 0 being the original variables. To determine which AR order is suitable for the dataset at hand, we can use the function ARorder() by typing: R> order <-ARorder(semdata) R> order The function ARorder fits all AR models from order min to order max and returns the order which produces the lowest AIC (Akaike 1973). The returned model order balances the variance accounted for, against the number of coefficients to be estimated. In our example, it appears that an AR(3) model best fits the data. We can test this USEM by typing: R> fit3 <-ARsem(model0, data = semdata, order = 3) R> summary (fit3) The summary function produces the model fit and the estimates of all connections. In Appendix B we show the summary of the AR(3) model fit. It immediately shows that extending a simple model with only 3 nodes and 2 connections to an AR(3) model results in a lot of extra nodes (9) and connections (36). However, using the function ARsem there is no need to manually construct the lagged dataset nor the extended AR model.
We now also have two types of regression coefficients in the model. Next to the contemporaneous correlations (e.g., from X 0 to Y 0 ) as was the case for the classical SEM, the model also estimates autoregressive relations (e.g., X 1 to X 0 ). These estimates express the autocorrelation within one brain region up to the specified order (in this case 3).
If we want to compare the classical SEM and the AR(3) SEM in terms of model fit, we can compute both AICs:

R> AIC(fit3)
[1] 142437.3 and we see that the likelihood of the classical model (fit0) in this case is higher than that of the AR model (fit3).

Granger causality
4.1. Theoretical background Wiener (1956) proposed a way to measure the causal influence of one time series on another by conceiving the notion that the prediction of one time series could be improved by incorporating knowledge about the other. Granger (1969) formalized this notion in the context of linear vector autoregression (VAR) modeling of stochastic processes (Guo et al. 2008). Although Granger causality (GC, sometimes called Wiener-Granger causality) was developed in the field of econometrics, it has recently received a lot of attention in the neuroscience community (Roebroeck et al. 2005;Guo et al. 2008;Deshpande et al. 2010;Bressler and Seth 2011).
Suppose X, Y , Z are univariate time series of length t. The joint autoregressive representation of Y and Z can be written as Checking whether X Granger causes Y given Z we can extend the model via and the Granger causality from X to Y conditioned on Z may be expressed as As can be seen, the returned value is a type of F -statistic. If the information in the previous time points of X helps predicting the current time point in Y , we expect VAR( 2t ) to be smaller than VAR( 1t ), resulting in a F 2 value larger than zero. If X does not help predicting Y , VAR( 2t ) and VAR( 1t ) will be of the same magnitude, producing F 2 statistics around zero.
Only very recently, this GC measure has been defined for multivariate time series. If the dependent variable Y t is no longer univariate, but a vector of observed responses measured at time t (Y t = Y 1t , Y 2t , . . . , Y pt ) then (12) is no longer useful. There have been two attempts to extend this formula to the multivariate case. The first was proposed by Ladroue et al.
and uses the multivariate mean squared error or trace of the error covariance matrix, tr[Σ( t )], leading to: Another approach, originally proposed by Geweke (1982), is to use the generalized variance or determinant of the error covariance matrix, |Σ( t )|, which leads to: The advantage of (14) over (13) is that it asymptotically approximates a χ 2 distribution for large samples (Barrett et al. 2010). For small samples however, the distribution is unknown and resampling techniques are more appropriate to construct the null distribution.
When all variables are Gaussian, this solution is also fully equivalent to the transfer entropy T x→y|z , an information-theoretic notion of causality (Barnett et al. 2009). For more properties and advantages of solution (14) over (13), we refer to Barrett et al. (2010). Conditional multivatiate GC (CMGC) is computed in FIAR based on (14) and is implemented as condGranger.
Next to conditional MGC, FIAR allows to compute partial MGC (PMGC), introduced by Guo et al. (2008). PMGC differs from CMGC in the inclusion of the present conditioning variables Z in the joint autoregressive representation of Y and Z Checking whether X partial Granger causes Y , given Z leads to and the partial Granger causality from X to Y conditioned on Z may be expressed as The F 1 statistic is again a ratio of two error variance/covariance matrices. If the logratio is significantly larger than zero, we conclude that the variables X partial Granger cause the variables Y , conditioned on the variables Z. As was the case for CMGC, the H 0 distribution of no causality is unknown. Therefore, resampling techniques must be used to construct confidence intervals around the F 1 estimator.
CMGC already controls for latent and/or exogenous influences to some extent, because the determinant of the error covariance matrix is sensitive to residual correlations. However, PMGC takes into account even more correlations, specifically to minimize the influence of sources of variation outside the model. When such influences are expected to be strong and uniform over all variables in the model, PMGC is to be preferred over CMGC (Barrett et al. 2010). PMGC is implemented in FIAR as partGranger. Notice that, when no conditional variables are included in the model, PMGC and CMGC will produce the same MGC measure.
Yet another type of GC test was presented in the work of Roebroeck et al. (2005). Geweke (1982) proposed a measure of linear dependence F xy between x and y which implements GC in terms of VAR models. F xy is the sum of three components where F x→y is the directed influence from x to y, F y→x is the directed influence from y to x, and F x.y is the instantaneous influence. Roebroeck et al. (2005) used the difference between F x→y and F y→x as a GC measure for the influence region x has on region y. A drawback of this approach is that feedback connections are not modeled. On the other hand, studies have shown that, with sufficiently short repetition time (TR), this measure better controls for the loss of information that arises from the low-pass filtering introduced by the HRF than conditional GC (Roebroeck et al. 2005). This difference measure is implemented in FIAR as diffGranger and its partial counterpart (based on PMGC) as pdiffGranger.
GC is yet another type of causality than the DCM and SEM notion of causality. GC is based on temporal precedence in the sense that previous observations in one time series should help predicting the current observation in a second time series. Furthermore, this additive information should reach statistical significance before we have evidence that time series one Granger causes time series two (Bressler and Seth 2011).
The biggest difference between MGC and the other two methods is that MGC can be used in an exploratory way. There is no need to construct an a priori connectivity model. In the case there is very little known about the connectivity model under investigation, this can be seen as an advantage over DCM and SEM.
The drawbacks of this method for use with fMRI data are similar as is the case for SEM. That is, there is no (hemo)dynamic model involved, so correlations are computed between BOLD signals as if they were neural signals. Also, in their current implementation, USEM and GC are restricted to modeling linear relations, which may seem a strong assumption when modeling brain connectivity. However, nonlinear systems often have extensive linear regimes, and when dealing with large-scale interactions linear approximations are found to work extremely well (Bressler and Seth 2011).

Software implementation
The function condGranger computes the conditional GC of a set of predictor variables X (one or more) on a set of dependent variables Y (one or more), conditional on a third set of variables Z (zero or more).
There is a logical argument boot in the function controlling the bootstrap procedure. When boot = FALSE, the GC function merely returns the F -like statistic. When boot = TRUE, the F -statistic is returned together with a bootstrap approximation of F and a bootstrap bias and standard deviation (SD). This allows to make inference on the F statistic based on the bootstrap H 0 distribution. When the original F -statistic falls within ±2SD around the bootstrap mean the newly created field fit$sig will take the value 0, and 1 otherwise.
The bootstrap is performed by the tsboot function from the package boot (Canty and Ripley 2011). The block length of the resampling chuncks of the time series is optimized by the b.star function from the package np (Hayfield and Racine 2008). The argument bs is by default set to 100 and represents the number of bootstrap samples that are taken.
It is very important how the data are entered in the function. The first nx columns (default nx = 1) are the predictor variables, followed by the ny (default ny = 1) dependent variables, and finally the variables we want to condition on. It is important to respect this order. After the data matrix has been constructed correctly, the number of predictor (nx) and dependent (ny) variables need to be specified. R> condGranger(grangerdata, nx = 1, ny = 2, order = 3) [1] 0.6457208 This tests whether variable x multivariate Granger causes variables y and z, conditioned on variables q and w. When boot = TRUE, not only the original F -value is returned, but also a bootstrap bias and standard error. This allows to construct bootstrap confidence intervals.

STATIONARY BOOTSTRAP FOR TIME SERIES Average Block Length of 70
Call: tsboot(tseries = data, statistic = condGranger, R = bs, l = l, sim = "geom", nx = nx, ny = ny, order = order) Bootstrap Statistics : original bias std. error t1* 0.6457208 -0.4498475 0.04592638 This bootstrap results in an F 2 statistic that is 0.65 − 0.45 = 0.2 with a standard deviation of 0.046. The probability that the original statistic of 0.65 comes out of this null distribution is therefore very small. Hence, we reject the null hypothesis of no effect and conclude that variable x Granger causes variables y and z, conditional on variables q and w. This is also reflected in the value of the created field fit$sig: Calculating partial MGC is very analogous. Again, the function can be used with or without the bootstrap. For example: R> partGranger(grangerdata, nx = 1, ny = 2, order = 3) [1] 0.8596885 R> set.seed(33333334) R> fit <-partGranger(grangerdata, nx = 1, ny = 2, order = 3, boot = TRUE) R> fit

Average Block Length of 70
Call: tsboot(tseries = data, statistic = partGranger, R = bs, l = l, sim = "geom", nx = nx, ny = ny, order = order) Bootstrap Statistics : original bias std. error t1* 0.8596885 -0.6554322 0.05589294 The bootstrap tells us that it is very unlikely that the F 1 statistic of 0.86 comes from the null distribution, so we conclude that variable x partial Granger causes variables y and z, conditional on variables q and w. Finally, diffGranger computes the difference between the MGC measures F x→y|z and F y→x|z . Its partial counterpart is implemented as pdiffGranger. Both functions can be used with the bootstrap option.

Average Block Length of 70
Call: tsboot(tseries = data, statistic = diffGranger, R = bs, l = l, sim = "geom", nx = nx, ny = ny, order = order) Bootstrap Statistics : original bias std. error t1* 0.6432665 -0.46078 0.05272068 The interpretation of this measure is somewhat different than for the other two types of GC. The measure of 0.643 is a relative measure and expresses the difference between the CMGC from variable X to variables Y and Z minus the CMGC from variables Y and Z to X. Since we know that the CMGC measure for the first causal path is 0.645, this implies that there is almost no CMGC in the other direction. Hence, X Granger causes Y and Z far more than the other way around.

Application to attention to visual motion data set
In this section we demonstrate the package FIAR on the example data set "Attention to visual motion" from the SPM website (http://www.fil.ion.ucl.ac.uk/spm/data/). This data was obtained by Buchel and Friston (1997). The experiment consisted of four conditions: (i) fixation (F), (ii) static (S), non-moving dots), (iii) no attention (N, moving dots but no attention required), and (iv) attention (A). The GLM analyses showed that activity in area V5 was not only enhanced by moving stimuli, but also by attention to motion. This effect in V5 was modeled using a DCM. Details about the experiment and the design parameters can be found in the SPM8 manual (Ashburner et al. 2008).

R> ts <-attentiondata
The last field (DCM$X0) need to be imported from the SPM analysis. This field represents the filtered and whitened design matrix and can be found in the field xY.X0 from the extracted volumes of interest (VOI) in SPM. The field ts contains the data from the extracted VOIs and can be found in the fields Y from the extracted VOIs in SPM. Every column of ts represents the time series in one region of the network and the order should correspond to the order in which the model parameters DCM$a, DCM$b, and DCM$c are specified. Here the order is V1, V5, SPC. After specifying the model, we fit the data to the model: For example, we find a connection from region V1 to region V5 of 0.15 Hertz and the probability that this is significantly different from zero is 0.97.
In order to test how likely this model is compared to other possible models, we can define a second model, say DCM2, which has no feedback connection between V5 and V1, and then compare the fit of the two models: R> DCM2 <-DCM R> DCM2$a <-c(0, 0, 0, 1, 0, 1, 0, 1, 0) R> DCM2 <-dcmEstimate(DCM2, ts = ts) EM-step(-): The Bayes factor (BF) expresses the evidence for one model over another. In general, a BF larger than 20 is seen as strong evidence for the first model entered, while a BF lower than 0.05 is seen as strong evidence for the second model entered (Raftery 1995). In our example, we find a very large BF, so we conclude that the model with the feedback loop from V1 to V5 (DCM) is to be preferred over the other model (DCM2). This does not imply that the model with the connection from V5 to V1 is the "true" model. It merely implies that the data suggest a connection, rather than no connection, given the priors of the DCM forward model. In practice, many models are theoretically possible and could be compared to each other using the BF. One way to find further evidence for this model is to test it with a different statistical method.

USEM analysis
If we want to test this connectivity model using USEM, we have to specify it in a different way.
USEM is a stochastic model, so we can not incorporate the input connections and modulatory connections directly in the model. Only the anatomical connections can be specified. In order to estimate the above mentioned model with USEM, we have to specify the connectivity matrix as follows: R> Model1 <-c(0, 1, 0, 1, 0, 1, 0, 1, 0) where the rows and colums correspond to region V1, V5 and SPC respectively.
The time series need to be transformed as well. When testing correlations between brain regions, only those time points corresponding to the experimental condition should be extracted (Honey et al. 2002). The function SEMextract has been provided for this purpose. The function takes as arguments ts, the time series, ons, the onsets of the experimental condition to be tested, and dur, the duration of the activation condition of interest. For convenience, we use the vector of onsets and duration specified earlier in the DCM object: R> ts_photic <-SEMextract(DCM$attentiondata, ons = DCM$ons$photic, + dur = DCM$dur$photic) We start with fitting the data to a classical SEM: R> fit1 <-ARsem(Model1, data = ts_photic) However, when we fit the model, lavaan produces an error. This is to be expected given the data at hand. The three time series form a 3×3 covariance matrix with six data points (3 · (3 + 1)/2). However, in the model, we need to estimate four anatomical connections and three error variances (all regions are endogenous), leading to seven free parameters for only six data points. This means this specific model is underdetermined (i.e., negative degrees of freedom) and can not be estimated with SEM.
One way to resolve this is to specify the feedback connections as one connection. Remember that SEM computes covariances between time series. If there is a feedback loop between region V1 and region V5, this should be reflected in a covariance between both regions. The same holds for the feedback loop between V5 and SPC. Therefore, we specify the following model: R> Model2 <-c(0, 0, 0, 1, 0, 0, 0, 1, 0) In this model, we only need to estimate two connections and two error variances (for V5 and for SPC) with six data points. Again, we fit the data to the model: The model fit is satisfactory (χ 2 1 = 2.3, p = 0.13) and we see strong correlations between region V1 and V5 (0.70) and between region V5 and SPC (0.62). The information on the direction of activation is lost (no feedback loops), but the correlations between regions under influence of visual presentation are preserved in the model. If we want to extend this classical SEM in order to estimate autoregressive connections, we can type: R> fit3 <-ARsem(Model2, data = ts_photic, order = 1) R> summary(fit3) This fits a SEM of order AR(1). We see that the model fit is worse than in the classical SEM (χ 2 6 = 34, p < 0.05). This is also reflected in a much higher AIC for the AR(1) model than for the classical model:

R> AIC(fit3)
[1] 4099.306 Although this simplified model is not identical to the feedback model tested with DCM, we again find evidence for connections between V1 and V5, and V5 and SPC respectively.

Granger analysis
Finally, we can investigate the functional integration between the three regions using GC. Since GC is an exploratory method, we do not need to specify an a priori connectivity model. It suffice to test the connections one at a time. For example, testing whether region V1 Granger causes region V5, conditioned on region SPC is done like this: We find no evidence that region V1 Granger causes region V5. Just as was the case for USEM, computation of GC requires specification of the model order. Too low order can lead to a poor representation of the data, whereas too high order can lead to problems of model estimation. In fMRI data this order is usually very high. The exact AR order of the data may be computed using ARorder. Notice that in all GC analyses of the real data set we used the default model order of one. We can compute that the actual order of the data is 89. However, testing GC for model order 89 leads to a singular solution and can not be tested. In order to test if we can find a connection from region V1 to region V5 for a much higher model order, we tested a model of order 50. Even with model order 50, we find no evidence that region V1 Granger causes region V5.

Overview of results
DCM, SEM, and GC were applied to a real-world application. The example data set "Attention to visual motion" from the SPM website was analyzed with the three methods as a practical example. We focused on the feedback loop between V1 and V5.
The DCM results showed strong evidence for a model with feedback loops between V1 and V5, and V5 and SPC. Comparison between this model and a competing model with no connection between V5 and V1, confirmed the evidence for the feedback model. This does not imply that the feedback model is the "true" model. It merely implies that the data suggest such connections, given the assumptions of DCM. Other biologically plausible models could be tested and compared to each other in the same way.
The feedback model could not be tested with SEM, due to negative degrees of freedom. Simplifying the model by removing the feedback connections lead to similar conclusions as with the DCM analysis. Regions V1, V5, and SPC are exchanging information when visual stimuli are presented.
Finally, we tested the relations between the three regions with conditional, partial and difference GC. The bootstrap results of all three GC methods indicated that there is no evidence that region V1 Granger causes region V5, irrespective of the model order used. This is opposite to the DCM and SEM results.
These results show what often will happen when analyzing the same data set with different methods. In this example, a DCM user will find strong evidence for the existence of a causal path from V1 to V5 under presentation of a visual stimulus. She will even find evidence for a feedback loop going from V5 to V1 based on a model comparison. The GC user will find no evidence that V1 Granger causes V5 with a model of order 1. He might start testing more models with higher orders, but would still not find evidence that V1 adds information to the prediction of V5, given past information of V5. The exclusive SEM user would have to solve the problem of an underdetermined model, which means testing a simplified and thus different covariance structure.
That these different methods lead to different results is partly due to the different types of causality that are being measured. We believe that testing these different notions of causality and finding converging evidence between them should become standard practice when studying causal networks in the brain.

General conclusion
We presented the package FIAR which allows one to perform functional integration with many of the latest techniques such as unified SEM, multivariate GC, and DCM. All of these techniques are widely used to investigate functional brain connectivity, but none were implemented in R to date.
In the paper we discussed some similarities and differences between the methods. DCM is probably the most sophisticated of the three, but is computationally demanding. In its present implementation in FIAR it is a deterministic model, assuming all dynamics are captured without error by the design matrix. USEM and GC are computationally simpler and the parameters are easier to interpret. However, the models do not take into account the hemodynamics of the data, which may lead to false inference.
Another issue with USEM and GC is that a model order needs to be specified. This model order is a trade off between variance accounted for and the number of parameters to be estimated. FIAR allows to calculate the order that optimizes this trade off.
We implemented three types of GC in FIAR. GC differs from DCM and USEM in that the latter are confirmatory methods, testing an a priori causal network, while the former can be used in an exploratory way. All types of GC can be used in a univariate or multivariate way and with or without conditioning variables. As such, the univariate conditional GC with only two time series reduces to the well established bivariate GC.
All methods were applied to a real-world data set. The results illustrate that different methods for studying functional integration can lead to completely different conclusions. This is partly due to the different types of causality that are being measured. We believe that testing these different notions of causality and finding converging evidence between them should become standard practice when studying causal networks in the brain. The brain is a very complex system that is neither linear, bilinear nor deterministic; neither bivariate, nor predictable. The abstractions and choices to be made in useful models of brain connectivity are therefore unlikely to be accommodated by one single "master" model (Roebroeck et al. 2011). Since "Essentially, all models are wrong, but some are useful" (Box and Draper 1987), we hope we have brought a selection of useful models to R with the package FIAR.