The SEMs applied in this field are based on the so-called type-III models and are identified by the general expression

where y_t are the values of the ROIs selected, β the parameter estimations, and ζ_t the structural errors associated with each endogenous variable y_t. Each ROI is the result of generating a score (y_t) for each brain volume. This score is calculated through the regression equation established with the values of the voxels which define it by means of a unidimensional principal components analysis. In these terms, we should bear two aspects in mind. Whatever the parameter estimation technique (usually Maximum Likelihood ML), the statistical problem involves estimating the parameters of the β matrix that encompass the effects between ROIs, and the ψ matrix's parameters that encompass the matrix of the variances/covariances between the ζ_t structural errors, so that ψ = E(ζζ′). The specific form that β adopts derives from the aforementioned effects between ROIs, and the form that ψ adopts summarizes the assumptions specified according to the distributions of the structural errors. Generally, the classical assumptions of SEM would involve the initial assumption that E(ζζ′) = E(ytζ′) = 0 and, consequently, the errors should be uncorrelated between themselves in relation to the endogenous variables, except for the possibility that the β matrix considers non-recursive effects (recursive models imply that the connection between two ROIs can only go in one direction, while non-recursive models facilitate reciprocal connections between ROIs). Evidently, this enables the error distribution to be independent of the β estimations. This model also involves the assumption that the variables, i.e., the values for each ROI, are observed continuous variables of a multinormal distribution. The truth is that each value of y_t representing the f MRI of a ROI is estimated through Principal Component Analyses (less often a peak voxel and averaging within sphere approach is also used) according to the selection of a specific number of voxels convoluting under a geometric form (generally a sphere) defined around a voxel of maximum statistical significance, univariate or multivariate, under the statistical assumptions of the massive general linear model. Therefore, every y_t extracted could be considered a latent variable (η_t) defined based on the actual observed values in each voxel and estimated according to the type of design used (Block Design or Event-Related Design) or in some cases, from a series of f MRI data registered in a resting state paradigm.

This issue is rather controversial, and there has been discussion about whether functional or effective connectivity is established through Path Analysis (PA) for observable variables and which does not allow non-recursive effects (reciprocal influence between ROIs), or whether it is strict SEM with latent variables (though not specified as such) and which allows non-recursive effects. To clarify this matter, we should note that the techniques based on PA only admit observable variables, whereas SEM accept both observable and latent variables. Therefore, the final scores of each ROI can be considered as latent. Thus, complex SEM models are essential to the study of connectivity, given that it would be sensible to include the existence of reciprocal effects between ROIs and to estimate what Berry (1984) called complete reciprocity models, that is, specifying all the possible reciprocal effects between variables, or ROIs in this case. This idea is indirectly repeated in the scheme by Inman et al. (2012).

From the point of view of SEM, there are basically two possible structures that can be submitted to an analysis to represent connectivity models: non-recursive and recursive structures, in the latter case either complete or incomplete. Figure 1 shows a simple diagram (with only three ROIs) of both possibilities, with an indication of the structural equations associated with each model.

From this figure, it is inferred that the β matrix can adopt several forms and that both the decomposition of effects and the resulting value of the estimations depend on it. Likewise, the above comment on the E(ζζ′) = 0 assumption yields a ψ diagonal matrix if assumed, or a ψ symmetrical matrix if unassumed, given that E(ζ_iζ_j) = E(ζ_jζ_i). Also, the β matrix can adopt a triangular form in the recursive models or any form in the non-recursive models, as long as the model identification problem is solved by justifying that E(ζζ′) = 0.

