There has been a surge of network models being applied to psychological data sets in recent years. This is consistent with a general call to conceptualize observed psychological processes not merely as indicative of latent common causes but rather as emergent behavior of complex, dynamical systems in which psychological, biological, and sociological components directly affect each other (Borsboom, Cramer, Schmittmann, Epskamp, & Waldorp, 2011; Cramer et al., 2012; Cramer, Waldorp, van der Maas, & Borsboom, 2010; Schmittmann et al., 2013; van der Maas et al., 2006). Such relationships are typically not known, and probabilistic network models (Koller & Friedman, 2009) are used to explore potential dynamical relationships between observables (Epskamp, Maris, Waldorp, & Borsboom, in press; van Borkulo et al., 2014). In this paper, we aim to provide a methodological introduction to a powerful probabilistic network model applicable in exploratory data analysis, the Gaussian graphical model (GGM), and to propose how it can be used and interpreted in the analysis of both cross-sectional and time-series data.

Two lines of network research in psychology. We can currently distinguish two distinct and mostly separate lines of research in which networks are utilized on psychological data sets: the modeling of cross-sectional data and the modeling of intensive repeated measures in relatively short time frames (e.g., several times per day during several weeks). In cross-sectional modeling, a model is applied to a data set in which multiple subjects are measured only once. The most popularly used methods estimate undirected network models—so-called pairwise Markov random fields (Epskamp et al., in press; Murphy, 2012). When the data are continuous and assumed normally distributed, the GGM can be estimated. The GGM estimates a network of partial correlation coefficients—the correlation between two variables after conditioning on all other variables in the data set (Epskamp, Borsboom, & Fried, 2017a). This model is applied extensively to psychological data (e.g., Cramer et al., 2012; Fried, Epskamp, Nesse, Tuerlinckx, & Borsboom, 2016; Isvoranu et al., 2017; Kossakowski et al., 2015; McNally et al., 2015; van Borkulo et al., 2015).

Researchers can obtain time-series data by using the experience sampling method (ESM; Myin-Germeys et al., 2009), in which subjects are asked several times per day to fill out a short questionnaire using a device or smartphone app. Also, time-series data can arise from diary studies (e.g., a questionnaire completed at the end of the day) or physiological measurements, among other methods. Often, repeated measures of one or multiple participants are modeled through the use of (multilevel) vector autoregressive (VAR) models, which estimate how well each variable predicts the measured variables at the next time point (Borsboom & Cramer, 2013). These models are increasingly popular in assessing intraindividual dynamical structures (e.g., Bringmann et al., 2013; Bringmann, Lemmens, Huibers, Borsboom, & Tuerlinckx, 2015; Wigman et al., 2015).

Estimating the GGM is not limited to cross-sectional data; the model merely does not take temporal information into account. As such, the lines of research on network modeling of cross-sectional data and time-series data can naturally be combined. First, GGM models can readily be estimated on repeated measures, if these can be assumed to be temporally independent. Second, as the VAR model can be seen as a generalization of the GGM that takes violations of independence between consecutive cases into account; the GGM can be used to model the residual (innovation) structure of a VAR model to gain insight in the contemporaneous time level of a time-series analysis. Finally, the between-subjects effects of n > 1 studies can also be modeled through the use of the GGM.

Outline.We show that in time-series modeling the GGM allows researchers to extend the modeling framework to incorporate contemporaneous and between-subjects effects. We do this by building up the model complexity in three steps: (1) when cases can be assumed to be independent (e.g., cross-sectional data or repeated measures in which no auto-regression is assumed), (2) temporally ordered data (e.g., n = 1 time-series data or n > 1 time-series data where no individual differences are assumed), and (3) temporally ordered data from multiple subjects (e.g., n > 1 time series). The final level of model complexity leads to a novel contribution of this paper: separation of variance into contemporaneous, temporal, and between-subjects network structures. We propose novel estimation procedures to estimate these models, which we have implemented in two free software packages: mlVAR,¹ and graphicalVAR.² We furthermore expand on existing literature by providing a comprehensive methodological discussion of the GGM, by comparing the GGM to structural equation modeling (SEM; Kaplan, 2000; Wright, 1921), by providing overviews of estimation methods and software packages useable in each kind of data set and by discussing the interpretation of networks estimated at the contemporaneous and between-subjects levels. We showcase network models estimated from n > 1 time-series data in two empirical examples by reanalyzing existing data sets (Bringmann et al., 2013; Geschwind et al., 2011; Mõttus, Epskamp, & Francis, 2017). In the supplementary materials, we provide codes to perform the analyses and we assess the performance of these methods in large-scale simulation studies. To aid the reader in the various different terms used in this paper, we have included a glossary of terms in the Appendix.

1. The Gaussian graphical model

Let denote a random vector with as its realization.³ We assume is centered⁴ and normally distributed with some variance-covariance matrix : (1) The subscript C denotes a case (a row in the spreadsheet). We currently do not define the nature of the observed variables. Thus, can consist of variables relating to one or more subjects, could contain repeated measures on one or more variables, could contain variables of a single subject that do not vary within-subject, and so forth. Consider three examples: (1) Y₁ could represent the level of anxiety of subject p on day 1 and Y₂ the level of anxiety of subject p on day 2, (2) Y₁ could represent the length of subject p and Y₂ the number of times subject p bumps his or her head, and (3) Y₁ could represent the number of cigarettes subject p smokes per day and Y₂ the number of cigarettes another subject p + 1 smokes per day (case C then represents a dyadic pair).

Partial correlation networks. Assuming multivariate normality, encodes all the information necessary to determine how the observed measures relate to one another. However, we will not focus on in this paper but rather on its inverse—the precision matrix : Of particular importance is that the precision matrix can be standardized to encode partial correlation coefficients of two variables, given all other variables (dropping subscript C for notational clarity; Lauritzen 1996):⁵ (2) in which κ_ij denotes an element of , and denotes the set of variables without i and j. These partial correlations can be graphically displayed in a weighted network, in which each variable Y_i is represented as a node, and connections (edges) between these nodes represent the partial correlation between two variables. When the partial correlation (thus the corresponding element in ) equals zero, no edge is drawn. Thus, modeling the inverse variance-covariance matrix, such that every nonzero element is treated as a freely estimated parameter, allows for a sparse model for (i.e., every element in may be nonzero while some elements in are zero; Epskamp, Rhemtulla, & Borsboom 2017d). Such a model is termed a GGM (Lauritzen, 1996). Of note, when the sample-variance-covariance matrix is inverted and standardized, no partial correlation will be exactly equal to zero and the GGM will therefore be saturated. To obtain a sparse model with testable implications, in this paper partial correlations are forced to zero either by using thresholding rules or regularization techniques.

When drawing a GGM as a network (often termed a partial correlation network), positive partial correlations are typically visualized with blue or green edges and negative partial correlations with red edges,⁶ and the absolute strength of a partial correlation is represented by the width and saturation of an edge (Epskamp et al., 2012). When a partial correlation is zero, we draw no edge between two nodes. As such, the GGM can be seen as a network model of conditional associations; no edge indicates that two variables are independent after conditioning on all other variables in the data set. This allows us to model conditional associations, which we might expect to be zero, rather than marginal associations, which we rarely expect to be zero (Meehl, 1990).

To exemplify the above, suppose for three variables “fatigue,” “concentration problems,” and “insomnia,” the true variance-covariance matrix is To model this matrix, we need six parameters (three covariances and three variances). The corresponding true precision matrix becomes Similar to SEM, a model can be devised that perfectly explains this pattern using only five parameters, because one of the elements in can be constrained to be zero (Epskamp et al., 2017d). We can now standardize this matrix and make the off-diagonal elements negative (Equation (2)) to obtain the partial correlation matrix, which we will denote : This matrix can be used to draw a network as is shown in Figure 1. This figure shows that someone who is tired is also more likely to suffer from concentration problems and insomnia. Furthermore, this network shows that the correlation between insomnia and concentration problems can be explained by the relationships of both variables with fatigue: concentration problems and insomnia are conditionally independent given the level of fatigue.

The Gaussian Graphical Model in Cross-Sectional and Time-Series Data

All authors

Sacha Epskamp, Lourens J. Waldorp, René Mõttus & Denny Borsboom

https://doi.org/10.1080/00273171.2018.1454823

Published online:

16 April 2018

Figure 1. A hypothetical example of a GGM on psychological variables. Nodes represent someone’s ability to concentrate, someone’s level of fatigue, and someone’s level of insomnia. Connections between the nodes, termed edges, represent partial correlation coefficients between two variables after conditioning on the third. Blue edges indicate positive partial correlations, red edges indicate negative partial correlations, and the width and saturation of an edge corresponds to the absolute value of the partial correlation.

Display full size

Interpreting GGMs. This paper concerns the exploratory estimation of GGMs from various sources of data, without prior knowledge on the model structure. Such undirected network models can be interpreted in strikingly different ways, ranging from a no causal interpretation to a strong causal interpretation:

In addition, the GGM is closely tied to factor analysis, allowing for extensions to factor modeling through the use of network modeling (Epskamp et al., 2017d). The main focus of this paper is discussing the second interpretation, while also describing how the GGM may be used to show predictive effects. We detail these first two points below by first showing how partial correlation coefficients correspond to multiple regression coefficients and next discussing the relationships between the GGM and SEM. Point 3 follows from observing that the GGM is directly related to similar undirected models such as the Ising model (Ising, 1925). A discussion on the causal interpretation of such models is beyond the scope of this paper, and we refer the reader for this topic to Epskamp et al. (in press) and van Borkulo et al. (2014).

1.1. The Gaussian graphical model and multiple regressions

An edge in a GGM indicates that one node predicts a connected node after controlling for all other nodes in the network. This can also be shown in the relationship between coefficients obtained from least-squares prediction and the inverse variance-covariance matrix. Let represent an k × k matrix with zeros on the diagonal. Furthermore, let represent the ith row of without the ith element (as the diagonal is set to zero), which contains the regression coefficients obtained in a multiple regression model: As such, γ_ij encodes how well the jth variable predicts the ith variable. This predictive effect is naturally symmetric; if knowing someone’s level of insomnia predicts his or her level of fatigue, then conversely knowing someone’s level of fatigue allows us to predict his or her level of insomnia. As a result, γ_ij is proportional to γ_ji. There is a direct relationship between these regression coefficients and the inverse variance-covariance matrix (Meinshausen & Bühlmann, 2006). Let denote a diagonal matrix on which the ith diagonal element is the inverse of the ith residual variance: d_ii = 1/Var(ϵ_Ci). It can then be shown (Pourahmadi, 2011) that⁷ (3) Thus, κ_ij is proportional to both γ_ij and γ_ji; a zero in the inverse variance-covariance matrix indicates that one variable does not predict another. Consequently, the network tells us something about the extent to which variables predict each other. This predictive quality is the cornerstone for how such network models are often applied (Hastie, Tibshirani, & Wainwright, 2015), for example, in recommender-systems that recommend users on products they might like depending on which products the user already liked (Marsman, Waldorp, & Maris, 2017b). In addition to these applications and aiding the interpretation of GGM models, this relationship between multiple regression and undirected network edges plays a crucial role in many network estimation procedures (Haslbeck & Waldorp, 2016b; Meinshausen & Bühlmann, 2006; van Borkulo et al., 2014), including the methods discussed below in this paper.

1.2. The Gaussian graphical model and structural equation modeling

Let represent a set of unobserved variables, which we assume to be jointly normally distributed with . Then, we can form an encompassing framework for several possible generating models:⁸ (4) in which is a diagonal matrix,⁹ indicating that after conditioning on all causes the variables are independent, is a square matrix with zeros on the diagonal of causal effects between observed variables, and is a factor-loading matrix. The variance-covariance matrix of may in turn be modeled in various ways to achieve complicated model setups. The expression above is well-known in SEM, which allows for confirmatory testing of causal models. In exploratory estimation, one could assume no latent variables exist and aim estimate (causal models), or one could assume no relationships between observed variables exist and aim to estimate (factor models). We contrast both to the GGM below.

1.2.2. Factor models

Suppose that, instead of assuming no unobserved causes as in (5), we take the generating model of (4) and only allow for unobserved causes of the observed variables. Then, (4) reduces to the well-known factor model (Brown, 2014). The corresponding model for now becomes: which can subsequently be inverted to obtain an expression for the equivalent GGM. Golino & Epskamp (2017) provide a detailed derivation of this inverted expression and show that a factor in the factor model will lead to its indicators to cluster (all nodes connected to each other with strong edges) in the GGM. This result is in line with mathematical equivalences between factor models and network models of binary variables (Epskamp et al., in press; Kruis & Maris, 2016; Marsman et al., 2017a; Marsman, Maris, Bechger, & Glas, 2015). As there is only one unique inverse to , there is only one unique GGM for every factor model. Conversely, however, one GGM may be equivalent to many different factor models (e.g., all possible rotations of ).

Due to these equivalences, network modeling and factor modeling are closely connected. A natural first step in performing an exploratory factor analysis would be to estimate and draw a GGM model and investigate if the nodes cluster as would be expected by a factor model. Cluster-detection algorithms on the GGM could even be performed to investigate the number of factors to extract (Golino & Epskamp, 2017). Of note, however, is that many GGM estimation methods will always aim to estimate sparse GGM (i.e., contains exact zeroes), which is not expected given a factor model (except when latent variables are orthogonal). As such, estimating a sparse network does not provide evidence that a latent variable model could not have generated the data (Epskamp, Kruis, & Marsman, 2017c; Epskamp et al., in press). GGM modeling can further be used to augment factor analysis by modeling the latent variable variance-covariance matrix or the residual variance-covariance matrix as a GGM (Epskamp et al., 2017d). Modeling as a GGM leads to a latent network model, which can be used in exploratory estimation of relationships between latent variables. Modeling as a GGM leads to a residual network model, which may be used to estimate factor models while local independence is structurally violated (Epskamp et al., 2017d; Pan, Ip, & Dubé, 2017).