These and other matters are not usually dealt with in papers on connectivity and have been solved in different ways. In recent years, many papers have been published in this field applying different statistical approximations to the topic. Sometimes they are not even mentioned, and sometimes the information provided is only partial. The first consequence of this is an extraordinary variety of approaches and statistical treatments, a diversity which occasionally compromises the comparability of results between papers and even a reasonable deduction regarding the functional conception of connectivity. In fact, SEMs present a series of restrictions that are linked to their own statistical structure and that of the associated assumptions which, in most of the fields in which they are applied, are ignored or, at least, no evidence is provided to the contrary. For example, in the field of Social Science and Health Science, in general, there are no data that prove that the model's assumption was complied with: for instance, the distributions of the observed variables, or the order or range conditions, especially in non-recursive models where both conditions could be seriously compromised. This is a situation that can be resolved by the appearance of computer software (Amos, MPlus, or sem library of R), since it incorporates several guarantees. However, little is said about it in model presentations, just as little, or nothing at all, is said of the statistical assumptions. For instance, as shown above, assuming or not assuming that E(ζζ′) = 0 leads us to a very different consideration of the β matrix and of the possible effects to assess. Some reviews have been published on this topic, so we know a certain amount about the limitations of the concept of functional or effective connectivity and its results in relation to its neurobiological parallel. If we look at some papers (James et al., 2009; McCormick et al., 2010; Bianchi et al., 2012; Deshpande and Hu, 2012; Murray et al., 2012; Sawyer et al., 2012; Voineskos et al., 2012; Yang et al., 2012; Bringmann et al., 2013) we find that all of them account for the limits of the concept of connectivity, of its possibilities in relation to effective connectivity, and the derivatives obtained from using it from an applied perspective; but little, if anything, is said about good practices in the use of SEMs and their adaptations for the statistical estimation of brain connectivity.

In view of all of the above, the present paper aims to review the application of SEM for the estimation of connectivity models in work published since 2001. By doing so, we mean to establish the effect of several variables pertaining to studies on connectivity with f MRI signal on the estimation of the R² (Coefficient of Determination representing the proportion of explained variance) each model presents. This way we intend to break through in the systematization of some of the statistical properties in their application and some of the characteristics of SEMs in the field of Computational Neuroscience when generating estimations that lead us to the consideration of a globally-analyzed functioning brain, and when generating models from complexity. At the moment we have no evidence of the possible effect on R² of variables like the type of design used, the number of ROIs defined, the parameter estimation technique used, or the type of sample analyzed, among others. We also aim to offer some recommendations to future users of SEM in this field in order to generate, in the near future, a good mechanism for comparing the results for functional or effective connectivity obtained in different studies.

To be included in the present meta-analysis, the articles had to comply with the following criteria: (a) they had to be original f MRI papers approaching a topic of brain connectivity using a data analysis technique directly related to SEM (so we selected all the papers whose data had been analyzed through SEM, PA, uSEM, and euSEM; (b) each model had to be estimated in a different group or sample, so that even if one paper presented multiple models, each one would be estimated in different samples; (c) they had to have been published between 2001 and 2016; (d) they had to be indexed in Journal Citation Reports, PsycInfo, Pubmed, or Scopus; and (e) they had to explicitly offer β matrix standardized parameter values for each model analyzed and the initial matrix R (Correlation Matrix) or S (Covariance-Variance Matrix) (some authors kindly provided us with this information, since it did not appear in the original published papers). The search for papers was conducted by means of a Boolean algorithm using the following keywords: “f MRI,” “Structural Equation Models,” and “Connectivity.” Studies listed in more than one of the aforementioned sources were not duplicated. With the general selection of the aforementioned keywords, we found 100 papers, 25 of which were rejected due to a lack of sufficient data on basic matters regarding the construction of SEM such as, for example, presenting the values of the specific parameters or not offering multiple equation models, using uni-equation models similar to multiple regression models. Eventually, after two independent reviews of this process, 75 papers were included (available in Data Sheet 1 of the Supplementary Material) containing a total of 160 models to analyze, given that most of them included more than one model (M = 1.47; SD = 1.06), ranging between one and eight models per paper]. The process of including and excluding papers is shown in Figure 2.

To each of the papers selected, we applied a template with which we obtained the values of the different variables assessed in this study by two independent researchers. As has just been mentioned, it should be noted that the majority of papers included more than one structural model, so the data were generated one model at a time; therefore, the number of models is much higher than the number of selected articles. For each model, we registered the variables listed in Tables 1, 2.

As mentioned, in each of the structural models selected, we identified the R²-value. In some cases the papers did not report this value explicitly but included the values of the initial R or S matrix (R is the input correlation matrix and S is the variance-covariance matrix, both as possible inputs data) and the structure and values of the β matrix; in these cases we estimated R² from the available values. In other cases, the contact authors sent us the information mentioned. Following the classical general scheme, and using the classical LISREL notation:

where Y is the matrix of values of specific ROIs, β is the matrix of the effects between ROIs, and ζ is the vector of errors associated with each ROI. It should be noted that the estimations of the β_ij free parameters respond to the so-called path-rule, so that each of the coefficients of correlation between ROIs is decomposed according to the effects established between them in the Structural Model specified. Therefore, the fitting system of a structural model like the one presented here responds to a component of the Σ partitioned matrix that can be represented as follows:

one of the basic expressions in the fit of SEMs, so each r original correlation is compared to the r^* reproduced value derived from Σ. So each residual is estimated by means of (r − r^*), which are evaluated through a χ²-test of fit (depending on the parameter estimation technique used). Accordingly, the residuals can be defined by means of (R-Σ) and therefore

An essential aspect to bear in mind in models of this type is the possibility of all the structural errors being either correlated or uncorrelated. The SEM general scheme requires that if E(ζζ′) ≠ 0, then the specification possibilities of the β matrix become restricted; whereas if we assume that E(ζζ′) = 0 and, therefore, the structural errors are uncorrelated, β can adopt several forms and assume non-recursive effects—which, in the case of functional connectivity, is essential. On the other hand, whatever the structure of the β matrix, it does not consider the specification of free parameters in the principal diagonal, so it is assumed that the β_ii-values are fixed. This is one of the aspects that presents statistical differences with regard to the effective connectivity models based on Dynamic Causal Models. So, if we assume, as is usual in these models, that E(ζζ′) = 0 and, by extension, that E(βζ′) = 0, the above expression can be rearranged so that

thus obtaining the expression of the reproduced matrix Σ as follows:

Therefore, the estimation of R² for each model was obtained by calculating the proportion of variance explained by the following simple calculation, standardizing all the values of the initial var-covar matrix (S)

We consider this to be a robust indicator of the effect between ROIs and therefore available for use as an effect size estimator (Vesterinen et al., 2014). Finally, in the models that do not include the p-value associated with the χ² contrast of fit, this value was reproduced based on the distribution model and the degrees of freedom reported for each model analyzed.

Finally, these values were analyzed following the scheme used by Redondo et al. (2002) adapted to our main objective, with the exception that, in our case, we eventually opted for a random-effect model, given the high variability of the observed distributions of the parameters considered as effect sizes. All the analyses were conducted with the IBM-SPSS software, version 23.0, and with some of the R software routines—more specifically, the Meta library (Schwarzer, 2013) and Mplus version 7.4.

Tables 3, 4 show the values of observed distributions for each variable according to the lists in Tables 1, 2 above.

As can be observed in the above tables, the SEMs studied present some interesting characteristics in accordance with the Carp structure (2012). If we try to define a prototypical SEM model for the study of functional connectivity, it would appear in a paper published between 2009 and 2012 (77 models, representing 48.1% of the total studied), published in the journals NeuroImage or Human Brain Mapping (83 models, i.e., 51.9%), based on the type-III SEM general model (135 models, i.e., 84.4%), extracting the ROIs after a cognitive paradigm based on a block design (Box-Car) or on some of its modifications (128 models, accounting for 80.0%), and which rarely use complete factorial designs (only 30 models, 18.8% of the total). It would be a SEM generated with real samples (only 5% simulate data), and the model specification would as likely come from statistical significances (66 models used a Data-Driven strategy, i.e., 41.3%) as from hypotheses with a neurobiological basis (63 models used that strategy, i.e., 39.4%). The parameter estimations would have been conducted through ML (114 models, 71.3%), and we would have no information on whether the study complies with the condition of application of multinormality for SEM (142 models offer no data on this aspect, i.e., 88.8%). The initial solution in the estimation process would stem from a matrix (S) of second-order centered moments. We would have no data on whether the study complies with the conditions of application of SEM (for example, range or order) given that 142 models offer no data on this (88.8% of the total analyzed). Additionally, from Table 4 we could infer that the general sample would comprise about 25 subjects (± 36.65), which implies a very wide variability in the sample sizes used and which, in mean values, is slightly above Friston's recommendation (Friston, 2012)—around 16 subjects per group analyzed, disregarding some recent clarifications on this topic (Lindquist et al., 2013). Moreover, the specific model would include around 6 or 7 ROIs (±3.499) and between 10 and 11 effects would have been specified (±9.311). Both results indicate a very high variability, which implies dispersion in the formulations and, given that we can expect many more effects specified than the number of ROIs involved, recursive effects would come as no surprise. As regards the fit data, we would have the usual indexes of fit in SEM, like the Goodness of Fit Index (GFI), or the Adjusted Goodness of Fit Index (AGFI), or the Comparative Fit Index (CFI) or the Akaike Criteria (AIC), or the Bayesian Criteria (BIC) among many others (124 models offer this type of information, 77.5%). Likewise, we would obtain some evidence of the global fit since the p-value associated with χ² would be >0.10 (51 models offer this fit, which means 58% of the 88 models with this information and represent 31,9% of the totals of the analyzed models); although it would not be unusual to have SEM models with p-values associated with χ² below the usual criterion, as this is perfectly comparable to the majority of SEM published in other fields. It is unusual to show the percentage of variance explained (R²) by the model (only 2.5% of the models analyzed offer this information) and the χ²-values offered would have a mean value of 50.18 (±162.87), which shows a very high variability and low consistency. From the above, only four of the total of 160 models include the explicit value of R² associated with each model; in the other 156 models the value of R² was estimated from the procedure discussed above.