A note on spurious edges. When interpreting edges in the GGM as indicative of potential causal effects, it is important to note that edges in a GGM may also result from latent variables. Such edges are termed spurious, and cannot be accounted for unless the latent variable is explicitly modeled (e.g., by using the residual network model described above). The same problem occurs in exploratory DAG estimation, in which case a latent variable may induce a directed edge in the causal network. Furthermore, such spurious associations may arise in any statistical model, to the extent that unmeasured latent variables are involved. Here, the downside that GGM loses information on direction of effect turns into an upside: when an edge is indicative of a causal effect, GGMs do not retrieve the direction of effect, however, when an edge is spurious due to the influence of a latent variable, the GGM also does not introduce a strong causal hypothesis on what would happen under intervention.

2. Estimating GGMs from different sources of data

2.2. Temporally ordered data of a single subject

In line with a call for more intraindividual and person-based research (Molenaar, 2004), an increasingly popular form of data pertains to n = 1 time series, in which a single individual is measured repeatedly over a period of time. One such situation is in clinical practice (Epskamp et al., 2018; Kroeze et al., 2017), where a patient can be measured several times per day over a period of a few weeks. We will limit our discussion to data obtained in a relatively short time-frame so that we can reasonably assume the model will remain stable over time. Then, we can apply the methodology above to obtain a GGM for the n = 1 data. However, such an analysis does not take temporal ordering of data into account (i.e., relationships between measurement occasions) and only investigates contemporaneous relationships between variables (e.g., within the same measurement occasion). This is important for several reasons. First, valuable information, especially in the context of dynamical relationships, might be contained at the temporal level rather than at the contemporaneous level. Second, not taking temporal ordering into account might bias the estimated contemporaneous relationships (see Section 4 of the supplementary materials). For example, if one variable causes itself and another variable at the next time, then not taking temporal ordering into account turns that variable into a latent cause, which would produce an edge in the GGM. Third, temporal information is needed when constructing the joint likelihood over time (e.g., to obtain the information retained in a system over time; Epskamp, 2017a; Quax, Kandhai, & Sloot, 2013). Finally, temporal information can aide in distinguishing reciprocal and cyclic effects by regarding these as acyclic effects unfolding over time.

Vector Auto-regression. The simplest way to deal with temporal ordering of cases is to incorporate the effect between consecutive measurements (Chatfield, 2016; Hamilton, 1994; Shumway & Stoffer, 2010). This is called a Lag-1 model because it includes both measurements at the current time point t as well as measurements from the previous time point t − 1. We will focus our discussion on Lag-1 models, noting that everything below also generalizes to more complicated models (e.g., Lag-2 models). In intraindividual analysis, VAR (Brandt & Williams, 2007; Rosmalen, Wenting, Roest, de Jonge, & Bos, 2012) has gained substantive footing in visualizing temporal information through networks. The lag-1 VAR model can be denoted as a regression model on the previous measurement occasion: (7) The model matrix encodes temporal predictive effects from variables on variables in the next measurement occasion, and can be used to obtain a directed network, which we term the temporal network. The variance-covariance matrix can be inverted () to obtain a GGM modeling effect within the same measurement occasion, after controlling for temporal effects. These can be displayed again as a network, which we term the contemporaneous network.

Temporal networks. Temporal networks, encoded by , have grown popular in recent psychological literature (e.g., Bringmann et al., 2013, 2015; Bos et al., 2017; Klippel et al., 2017; Snippe et al., 2017; Wigman et al., 2015). A temporal network is formed by combining a lagged variable y_{t − 1} and current variable y_t into a single node, connected with directed edges which are weighted according to the regression parameters contained in .¹² Thus, an edge in the temporal network indicates that a node predicts another node (or itself in the common case of self-loops) at the next measurement occasion, after controlling for all other variables at the previous measurement occasion. Temporal prediction is central to the concept of Granger causality in the economic literature (Eichler, 2007; Granger, 1969), and it satisfies at least the temporal requirement for causation (i.e., the cause must precede the effect). Temporal networks may thus highlight potential causal pathways. While temporal networks are typically cyclic, they can also be interpreted as summarizing a DAG unfolding over time.

Contemporaneous networks. In addition to temporal effects, VAR analyses also include contemporaneous effects, which can be modeled as a GGM. We will term this modeling framework (a VAR model with contemporaneous effects explicitly modeled and portrayed as a GGM) graphical VAR (GVAR; Wild et al., 2010).¹³ A useful equivalent way to denote a GVAR model is by using a conditional Gaussian distribution: Which is equivalent to Equation (7). Now, it becomes evident that if consecutive cases can be assumed to be independent, and thus , the GVAR model is exactly the same as the GGM model described above for independent cases. Thus, the GVAR model can be seen as a generalization of the GGM model to temporally ordered data. GVAR only differs from regular VAR in that the contemporaneous structure is modeled and represented as a GGM, instead of being saturated. This leads to a strikingly different interpretation of the VAR model; the VAR model can be seen as an inclusion of temporal effects on a GGM.

Temporal and contemporaneous information. Figure 2 shows a hypothetical example of the two network structures obtained in a GVAR analysis and shows how they might plausibly differ. The left panel shows the temporal network. The self-loop shows that whenever the subject in question felt energetic (or tired), this person also felt more (or less) energetic in the next measurement. The temporal network also shows us that after exercising, this person felt less energetic. The contemporaneous network in the right panel shows a plausible reverse relationship: Whenever this person exercised, he or she felt more energetic in the same measurement occasion. In psychology, there will likely be many causal relationships that occur much faster than the lag interval of a typical ESM study; in this case, these pathways will be captured in the contemporaneous network. For example, if someone is experiencing bodily discomfort, that will immediately negatively affect that person’s ability to enjoy him or herself (Epskamp et al., 2018). Especially when the measurement is on blocks of time (e.g., “since the last measurement did you feel...”), such effects are likely to be caught in the contemporaneous network.

The Gaussian Graphical Model in Cross-Sectional and Time-Series Data

All authors

Sacha Epskamp, Lourens J. Waldorp, René Mõttus & Denny Borsboom

https://doi.org/10.1080/00273171.2018.1454823

Published online:

16 April 2018

Figure 2. A hypothetical example of two network structures obtained from a GVAR analysis. The network on the left indicates the temporal network, demonstrating that a variable predicts another variable at the next time point. The network on the right indicates the contemporaneous network, demonstrating that two variables predict each other at the same time point.

Display full size

2.2.1. Estimation

Estimating saturated (fully connected temporal and contemporaneous networks) GVAR models is straightforward. First, one needs to estimate temporal effects of a regular VAR model by performing multivariate multiple regression of all variables on the previous measurement occasion, or by estimating univariate models for every variable, in which denotes the ith row of . Next one can invert the variance-covariance matrix of the residuals to obtain a GVAR model. Step-wise model selection in latent network models (Epskamp et al., 2017d) can also be used to estimate sparse GVAR models. Missing data can be handled in default ways of SEM or regression models (e.g., listwise deletion and full-information maximum likelihood), or by using more sophisticated techniques such as Bayesian estimation (Schuurman, Grasman, & Hamaker, 2016b) or the Kalman filter (Harvey, 1990; Kim & Nelson, 1999).

Novel estimation methods. A promising recent method for estimating VAR models is the Bayesian dynamical SEM implementation in version 8 of Mplus (Asparouhov, Hamaker, & Muthén, 2017; Muthén & Muthén, 2017), which includes handling of missing data, measurement invariance, and latent variables. Mplus can be used to estimate saturated GVAR models, and to perform model selection in the temporal network of a GVAR model. Model selection in the contemporaneous network of a GVAR model is not yet implemented, but credibility intervals around contemporaneous effects can be obtained by manually inverting each sampled residual variance-covariance matrix (these can be stored using the BPARAMETERS option).

When estimating GVAR models, regularization methods can be used similar to the estimation of GGMs on nontemporally ordered data. Abegaz and Wit (2013) proposed to apply LASSO estimation to jointly estimate the temporal and contemporaneous network structures using the multivariate regression with the covariance estimation (MRCE) algorithm described by Rothman, Levina, & Zhu (2010). MRCE involves iteratively optimizing , using cyclical-coordinate descent, and , using the GLASSO algorithm (Friedman et al., 2008, Friedman et al., 2014). EBIC model selection can be used to obtain the best performing model. This methodology has been implemented in two open source R packages: sparseTSCGM (Abegaz & Wit, 2015), which aims to estimate the model on repeated multivariate genetic data, and graphicalVAR (Epskamp, 2017c), which was designed to estimate the model on the psychological data of a single subject. The graphicalVAR package also allows for unregularized multivariate estimation.

An alternative to estimating GVAR models is to estimate structural VAR (SVAR; Chen et al., 2011) models, also called unified SEM (Gates, Molenaar, Hillary, Ram, & Rovine, 2010). In SVAR, the contemporaneous effects are modeled using a directed network instead of an undirected network. The sparsity of the undirected GVAR contemporaneous network corresponds in the same way to the sparsity of the directed contemporaneous network in SVAR as how the GGM corresponds to causal models (edges arise in the GGM due to edges in the causal network or conditioning on common effects). The temporal SVAR network is sparser than the temporal GVAR network, as contemporaneous mediators can be controlled for in SVAR but not in GVAR. A saturated SVAR model can be obtained by using regressions on the previous time-point as mentioned above, followed by transforming the contemporaneous variance-covariance matrix (e.g., by using a Cholesky decomposition on its inverse; Lütkepohl, 2005) and subsequently transforming the temporal effects to take contemporaneous mediators into account. This technique of obtaining an SVAR model leads to multiple solutions (Beltz & Molenaar, 2016). Step-wise model selection can also be used to estimate sparse SVAR, for example, by using model selection in (unified) SEM (Gates et al., 2010) or Bayesian dynamical SEM models (Asparouhov et al., 2017; Muthén & Muthén, 2017).

2.3. Temporally ordered data of multiple subjects

A type of data that is increasingly common due to the emergence of ESM studies is time series of multiple subjects (e.g., Bringmann et al., 2013, 2015; Mõttus et al., 2017; Schmiedek, Lövdén, & Lindenberger, 2010; Wigman et al., 2015). Such data sets pose a promising gateway to study both intraindividual dynamics and between-subjects overlap as well as their differences. Here, we assume that the number of time points might differ per person and that measurement occasions are nested in people. We can model the temporal data of every person with an individual GVAR model: in which indicates the stationary mean vector of subject p (which enters the model because we can no longer assume within-subject means are zero without loss of generality), encodes the person-specific temporal network, and encodes the person-specific contemporaneous GGM.

Multilevel modeling. To gain insight in the general network structure over subjects, we can investigate the individual networks at a second level. Doing so is termed multilevel modeling, explained in more detail in Section 2.1 of the supplementary materials. Let and encode the expected temporal and contemporaneous network when selecting a person at random. Furthermore, we can assume without loss of generality that data are grand-mean centered. We then obtain: Here, and now encode the average parameters in the population: the fixed effects. Deviations from these fixed effects, such as , are often called random effects. Besides the individual network structures, researchers often aim to estimate the structure and parameters of these fixed effects because these tell us something about the average intraindividual effect. Researchers also aim to estimate the variance-covariance structure of the random effects because it tells us something about individual differences (Bringmann et al., 2013).

The random effects can be modeled by assuming a second level normal distribution on all the parameters. This can be complicated, however, especially when modeling partial correlation coefficients in such a way (e.g., any hierarchical model for needs to take into account that this matrix must remain positive definite). The interpretation of, for example, correlations between different temporal or contemporaneous edges is also difficult. Therefore, we only focus here on a subset of the parameters where we can easily interpret the second-level model: the mean structure. As a result, if a multivariate normal is assumed for all parameters, then it is also assumed for the marginal distribution of the means—regardless of other parameters: Again, we can invert the variance-covariance matrix to obtain a GGM, which we will term the between-subjects network, a network between stationary means of different subjects.¹⁴ As such, estimating the GVAR model on n > 1 time-series analysis allows for the separation of variance into three distinct network structures: temporal networks, contemporaneous networks, and the between-subjects network.

2.3.1. Estimation

In this section, we^{Table 2}outline several different ways in which individual network structures as well as fixed effect network structures may be estimated. We first discuss applying the methodology of estimating n = 1 GVARs discussed above to both pooled data as well as data of each subject separately, followed by a discussion of different multilevel estimation procedures that take clustering of the data into account. An overview of these methods is also included in Table 2.

Pooled and individual LASSO estimation. First, we can estimate a GVAR model for every subject to obtain subject-specific estimates for the temporal and contemporaneous networks. Similarly, we can estimate fixed-effects networks by estimating a GVAR model on the entire within-subjects centered data set, using the sample means of every subject on every variable as a plug-in for the within-subject means. Consequently, we can estimate the between-subjects network by estimating a GGM on the sample means of each subject on all variables. We can readily apply the LASSO regularization methods described earlier for this purpose: the methodology outlined by Abegaz and Wit (2013) to estimate temporal and contemporaneous networks and the methodology outlined by Foygel & Drton (2010) to estimate the between-subjects GGM. We term this framework pooled and individual LASSO estimation and have implemented it in the R package graphicalVAR (Epskamp, 2017c). The performance of pooled and individual LASSO estimation is assessed in simulations reported in Section 3 of the supplementary materials.