Finally, in this case, the estimates of R² in each model can be considered independent of the sample size used since the sample estimates of R (S) and that of Σ do not depend on this value. In the same way, we can assume that the individual estimate of R² is equal to the common weighted estimate since the estimation procedure was applied equally in all the models, thus guaranteeing the homogeneity of R² in all the models studied.

For this section, we decided to use the estimations derived from the fit of a General Linear Model, as is usual in classical meta-analyses, in order to establish the differential effects between each of the categories of the moderating variables used, or to estimate the correlations between the distributions of the quantitative moderating variables. In both cases, the corresponding parameters were estimated through ANOVA for the categorical variables and through Linear Regression for the quantitative variables.

As is well-known, estimations of effect size in any meta-analysis are subject to several perverse effects which may generate bias in the process: in this case, the estimation of R² in each of the models which did not report it (i.e., 154 of the 160 studied). Some of the usual problems are related to the different number of variables (ROIs) involved in each model and dissimilar sample sizes. In line with Cheung (2015), we approached these two issues based on the conceptions of the SEMs for the study of meta-analyses (Meta-Analysis Structural Equation Models, MASEM).

In this case, we applied the strategy described by Cheung (2015), using SEMs to estimate the homogeneity of effect sizes (Q index) and the percentage of variation attributable to the models analyzed (I² index). To estimate both indicators, we used the following expressions:

where σ^ẽi2= ∑i=1k(ỹi-wi β^F)2k , which involves the estimation of the error variance ẽ_i in the k models, using (k) instead of (k − 1) for a better adjustment to the demands of the maximum likelihood estimations (ML) and following a χ² distribution. Along with this premise, the estimation of I² was conducted with the simple expression

Both values showed the heterogeneity of the effect sizes (Q² = 1032.46; df = 160; p < 0.001) and the variation can be explained by the effects between models rather than by intra variability (I² = 0.846). In the light of these results we estimated the effects corresponding to each moderator variable by means of Mplus (MASEM). Tables 5, 6 summarize the information from both analyses with specification of the β_ij parameters (impact of each variable on R²) for each variable or category, depending on the case.

When studying the effect of the moderating variables on the value of R² for each structural model analyzed, the first point to note is the larger size of the effect generally observed in the different values of the qualitative variables (Table 5). Additionally, the type of design with the greatest effect is Box Car, whether it is one group (β = 0.824; p < 0.001, η² = 0.462) or more than one group (β = 1.100; p < 0.001, η² = 0.482). The Resting State design category is also statistically significant, although with a very low effect (β = 0.589; p = 0.004, η² = 0.073). The models in which the tasks are compared present the greatest effect and their intensity is high (β = 0.907; p < 0.001, η² = 0.524). As regards the type of population studied, in studies based on healthy subjects that present the greatest effect (β = 0.844; p = < 0.001, η² = 0.488), although the effect is similar to that found in studies based on both normal and clinical populations (β = 1.050; p = < 0.001, η² = 0.432). The models that do not account for recursive effects present the greatest effect, and with a high-intensity effect size in this case (β = 0.942; p < 0.001, η² = 0.508). Finally, as regards the p-value associated with the chi-square value, the models presenting a p > 0.10 yield a high-intensity effect size measurement (β = 0.961; p < 0.001, η² = 0.549).

With regard to the quantitative variables of the value of R², the significant effect is the number of paths established in the model (β = 0.021; p = 0.001, r² = 0.106). Though a low-intensity effect, it is important to remember that the association with the R² of the SEM model is positive; this suggests that a complex effect could occur according to which, with more complex models, we would obtain some more explained variance, i.e., higher in statistical terms, but not necessarily yielding a higher number of significant parameters. This effect, widely described in other areas, is equally important here with regard to the idea of complexity associated with a possible network of functional connectivity.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.