Multilevel estimation. The second and third procedures described in Table 2 make use of multilevel modeling (Hamaker, 2012). Two main benefits of this approach are (1) instead of estimating the VAR model in each subject, only the fixed effects and variance-covariance of the random effects need to be estimated, and (2) afterward, estimates of subject-specific parameters can be obtained, which are somewhat pulled together (termed shrinkage). Shrinkage allows the estimation of the model for one subject to borrow information from other subjects. Multilevel estimation can be performed by specifying the multivariate model using hierarchical Bayesian Monte Carlo sampling methods or by integrating over the distribution of the random effects (Gelman & Hill, 2006; Schuurman et al., 2016b).

Multivariate Bayesian multilevel. Bayesian multivariate estimation has proven to be powerful in estimating multivariate multilevel models, especially given its flexibility in adding measurement error, latent variables, and in handling missing data (Schuurman, Houtveen, & Hamaker, 2015). Recently, the dynamic SEM methodology implemented in Mplus version 8 (Asparouhov et al., 2017; Muthén & Muthén, 2017) has made estimation of multivariate multilevel VAR models much faster and more user-friendly than other Bayesian software routines. Specifying a temporal VAR model with correlated random effects is straightforward and relatively fast to compute with a moderate number of variables (e.g., 6). At the time of writing, Mplus does not return partial correlations by default, but these can be obtained by using the BPARAMETERS option and manually inverting the sampled variance-covariance matrices. Mplus allows for specifying random effects on the contemporaneous covariances and thus, by extension, allows for estimating random contemporaneous networks in addition to random temporal networks. Specifying such a model can be done by specifying dummy latent variables for the residual covariance between each pair of variables (a prior guess on the sign of the covariance is needed). Doing so, however, can significantly increase computation length especially when all random effects are allowed to correlate. To facilitate estimation, we have implemented a function generating Mplus code for a multilevel GVAR model and subsequently running the model using the MplusAutomation package (Hallquist & Wiley, 2017) in version 0.4 of the mlVAR package, which can be called using estimator = ''Mplus'' and requires the Mplus program to be installed.

Two-step multilevel VAR. A downside of multivariate estimation is that the number of random effect covariances to be estimated increases quadratically with the number of variables. Forcing random effects to be uncorrelated helps, but places strict assumptions on the model. Bringmann et al. (2013) proposed to estimate multilevel VAR models using univariate models instead, using a frequentist estimation procedure. In this work, the multilevel VAR model is estimated by sequentially estimating univariate multilevel regression models of one variable given all lagged variables. Doing so ignores several correlations of random effects because many parameters are not estimated in the same model, simplifying the analysis: only correlations between incoming edges to the same node and the intercept of that node are included in an univariate model. This method scales up well to approximately eight variables when estimating correlated random effects and around 20 variables when estimating orthogonal random effects (or by using a moving window approach; Bringmann et al., 2015). Of note, when specifying orthogonal random effects not all random effects are assumed to be uncorrelated, merely the ones used in the same univariate model.

The methodology of Bringmann et al. (2013) does not estimate contemporaneous or between-subjects networks. To this end, we extended the algorithm in a framework we term two-step multilevel VAR. The details of this estimation procedure are explained in Section 2 of the supplementary materials. In short, we extend the methodology of Bringmann et al. (2013) by within-subject centering and by adding subject sample means as between-subjects predictors (as discussed by, e.g., Curran & Bauer, 2011; Hoffman & Stawski, 2009; Hamaker & Grasman, 2014). This allows us to estimate between-subjects networks by collecting regression coefficients as in Equation (3) and symmetrizing the resulting matrix.¹⁵ In a second step, we take the residuals of the first analysis and again perform sequential univariate multilevel regression models to predict each residual from all other residuals in the same measurement occasion. Again, these can be collected, as in Equation (3), and symmetrized to obtain contemporaneous networks. Networks can be thresholded by removing all effects that are not significant. For the between-subjects and contemporaneous networks, this results in two p-values for every edge—either both can be required to be significant (“and”-rule) or an edge can be included if one of the two p-values is significant (“or”-rule). Using the “and”-rule means erring more on the side of caution (sparser network), whereas using the “or”-rule means erring more on the side of discovery. We have implemented two-step multilevel VAR in the mlVAR package, which can be called using estimator = ''lmer'' (the default).

Choosing the estimation method. The choice of which estimator to use is not trivial and depends on the interests of the researcher. In Table 2, we list some pros and cons of each of the methodologies. In particular, multilevel estimation can be very complicated and is harder in high-dimensional settings. Assuming normally distributed parameters can also be problematic because doing so imposes that subjects cannot differ on the structure of the networks, merely on their parameterization. When a parameter (e.g., a temporal edge) is zero in some subjects but nonzero in others, then this parameter cannot be normally distributed (the distribution would peak at 0). Therefore, it is currently hard to estimate differently structured individual networks (different edges set to be exactly 0 between subjects) in multilevel estimation. Nonetheless, multilevel estimation particularly shines in that when estimating an individual network, researchers can borrow information from other subjects. We have performed simulation studies to assess the performance of the two proposed methods in this paper: pooled and individual LASSO estimation and two-step multilevel VAR. We report the results of these studies in Section 3 of the supplementary materials, which shows that both methods adequately detect the true fixed-effect network structures with increasing sample size. Having more time points per subject helps to estimate the contemporaneous and temporal networks, and having more subjects helps to estimate the between-subject networks. Two-step multilevel VAR performs well in estimating intraindividual networks when the number of observations is low, but does not perform subject-specific model selection: all estimated intra-individual networks are saturated and contain all edges. Pooled and aggregated LASSO estimation does estimate the structure of intraindividual networks, but performs poorer in intraindividual parameter estimation with fewer observations as no information is borrowed from other subjects.

GIMME. Finally, when analyzing n > 1 data, another option is to estimate SVAR models instead. A promising estimation procedure to estimate such models over many individuals, while dealing with potential heterogeneity, is “group iterative multiple model estimation” (GIMME; Gates & Molenaar, 2012), which is implemented in R using the gimme package (Lane, Gates, Molenaar, Hallquist, & Pike, 2016). In GIMME, no multilevel structure is imposed and subject-specific networks are allowed to differ in structure. Information from other subjects is borrowed, however, in that the structure of individual networks can be based on other subjects (e.g., an edge can be included because it is present in many other subjects). No shrinkage is induced on the parameter estimates that are nonzero (as would be the case in multilevel or hierarchical Bayesian modeling). A variant of GIMME that estimates the GVAR or a combination of structural and GVAR models has not yet been developed, and we note that this may be a promising avenue for further research.

3. Interpreting GGMs estimated from between-subjects data

This paper describes the estimation of network models on data from different subjects (both cross-sectional as well as person-wise average scores). Cross-sectional network modeling is often criticized for inappropriately taking cross-sectional results to be reflective of within-person causal processes (e.g., Bos & Wanders, 2016; Bos et al., 2017), as it can be shown that such results will not equal within-person processes except under strong assumptions (Molenaar, 2004). To this end, this section discusses the interpretation of GGMs estimated from data of different subjects. We first argue that cross-sectional data may be interpreted as a between-subjects analysis, assuming that between-subjects variance is dominant, and next discuss the interpretation of potential causal effects at the between-subjects level.

3.2. Within- and between-subjects effects

In contrast to prior work on multilevel VAR modeling (e.g., Bringmann et al., 2013, 2015; Pe et al., 2015; Wigman et al., 2015), in this paper between-subjects effects are conceptualized in addition to the within-subjects effects in a separate GGM. Furthermore, in contrast to prior work, cross-sectional networks are not interpreted to be reflective of within-subject effects, but rather to potentially reflect a between-subjects structure, assuming that observed scores are not dominated by state-like variance (but see Brose et al., 2013). This raises the question on how such models could be interpreted. In particular, if edges in the GGM are interpreted as generating hypotheses to potential causal pathways, the question is raised how such causal effects can occur at the between-subjects level. This section therefore discusses the topic of causation at the between-subjects level. Here, we interpret the stationary means as being locally stationary: the average of a subject in a relatively short time span of measurement (e.g., a few weeks). As such, we do not interpret the mean vector as a lifetime average. Instead, we assume it could change, potentially due to experimental intervention. As a result, we argue that the between-subjects network can also be indicative of potential causal pathways—regardless of whether it is estimated from a cross-sectional interview concerning variables that are not expected to vary much over time or obtained from estimating the means from time-series data. To simplify the argumentation below, we do not discuss separate temporal and contemporaneous networks but only general within-subjects networks (a GGM of within-subject data without taking temporal ordering into account).^{Figure 3}

Simpson’s paradox. Hamaker (2012) described an example of how within- and between-subject effects can strongly differ from each other. Suppose we let people write several texts, and we measure the number of spelling errors they make and the number of words per minute they type (typing speed). We would expect to see the seemingly paradoxical network structures shown in Figure 3, Panel (a). We would expect a positive relationship in the within-subjects network (e.g., typing faster than your average leads to making more errors). Conversely, we would expect a negative relationship in the between-subject network (e.g., people who type fast, on average, generally make fewer spelling errors). This is because people who type fast, on average, are likely to be more skilled in writing (e.g., a courtroom stenographer) and are less prone to make a lot of spelling errors, compared to someone who types infrequently. Panel (b) of Figure 3 shows another example in which the structures might differ (Hoffman 2015; provided by Hamaker 2017). These network structures show that when people exert more physical activity than their average they likely experience an elevated heart rate, while people who on average are often physically active likely have a lower average heart rate. Such a different effect depending on the level of analysis is well known in the statistical literature as Simpson’s paradox (Simpson, 1951).

The Gaussian Graphical Model in Cross-Sectional and Time-Series Data

All authors

Sacha Epskamp, Lourens J. Waldorp, René Mõttus & Denny Borsboom

https://doi.org/10.1080/00273171.2018.1454823

Published online:

16 April 2018

Figure 3. Two hypothetical examples of differing within- and between-subject networks. The networks on the left indicates the within-subject network, showing that personal deviations from the means predict each other at the same time point, and the networks on the right indicates the between-subjects network, showing how the means of different subjects relate to one another. (a) Example based on Hamaker (2012). (b) Example based on Hoffman (2015); Hamaker (2017).

Display full size

Interventionist accounts of causation. The different ways of thinking about the effects of manipulations in time-series models can be organized in terms recently developed from interventionist accounts of causation (Woodward, 2005). According to Woodward, causation is fleshed out in terms of interventions: X is a cause of Y if an intervention (natural or experimental) on X leads to a change in Y. Statistically, the interventionist account is compatible with, for example, Pearl’s 2000 semantics in terms of a “do-operator.” Here, an intervention on X is represented as Do(X = x), and the causal effect on Y is formally expressed as . Pearl distinguished this from the classical statistical association, in which no intervention is present, and we get the ordinary regression . This notation is useful here, because it can be used to show how different kinds of causal manipulations, each at the intraindividual level, can produce a signal in either the between-subjects or the within-subjects network.

Cashing out causal effects in terms of interventions is useful for understanding the intervention Do(X = x). We can think of this in terms of a random shock to the system, which sets X to value x at a particular time point and evaluates the effect on another variable Y shortly afterward. If we want to gauge this type of causal relationship, we might look at the within-subjects VAR model. Consider Hamaker’s (2012) example regarding typing errors: If a researcher forced a person to type very fast, that researcher would need to evaluate the within-subject data, which would show a positive association between typing speed and the number of errors. In this example, between-subjects data would be misleading because individual differences would probably yield a negative correlation between speed and accuracy—faster typists are more likely to make less errors.

Interventions at the mean level. However, we can also think of a manipulation that sets X to value x in a different way, for instance, by inducing a long-term change in the system that leads it to converge on X = x in expectation. To evaluate the effect of this type of intervention, it is important to consider the behavior of the system as it relates to the changes of the intercept of X. When analyzing time-series data gathered in a relatively short time-span, the within-subjects VAR network as discussed here cannot represent the relevant effects, because it assumes stationarity. However, such effects will be visible in the between-subjects network, which may thus contain important clues to the behavior of the system under potential changes in the intercept of one variable. In terms of Panel (b) of Figure 3, if we are interested in the effect of changing someone structurally—reducing the heart rate of a person on average—our preferred source of hypothesis generation would likely stem from the between-subjects model, as the corresponding within-subject model using the methods described in this paper only models deviations from the stationary mean. Such hypotheses could then be further investigated by using experimental design or lengthier longitudinal data analysis.

Many such examples can be envisioned, especially in the field of psychopathology. For instance, short-term deviations from the mean in abusing a substance might not immediately develop tolerance or lead to one suffering from work or life inferences, but a subject who abuses a substance on average over a long time period might develop these problems (example based on variables used by Rhemtulla et al., 2016). A between-subjects network could similarly show that loneliness mediates the effect of losing a spouse on depressive symptoms (Fried et al., 2015) or highlight the possible effects of childhood trauma and urbanization on psychotic symptoms (Isvoranu et al., 2016, 2017)—both cases in which within-subjects networks based on short-term deviations from the average seem less applicable. This analysis is important because it shows that, even though relevant causal interventions in psychology will typically operate at the intra-individual level, evidence for the effect of such interventions may arise at either the within- or the between-subjects level depending on the nature of the intervention.

4. Empirical examples

4.1. Reanalysis of Mõttus et al. (2017)

We reanalyzed the data of Mõttus et al. (2017) to provide an empirical example of the multilevel VAR methods described above. These data consist of two independent ESM samples, in which items tapping three of the five Five-Factor Model (neuroticism, extraversion, and conscientiousness; McCrae & John, 1992) domains were administered, as was an additional question that asked participants how much they had exercised since the preceding measurement occasion. Sample 1 consisted of 26 people providing 1323 observations in total, and Sample 2 consisted of 62 people providing a total of 2193 observations. Participants in Sample 1 answered questions three times per day, whereas participants in Sample 2 answered questions five times per day. In both samples, the minimum time between measurements was 2 hours. For more information about the samples and the specific questions asked, we refer readers^{Figure 4}to Mõttus et al. (2017).

To obtain an easier and more interpretable example, we first only analyzed questions aimed as measuring the extraversion trait and the question measuring exercise. This leads to five variables of interest: questions pertaining to feeling outgoing, energetic, adventurous, or happy and the question measuring participants’ exercise habits. We analyzed the data using the two-step multilevel VAR procedure as described in detail in the Section 2 of the supplementary materials. We used the mlVAR package, version 0.4, for the estimation of this model. Because the number of variables was small, we estimated the model using correlated temporal and contemporaneous random effects. We ran the model separately for both samples and computed the fixed effects for the temporal, contemporaneous, and between-subjects networks. Correlations of the edge weights indicated that all three networks showed high correspondence between the two samples (temporal network: 0.82, contemporaneous network: 0.94, between-subjects network: 0.70). Owing to the degree of replicability, we combined the two samples and estimated the model on the combined data.

Results. Figure 4 shows the estimated fixed effects of the temporal, contemporaneous, and between-subjects network. In these figures, only significant edges (α = 0.05) are shown. In the contemporaneous and between-subjects networks, an edge was retained if one of the two regressions on which the partial correlation is based was significant (the so-called “or” rule; van Borkulo et al., 2014). These results are in line with the hypothetical example shown in Figure 2: People who exercised were more energetic while exercising and less energetic after exercising. In the between-subjects network, no relationship between exercising,^{Figure 5, Figure 6}feeling energetic, and feeling adventurous was found. The between-subjects network, however, showed a strong relationship between feeling adventurous and exercising: People who, on average, exercised more also felt, on average, more adventurous. This relationship was not present in the temporal network and much weaker in the contemporaneous network. Also noteworthy is that people were less outgoing after exercising. Figure 5 shows the standard deviation of the random effects in the temporal and contemporaneous networks. The largest individual differences in the temporal network were found in the auto-regressions, and the largest individual differences in the contemporaneous network were found in relationship between exercising and feeling energetic.

The Gaussian Graphical Model in Cross-Sectional and Time-Series Data

All authors

Sacha Epskamp, Lourens J. Waldorp, René Mõttus & Denny Borsboom

https://doi.org/10.1080/00273171.2018.1454823

Published online:

16 April 2018

Figure 4. The estimated fixed effects of the three network structures obtainable in multilevel VAR. The model is based on ESM data of 88 people providing a total of 3516 observations. Due to differences in the scale of the networks, the temporal network was drawn with a different maximum value (i.e., the value indicating the strongest edge in the network) than the contemporaneous and between-subjects networks. Edges that were not significantly different from zero were removed from the networks.

Display full size

Figure 4. The estimated fixed effects of the three network structures obtainable in multilevel VAR. The model is based on ESM data of 88 people providing a total of 3516 observations. Due to differences in the scale of the networks, the temporal network was drawn with a different maximum value (i.e., the value indicating the strongest edge in the network) than the contemporaneous and between-subjects networks. Edges that were not significantly different from zero were removed from the networks.

The Gaussian Graphical Model in Cross-Sectional and Time-Series Data

All authors

Sacha Epskamp, Lourens J. Waldorp, René Mõttus & Denny Borsboom

https://doi.org/10.1080/00273171.2018.1454823

Published online:

16 April 2018

Figure 5. The networks showing the standard deviation of random effects in the temporal and contemporaneous networks. Due to scale differences, networks were plotted using different maximum values.

Display full size

Figure 5. The networks showing the standard deviation of random effects in the temporal and contemporaneous networks. Due to scale differences, networks were plotted using different maximum values.

In addition to using only the extraversion and exercise items, we also ran the model on all 17 administered items in the data set. In this analysis, we used orthogonal random effects to estimate the model because correlated random effects cannot be estimated with such a large number of variables. Figure 6 shows the estimated fixed effects of the three network structures; it can be seen that indicators of the three traits tend to cluster together in all three networks. Regarding the node exercise, we found the same relationships between exercise, energetic, and adventurous (also found in the previous example) in the larger networks. Furthermore, we noted that exercising was connected to feeling angry in the between-subjects network but not in the other networks. Finally, there was a between-subjects connection between exercising and feeling self-disciplined: People who, on average, exercised more also felt, on average, more self-disciplined.

The Gaussian Graphical Model in Cross-Sectional and Time-Series Data

All authors

Sacha Epskamp, Lourens J. Waldorp, René Mõttus & Denny Borsboom

https://doi.org/10.1080/00273171.2018.1454823

Published online:

16 April 2018

Figure 6. The estimated fixed effects of the three network structures based on all 17 variables administered. Only significant edges are shown. Legend: 1 = “Worried”; 2 = “Organized”; 3 = “Ambitious”; 4 = “Depressed”; 5 = “Outgoing”; 6 = “Self-Conscious”; 7 = “Self-Disciplined”; 8 = “Energetic”; 9 = “Frustrated”; 10 = “Focused”; 11 = “Guilty”; 12 = “Adventurous”; 13 = “Happy”; 14 = “Control”; 15 = “Achieved”; 16 = “Angry”; 17 = “Exercise.”

Display full size

Figure 6. The estimated fixed effects of the three network structures based on all 17 variables administered. Only significant edges are shown. Legend: 1 = “Worried”; 2 = “Organized”; 3 = “Ambitious”; 4 = “Depressed”; 5 = “Outgoing”; 6 = “Self-Conscious”; 7 = “Self-Disciplined”; 8 = “Energetic”; 9 = “Frustrated”; 10 = “Focused”; 11 = “Guilty”; 12 = “Adventurous”; 13 = “Happy”; 14 = “Control”; 15 = “Achieved”; 16 = “Angry”; 17 = “Exercise.”

4.2. Reanalysis of Bringmann et al. (2013)

To showcase additional information that can be obtained using the GGM model, we reanalyzed the data set used and made publicly available by Bringmann et al. (2013), which has been collected by Geschwind et al. (2011). This data set contains ESM measures of 129 participants, which was collected in two periods over 6 days each: a baseline period and a posttreatment period (mindfulness treatment and a control group). Participants answered 60 measurements per period. Similar to Figure 1 of Bringmann et al. (2013), we analyzed only the baseline data set on the six items selected by Bringmann et al. (2013). We estimated the networks using three modeling frameworks discussed in Table 2. First, we analyze data using multilevel Bayesian estimation using Mplus version 8 (model generated using the mlVAR package). We estimated correlated random effects for the temporal effects but only fixed effects for the contemporaneous effects (making these random led to slow convergence). The model was estimated using three chains that ran until convergence. Nights were handled by adding a row of missing values between consecutive days. Second, we analyzed the data using two-step multilevel VAR estimation as implemented in the mlVAR package, using an “and”-rule and estimating correlated random temporal and contemporaneous effects. Finally, we estimated the data using pooled and individual LASSO estimation using the graphicalVAR package, using γ = 0.25. In the final two analyses, we did not regress the first measurement of the day on the last measurement of the previous day, and removed all pairs of lagged and current variables that contained missing responses. The final sample size was 5927 observations. Edges were retained if they were significant at the α = 0.05 level, or if 0 was not included in the 95% credibility interval.

Results. Figure 7 shows the resulting network structures, and shows that all three methods are mostly aligned. Unsurprisingly, the temporal networks are very similar to those reported by Bringmann et al. (2013).¹⁷ Both the temporal and contemporaneous network are in line with what would be expected under a unidimensional auto-correlated latent variable model (many edges selected, low-rank structure, edges of expected sign) with the exception of the positive temporal edge from “fearful” to “pleasant” in the two-step multilevel network (which was not selected by the other methods). Of note is that Bayesian multilevel estimation resulted in a sparser temporal network. This difference in sparsity is possibly because the multivariate Bayesian multilevel approach more accurately represents the uncertainties in parameter estimation, while two-step multilevel VAR estimates the model piecewise and pooled LASSO estimation does not take the multilevel structure into account. Remarkable is the positive edge between “sad” and “relaxed” in the two-step multilevel between-subjects network, which is based on two significant positive Level 2 regression coefficients (β = 0.202, p = 0.046 and β = 0.151, p = 0.036) where the estimated between-subjects correlation is strongly negative (− 0.53). This edge is especially remarkable since both nodes are strongly connected to other nodes in the network. The Bayesian multilevel between-subjects network showed a similar positive edge between “cheerful” and “worry.” This is noteworthy because under a unidimensional factor model, we would not expect partial correlation coefficients to switch sign from marginal correlation coefficients (Holland & Rosenbaum, 1986; van Bork, Grasman, & Waldorp, 2016). A possible way the partial correlation coefficient switches sign is if it has been conditioned on one or more common effects between the two variables of interest (in this case, potentially “worry,” “pleasant,” or “fearful”). Of course, these effects must be interpreted with great care, especially given the high p-values; we did not control for multiple comparisons, and the same edges are not retained in the other methods. Still, it is noteworthy that if this edge is weak or nonexistent, the between-subjects structure is still not in line with a unidimensional factor model. In such a factor model, “sad” and “relaxed”(which feature the most connections) would be expected to have a strong negative edge between them (a depression factor would lead to “sad” having a strong positive factor loading and “relaxed” having a strong negative factor loading).

The Gaussian Graphical Model in Cross-Sectional and Time-Series Data

All authors

Sacha Epskamp, Lourens J. Waldorp, René Mõttus & Denny Borsboom

https://doi.org/10.1080/00273171.2018.1454823

Published online:

16 April 2018

Figure 7. Reanalysis of the Geschwind et al. (2011) data set used by Bringmann et al. (2013). (a) Fixed effect network structures estimated via multilevel Bayesian estimation. (b) Fixed effect network structures estimated via two-step multilevel estimation. (c) Fixed effect network structures estimated via pooled and individual LASSO estimation.

Display full size

5. Discussion

We discussed the Gaussian graphical model (GGM; Lauritzen, 1996), an undirected network model of partial correlation coefficients, and discussed its utility in the analysis of psychological data sets. The GGM presents a promising exploratory data analysis tool that allows for different levels of interpretation: (1) Edges in the GGM can be interpreted without reliance on a causal interpretation and merely used to show which variables predict each other. (2) Causal effects between variables result in an edge, whereas the lack of a causal effect results in no edge, except in the presence of latent variables or a common effect. The GGM can, therefore, be seen as hypothesis generating structures that highlight potential causal pathways. (3) Undirected models can be used and interpreted as causal data-generating process and have been used as such in several fields of research.

The GGM can readily be estimated on any data set that contains multiple observations of the same variables (e.g., multiple people in cross-sectional data or multiple responses in time-series data). LASSO regularization methods perform especially well in estimating such a GGM structure. In temporally ordered data (e.g., n = 1 time series), the graphical VAR (GVAR; Wild et al., 2010) model generalizes the GGM to incorporate temporal effects. We showed how two network structures can be obtained: a temporal network, which is a directed network of regression coefficients between lagged and current variables, and a contemporaneous network, which is a GGM describing the relationships that remain after controlling for temporal effects. In temporally ordered data of multiple subjects (e.g., n > 1 time series), the natural combination of cross-sectional and time-series data came by adding a third network structure: the between-subjects network, which is a GGM that describes relationships between the stationary means of subjects. We proposed two methods to estimate the three network structures: (1) two-step multilevel estimation, which we implemented in the open source R package mlVAR, and (2) pooled and individual VAR model estimations using LASSO regularization, which we implemented in the open source R package graphicalVAR.

6. Conclusion

This paper provides a methodological overview of how the GGM can be used in various different kinds of psychological data. The GGM can be used to map out unique variance in cross-sectional data or at the contemporaneous and between-subjects levels of time-series analysis. We contrasted this method to exploratory estimation of causal models. While losing information on the direction of effect, estimating GGMs offers an attractive alternative in that these models are uniquely identified, well parameterized, closely related to causal models and also offer exploratory insight on predictive effects between observed variables. When the aim is to discover psychological dynamics, the GGM can be used as a hypothesis generating technique inspiring future research or therapy directions (Epskamp et al., 2018; Kroeze et al., 2017). For example, an effect found in a cross-sectional analysis could inspire a time-series study, a contemporaneous effect could inspire a shorter time-lag time-series study and a between-subjects effect could inspire lengthy longitudinal studies. All network structures may inspire experimental design, or to gather a mixture of observational and experimental data (Magliacane et al., 2017). The GGM thus provides a powerful addition to the exploratory toolbox in behavioral research.

Article Information

Conflict of Interest Disclosures: Each author signed a form for disclosure of potential conflicts of interest. No authors reported any financial or other conflicts of interest in relation to the work described.

Ethical Principles: The authors affirm having followed professional ethical guidelines in preparing this work. These guidelines include obtaining informed consent from human participants, maintaining ethical treatment and respect for the rights of human or animal participants, and ensuring the privacy of participants and their data, such as ensuring that individual participants cannot be identified in reported results or from publicly available original or archival data.

Funding: This work was supported by Grant 406-11-066 from the NWO.

Role of the Funders/Sponsors: None of the funders or sponsors of this research had any role in the design and conduct of the study; collection, management, analysis, and interpretation of data; preparation, review, or approval of the manuscript; or decision to submit the manuscript for publication.

Acknowledgments: The authors would like to thank Laura Bringmann, Noémi Schuurman, Oisín Ryan, and Ellen Hamaker for helpful tips and invigorating discussions, and Katharina Jorgensen for valuable comments on earlier versions of this paper. An earlier version of this paper has been adapted as a chapter in the dissertation of the main author (Epskamp, 2017b).