Multivariate Behavioral Research

ISSN: 0027-3171 (Print) 1532-7906 (Online) Journal homepage: [Link]

The Gaussian Graphical Model in Cross-Sectional

and Time-Series Data

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

To cite this article: Sacha Epskamp, Lourens J. Waldorp, René Mõttus & Denny Borsboom (2018)
The Gaussian Graphical Model in Cross-Sectional and Time-Series Data, Multivariate Behavioral
Research, 53:4, 453-480, DOI: 10.1080/00273171.2018.1454823

To link to this article: [Link]

© 2018 Taylor & Francis Group, LLC

View supplementary material

Published online: 16 Apr 2018.

Submit your article to this journal

Article views: 2556

View Crossmark data

Citing articles: 3 View citing articles

Full Terms & Conditions of access and use can be found at

[Link]
MULTIVARIATE BEHAVIORAL RESEARCH, 2018
VOL. 53, NO. 4, 453–480
[Link]

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

Sacha Epskampa , Lourens J. Waldorpa , René Mõttusb , and Denny Borsbooma
a
Department of Psychological Methods, University of Amsterdam; b Department of Psychology, University of Edinburgh

ABSTRACT KEYWORDS
We discuss the Gaussian graphical model (GGM; an undirected network of partial correlation coeffi- Time-series analysis;
cients) and detail its utility as an exploratory data analysis tool. The GGM shows which variables predict multilevel modeling;
one-another, allows for sparse modeling of covariance structures, and may highlight potential causal multivariate analysis;
exploratory-data analysis;
relationships between observed variables. We describe the utility in three kinds of psychological data
network modeling
sets: data sets in which consecutive cases are assumed independent (e.g., cross-sectional data), tem-
porally ordered data sets (e.g., n = 1time series), and a mixture of the 2 (e.g., n > 1time series). In time-
series analysis, the GGM can be used to model the residual structure of a vector-autoregression anal-
ysis (VAR), also termed graphical VAR. Two network models can then be obtained: a temporal network
and a contemporaneous network. When analyzing data from multiple subjects, a GGM can also be
formed on the covariance structure of stationary means—the between-subjects network. We discuss
the interpretation of these models and propose estimation methods to obtain these networks, which
we implement in the R packages graphicalVAR and mlVAR. The methods are showcased in two empir-
ical examples, and simulation studies on these methods are included in the supplementary materials.

There has been a surge of network models being applied data and the modeling of intensive repeated measures
to psychological data sets in recent years. This is consis- in relatively short time frames (e.g., several times per
tent with a general call to conceptualize observed psy- day during several weeks). In cross-sectional model-
chological processes not merely as indicative of latent ing, a model is applied to a data set in which multiple
common causes but rather as emergent behavior of com- subjects are measured only once. The most popularly
plex, dynamical systems in which psychological, biologi- used methods estimate undirected network models—so-
cal, and sociological components directly affect each other called pairwise Markov random fields (Epskamp et al.,
(Borsboom, Cramer, Schmittmann, Epskamp, & Wal- in press; Murphy, 2012). When the data are continu-
dorp, 2011; Cramer et al., 2012; Cramer, Waldorp, van ous and assumed normally distributed, the GGM can
der Maas, & Borsboom, 2010; Schmittmann et al., 2013; be estimated. The GGM estimates a network of partial
van der Maas et al., 2006). Such relationships are typically correlation coefficients—the correlation between two
not known, and probabilistic network models (Koller & variables after conditioning on all other variables in the
Friedman, 2009) are used to explore potential dynamical data set (Epskamp, Borsboom, & Fried, 2017a). This
relationships between observables (Epskamp, Maris, Wal- model is applied extensively to psychological data (e.g.,
dorp, & Borsboom, in press; van Borkulo et al., 2014). In Cramer et al., 2012; Fried, Epskamp, Nesse, Tuerlinckx, &
this paper, we aim to provide a methodological introduc- Borsboom, 2016; Isvoranu et al., 2017; Kossakowski et al.,
tion to a powerful probabilistic network model applicable 2015; McNally et al., 2015; van Borkulo et al., 2015).
in exploratory data analysis, the Gaussian graphical model Researchers can obtain time-series data by using the
(GGM), and to propose how it can be used and inter- experience sampling method (ESM; Myin-Germeys et al.,
preted in the analysis of both cross-sectional and time- 2009), in which subjects are asked several times per day
series data. to fill out a short questionnaire using a device or smart-
Two lines of network research in psychology. We phone app. Also, time-series data can arise from diary
can currently distinguish two distinct and mostly sepa- studies (e.g., a questionnaire completed at the end of
rate lines of research in which networks are utilized on the day) or physiological measurements, among other
psychological data sets: the modeling of cross-sectional methods. Often, repeated measures of one or multiple

CONTACT Sacha Epskamp [Link]@[Link] Department of Psychological Methods, University of Amsterdam,  WX Amsterdam, Nether-
lands.
Color versions of one or more of the figures in this article can be found online at [Link]/hmbr.
Supplemental data for this article can be accessed on the publisher’s website.
©  Taylor & Francis Group, LLC
454 S. EPSKAMP ET AL.

participants are modeled through the use of (multilevel) (Bringmann et al., 2013; Geschwind et al., 2011; Mõttus,
vector autoregressive (VAR) models, which estimate how Epskamp, & Francis, 2017). In the supplementary mate-
well each variable predicts the measured variables at the rials, we provide codes to perform the analyses and we
next time point (Borsboom & Cramer, 2013). These mod- assess the performance of these methods in large-scale
els are increasingly popular in assessing intraindividual simulation studies. To aid the reader in the various differ-
dynamical structures (e.g., Bringmann et al., 2013; Bring- ent terms used in this paper, we have included a glossary
mann, Lemmens, Huibers, Borsboom, & Tuerlinckx, of terms in the Appendix.
2015; Wigman et al., 2015).
Estimating the GGM is not limited to cross-sectional
data; the model merely does not take temporal informa- 1. The Gaussian graphical model
tion into account. As such, the lines of research on net- Let yC = [YC1 YC2 . . . YCm ] denote a random vector with
work modeling of cross-sectional data and time-series y c as its realization.3 We assume yC is centered4 and nor-
data can naturally be combined. First, GGM models can mally distributed with some variance-covariance matrix
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 yC ∼ N(00,  ). (1)
that takes violations of independence between consecu-
tive cases into account; the GGM can be used to model The subscript C denotes a case (a row in the spreadsheet).
the residual (innovation) structure of a VAR model to gain We currently do not define the nature of the observed
insight in the contemporaneous time level of a time-series variables. Thus, yC can consist of variables relating to one
analysis. Finally, the between-subjects effects of n > 1 or more subjects, could contain repeated measures on one
studies can also be modeled through the use of the GGM. or more variables, could contain variables of a single sub-
[Link] show that in time-series modeling the ject that do not vary within-subject, and so forth. Con-
GGM allows researchers to extend the modeling frame- sider three examples: (1) Y1 could represent the level of
work to incorporate contemporaneous and between- anxiety of subject p on day 1 and Y2 the level of anxiety
subjects effects. We do this by building up the model of subject p on day 2, (2) Y1 could represent the length
complexity in three steps: (1) when cases can be assumed of subject p and Y2 the number of times subject p bumps
to be independent (e.g., cross-sectional data or repeated his or her head, and (3) Y1 could represent the number of
measures in which no auto-regression is assumed), (2) cigarettes subject p smokes per day and Y2 the number of
temporally ordered data (e.g., n = 1 time-series data or cigarettes another subject p + 1 smokes per day (case C
n > 1 time-series data where no individual differences then represents a dyadic pair).
are assumed), and (3) temporally ordered data from Partial correlation networks. Assuming multivari-
multiple subjects (e.g., n > 1 time series). The final level ate normality,  encodes all the information necessary
of model complexity leads to a novel contribution of to determine how the observed measures relate to one
this paper: separation of variance into contemporane- another. However, we will not focus on  in this paper
ous, temporal, and between-subjects network structures. but rather on its inverse—the precision matrix K :
We propose novel estimation procedures to estimate
these models, which we have implemented in two free K =  −1 .
software packages: mlVAR,1 and graphicalVAR.2 We
furthermore expand on existing literature by providing a Of particular importance is that the precision matrix can
comprehensive methodological discussion of the GGM, be standardized to encode partial correlation coefficients
by comparing the GGM to structural equation model- of two variables, given all other variables (dropping sub-
ing (SEM; Kaplan, 2000; Wright, 1921), by providing script C for notational clarity; Lauritzen 1996):5
overviews of estimation methods and software packages   κi j
useable in each kind of data set and by discussing the Cor Yi , Y j | y −(i, j) = − √ √ , (2)
κii κ j j
interpretation of networks estimated at the contempora-
neous and between-subjects levels. We showcase network  We use capitalized subscripts to denote random variables and lower case
models estimated from n > 1 time-series data in two subscripts to denote fixed variables. A variable can potentially be fixed with
empirical examples by reanalyzing existing data sets respect to one subscript but random with respect to another. Supplementary
materials’Section  contains a complete overview of the notation used in this
paper.
 CRAN link: [Link] Github link (develop-  Because we assume data to be centered, we do not need to model the (grand)
mental): [Link] mean vector. This simplifies notation.
 CRAN link: [Link] Github link  This relationship can be traced back much further. For example, Heiser ()
(developmental): [Link] traced this relationship back to the work of Guttman ().
 MULTIVARIATE BEHAVIORAL RESEARCH 455

in which κi j denotes an element of K , and y −(i, j) denotes

the set of variables without i and j. These partial corre-
Fatigue
lations can be graphically displayed in a weighted net-
work, in which each variable Yi 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 K ) −0.25 0.3
equals zero, no edge is drawn. Thus, modeling the inverse
variance-covariance matrix, such that every nonzero ele-
ment is treated as a freely estimated parameter, allows for Concentration Insomnia
a sparse model for  (i.e., every element in  may be
nonzero while some elements in K are zero; Epskamp,
Rhemtulla, & Borsboom 2017d). Such a model is termed Figure . A hypothetical example of a GGM on psychological
a GGM (Lauritzen, 1996). Of note, when the sample- variables. Nodes represent someone’s ability to concentrate, some-
variance-covariance matrix is inverted and standardized, one’s level of fatigue, and someone’s level of insomnia. Connec-
tions between the nodes, termed edges, represent partial correla-
no partial correlation will be exactly equal to zero and the tion coefficients between two variables after conditioning on the
GGM will therefore be saturated. To obtain a sparse model third. Blue edges indicate positive partial correlations, red edges
with testable implications, in this paper partial correla- indicate negative partial correlations, and the width and satura-
tions are forced to zero either by using thresholding rules tion of an edge corresponds to the absolute value of the partial
or regularization techniques. correlation.
When drawing a GGM as a network (often termed a
precision matrix becomes
partial correlation network), positive partial correlations
are typically visualized with blue or green edges and neg- ⎡ ⎤
1.18 0.28 −0.34
ative partial correlations with red edges,6 and the abso- K = = −1 ⎣ 0.28 1.07 0 ⎦.
lute strength of a partial correlation is represented by −0.34 0 1.11
the width and saturation of an edge (Epskamp et al.,
2012). When a partial correlation is zero, we draw no edge Similar to SEM, a model can be devised that perfectly
between two nodes. As such, the GGM can be seen as a explains this pattern using only five parameters, because
network model of conditional associations; no edge indi- one of the elements in K can be constrained to be zero
cates that two variables are independent after condition- (Epskamp et al., 2017d). We can now standardize this
ing on all other variables in the data set. This allows us to matrix and make the off-diagonal elements negative
model conditional associations, which we might expect (Equation (2)) to obtain the partial correlation matrix,
to be zero, rather than marginal associations, which we which we will denote R :
rarely expect to be zero (Meehl, 1990). ⎡ ⎤
1 −0.25 0.3
To exemplify the above, suppose for three variables R = ⎣ −0.25 1 0 ⎦.
“fatigue,” “concentration problems,” and “insomnia,” the 0.3 0 1
true variance-covariance matrix is
This matrix can be used to draw a network as is shown
⎡ ⎤ in Figure 1. This figure shows that someone who is tired
1 −0.26 0.31
is also more likely to suffer from concentration problems
 = ⎣ −0.26 1 −0.08 ⎦.
and insomnia. Furthermore, this network shows that the
0.31 −0.08 1
correlation between insomnia and concentration prob-
lems can be explained by the relationships of both vari-
To model this matrix, we need six parameters (three ables with fatigue: concentration problems and insomnia
covariances and three variances). The corresponding true are conditionally independent given the level of fatigue.
Interpreting GGMs. This paper concerns the
exploratory estimation of GGMs from various sources of
 Many publications make use of the default color setup used in qgraph data, without prior knowledge on the model structure.
(Epskamp, Cramer, Waldorp, Schmittmann, & Borsboom, ): green for pos- Such undirected network models can be interpreted
itive edges and red for negative edges. A later version of qgraph includes the
option theme = ”colorblind” using a more colorblind friendly col- in strikingly different ways, ranging from a no causal
oring scheme and setting the positive edge color to blue. This option has interpretation to a strong causal interpretation:
been used for all graphs in this paper. Note that some publications (e.g., Schu-
urman, ) also use blue and red edges but use red to denote positive and
(1) Predictive effects. The GGM can be interpreted
blue to denote negative effects akin to a heat map. without any causal interpretation and used merely
456 S. EPSKAMP ET AL.

as a tool to show which variables predict one- correspond to multiple regression coefficients and next
another. Interpreting the parameters associ- discussing the relationships between the GGM and SEM.
ated with the model A → B → C requires a Point 3 follows from observing that the GGM is directly
causal interpretation, while the predictive quality related to similar undirected models such as the Ising
between these nodes can directly be obtained from model (Ising, 1925). A discussion on the causal interpre-
the equivalent GGM A–B–C: only information tation of such models is beyond the scope of this paper,
on node B is needed when predicting A or C. As and we refer the reader for this topic to Epskamp et al. (in
such, the GGM can always be interpreted to show press) and van Borkulo et al. (2014).
predictive effects and offers a powerful exploratory
tool to map out multicollinearity.
(2) Indicative of causal effects. The GGM is closely
1.1. The Gaussian graphical model and multiple
tied to causal modeling. If a causal model between regressions
observed variables generated the data, then an An edge in a GGM indicates that one node predicts a
edge A – B appears in the GGM only if there is a connected node after controlling for all other nodes in
causal link between the variables (e.g., A → B or the network. This can also be shown in the relationship
A ← B), or if both variables cause a third variable between coefficients obtained from least-squares predic-
in the data (e.g., A → C ← B). Exploratory esti- tion and the inverse variance-covariance matrix. Let 
mation of such models relies on stringent assump- represent an k × k matrix with zeros on the diagonal. Fur-
tions (e.g., acyclicity), suffers from a problem of thermore, let γ i,−(i) represent the ith row of  without the
many equivalent models, and may lead to over- ith element (as the diagonal is set to zero), which contains
saturated models. The GGM, on the other hand, the regression coefficients obtained in a multiple regres-
is well identified and does not feature equivalent sion model:
models. Therefore, at the cost of losing informa-
tion on the direction of effect, exploratory search yci = τ + γ i,−(i)y c,−(i) + εci .
algorithms perform well in identifying a GGM. As such, γi j encodes how well the jth variable predicts
Because of this close tie to causal modeling, edges the ith variable. This predictive effect is naturally sym-
in the GGM may be interpreted as indicative of metric; if knowing someone’s level of insomnia predicts
potential causal pathways. his or her level of fatigue, then conversely knowing some-
(3) Causal generating model. Undirected network one’s level of fatigue allows us to predict his or her level of
models have a long history of being used as data- insomnia. As a result, γi j is proportional to γ ji . There is
generating models in diverse scientific fields such a direct relationship between these regression coefficients
as statistical physics (Murphy, 2012). For exam- and the inverse variance-covariance matrix (Meinshausen
ple, in a simple ferromagnetic Ising model of two & Bühlmann, 2006). Let D denote a diagonal matrix on
particles that tend to be aligned (Epskamp et al., which the ith diagonal element is the inverse of the ith
in press), A — B, intervening on A would impact residual variance: dii = 1/Var(εCi ). It can then be shown
B and intervening on B would impact A. To this (Pourahmadi, 2011) that7
end, undirected network models allow for a unique
causal interpretation: one of genuine symmetric K = D (II −  ) . (3)
effects. This interpretation is discussed often in the
Thus, κi j is proportional to both γi j and γ ji ; a zero in
literature on network psychometrics, and used in
the inverse variance-covariance matrix indicates that one
complexity research demonstrating emergent phe-
variable does not predict another. Consequently, the net-
nomena (e.g., the positive manifold or phase tran-
work tells us something about the extent to which vari-
sitions) that may occur in such a network of cel-
ables predict each other. This predictive quality is the cor-
lular automata (Cramer et al., 2016; Dalege et al.,
nerstone for how such network models are often applied
2016; Kruis & Maris, 2016; van der Maas et al.,
(Hastie, Tibshirani, & Wainwright, 2015), for example, in
2006).
recommender-systems that recommend users on prod-
In addition, the GGM is closely tied to factor analy-
ucts they might like depending on which products the
sis, allowing for extensions to factor modeling through
user already liked (Marsman, Waldorp, & Maris, 2017b).
the use of network modeling (Epskamp et al., 2017d). The
main focus of this paper is discussing the second interpre-  This expression may differ by a scalar, depending on the estimation method.
tation, while also describing how the GGM may be used For example, by default R computes the variance-covariance matrix by using
n − 1 in the denominator, but computes Var(εCi ) by using n − m in the
to show predictive effects. We detail these first two points denominator. This denominator is cancelled out in Equation () when stan-
below by first showing how partial correlation coefficients dardizing to partial correlation coefficients.
 MULTIVARIATE BEHAVIORAL RESEARCH 457

In addition to these applications and aiding the interpre-  contains m(m + 1)/2 elements, while contains m
tation of GGM models, this relationship between multiple parameters and B contains m(m − 1) parameters. As a
regression and undirected network edges plays a crucial result, the model above is underidentified without strin-
role in many network estimation procedures (Haslbeck gent restrictions on B . One assumption is that yC can
& Waldorp, 2016b; Meinshausen & Bühlmann, 2006; van be ordered such that B is lower triangular, indicating
Borkulo et al., 2014), including the methods discussed that if this matrix is used to draw a directed graph—a
below in this paper. graph in which A → B indicates that A causes B—that
graph does not contain any cycles, meaning that directed
edges cannot be traced from any node back to itself (e.g.,
1.2. The Gaussian graphical model and
A → B → A). Such a graph is called a directed acyclic
structural equation modeling
graph (DAG; Kalisch & Bühlmann, 2007; Pearl, 2000). If
Let ηC represent a set of unobserved variables, which we repeated measures are available at the correct time scale,
assume to be jointly normally distributed with yC . Then, reciprocal effects and cycles can often be adequately mod-
we can form an encompassing framework for several pos- eled as acyclic effects unfolding over time. Without such
sible generating models:8 information, cycles can be modeled and can be identi-
fied when exogenous variables are present (such as the
y c = By c + η c + ε c weather, time, or, depending on the modeling framework,
ε C ∼ N(00, ) lagged variables; Rigdon, 1995), but the interpretation
ηC ∼ N(00, ), (4) of such cycles is not without problems (Hayduk, 2009).
Several software packages exist that aim to find such a
in which is a diagonal matrix,9 indicating that after con- DAG (e.g., pcalg, Kalisch, Mächler, Colombo, Maathuis, &
ditioning on all causes the variables are independent, B Bühlmann, 2012; bnlearn, Scutari, 2010). The assumption
is a square matrix with zeros on the diagonal of causal of acyclicity, however, is debatable in the context of psy-
effects between observed variables, and is a factor- chological variables (Schmittmann et al., 2013) because
loading matrix. The variance-covariance matrix of η C may many effects can be plausibly assumed cyclic (e.g., fatigue
in turn be modeled in various ways to achieve compli- → concentration problems → stress → fatigue).
cated model setups. The expression above is well-known Second, the same structure for  can be obtained
in SEM, which allows for confirmatory testing of causal under various different specifications of B . Thus, many
models. In exploratory estimation, one could assume no equivalent models can lead to exactly the same fit. This
latent variables exist and aim estimate B (causal models), can be seen because several matrix decompositions of
or one could assume no relationships between observed  , such as a Cholesky decomposition or an eigendecom-
variables exist and aim to estimate (factor models). We position, can be used to produce equivalently fitting B .
contrast both to the GGM below. The problem of equivalent models is also well-known in
the literature on directed networks and SEM (MacCal-
1.2.1. Causal models lum, Wegener, Uchino, & Fabrigar, 1993; Pearl, 2000). For
example, the following three causal models are not statis-
Suppose there are no unobserved causes to any of the tically distinguishable:
variables in yC , and the variables in yC are only caused (1) Concentration → Fatigue → Insomnia
by other variables in yC . The corresponding model for  (2) Concentration ← Fatigue → Insomnia
becomes: (3) Concentration ← Fatigue ← Insomnia
 = (II − B )−1 (II − B )−1 . (5) All three models only imply that concentration and
insomnia are conditionally independent given fatigue.
In this expression, B can now be seen to encode the causal With more variables, the number of potential equiva-
model (Pearl, 2000). Table 1 summarizes the comparison lent models increases drastically, making it evident that
between such causal models and GGMs. Although use- model search is likely to fail. At best, exploratory esti-
ful for generating data and confirmatory testing, we can mation can result in a set of equally plausible DAGs (an
see two problems in exploratory estimation of B without equivalence class; Drton & Maathuis, 2017), each differ-
any prior knowledge. First, if m variables are included, ently parameterized and each leading to different strong
causal hypotheses.
 This expression should not be confused with Equation (), in which  is
obtained by performing univariate multiple regressions in which error terms Causal modeling and the GGM. The undirected
are not independent. GGM offers an attractive alternative to exploratory DAG
 This matrix is often denoted using the Greek letter instead of . We use
here to avoid confusion with the contemporaneous variance-covariance
estimation: the GGM is saturated rather than overi-
matrix used below, which is not diagonal. dentified if all edges are present (K K contains the same
458 S. EPSKAMP ET AL.

Table . Overview of causal models (directed networks) and Gaussian graphical models (undirected networks).
Causal model Gaussian graphical model

 −1 = (II − B ) −1(II − B) = K

A C A C

+ + + +
B B
A ⊥⊥ C | B

A C A C

+ + + +
B B
A ⊥⊥ C | B

A C A − C

+ + + +
B B
A ⊥⊥ C | B
R packages (confirmatory) Any SEM package lvnet (fit measures); qgraph (fit measures); ggm (estimation
only); GLASSO (estimation only)
R packages (exploratory) pcalg; bnlearn qgraph (EBICglasso function); GLASSO (no automatic
tuning parameter selection); huge; parcor; BDgraph;
lvnet (for GGM at latent or residual level of SEM)
Pros Causal interpretation; allows for confirmatory testing of No equivalent models; fast structure and parameter
causal hypotheses; can detect common effect structures estimation using LASSO; edges parameterizable as
partial correlation coefficients; edges interpretable as
predictive effects; latent variables result in clusters;
edges can be indicative of potential causal effects
Cons Exploratory estimation requires assumption of acyclicity; No direction of effect; common effect structure can induce
many equivalent models; direction of effect poorly or spurious edge; LASSO estimation assumes true model is
not identified; strongly depends on assumption of no sparse
latent variables

number of unique elements as  ), does not feature equiv- edge between node i and j (e.g., Yi → Y j or Yi ← Y j ) and
alent models10 (there is only one unique inverse for  ), if there is no common effect of node i and node j (e.g.,
does not suffer from a questionable direction of causal Yi → Yk ← Y j ; Koller & Friedman, 2009).11 Thus, assum-
effect, does not require the assumption of acyclicity, and is ing a causal model as in Equation (5) generated the data,
easily parameterized using partial correlation coefficients an edge in a GGM emerges as a result of a direct causal
(Epskamp et al., 2017d). These benefits come, however, at effect between the variables, or as a result of the fact that
the cost of losing information on the direction of effect. both variables have a common effect on a third variable.
To investigate the structure of a GGM under the causal Note that within the causal model of Equation (5), there
model of Equation (5), in which observed variables can are no latent common causes by assumption. Edges in
only be caused by other observed variables, we can invert the GGM can therefore be indicative of potential causal
that expression to obtain: effects.
A note on common effects. As mentioned above and
K = (II − B ) −1
(II − B ) , (6)
shown in Table 1, conditioning on a common effect may
in which −1 is still a diagonal matrix. It becomes evident induce a spurious edge in the GGM. In this case, the
that there is no longer a matrix inversion needed and that sign of the edge can be informative: two positive causal
the sparsity in B directly corresponds to the sparsity in K ; effects from two variables on a third lead to a negative
the GGM thus acts on the same level as causal modeling. partial correlation. As such, when observing an edge of an
We can derive that κi j equals zero if there is no directed unexpected sign in the GGM, this may be indicative of a
common effect, especially when the marginal correlation
 Note that the uniqueness of the GGM relates to the psychometric model: for
every  there is only one unique inverse K and vice versa. When estimating a
GGM from data, different estimation methods may lead to different estimated  A common effect node is also termed a “collider” in the literature on causal
GGMs modeling.
 MULTIVARIATE BEHAVIORAL RESEARCH 459

coefficient between the two variables was of the different modeling can further be used to augment factor analy-
sign. It should also be noted that conditioning on a com- sis by modeling the latent variable variance-covariance
mon effect might cancel out a weak effect between two matrix or the residual variance-covariance matrix
variables. In addition, because edges may be induced due as a GGM (Epskamp et al., 2017d). Modeling as a
to conditioning on a common effect, the GGM does not GGM leads to a latent network model, which can be
estimate a skeleton graph, a causal network with arrow- used in exploratory estimation of relationships between
heads removed (an edge may be in the GGM that is not in latent variables. Modeling as a GGM leads to a resid-
the causal model). Skeleton graphs can also be estimated ual network model, which may be used to estimate factor
from data (e.g., Kalisch, Maechler, & Colombo, 2017), but models while local independence is structurally violated
are not parameterized and rely on many separate con- (Epskamp et al., 2017d; Pan, Ip, & Dubé, 2017).
ditional independence tests, potentially leading to power A note on spurious edges. When interpreting edges
issues. 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,
1.2.2. Factor models and cannot be accounted for unless the latent variable
Suppose that, instead of assuming no unobserved causes is explicitly modeled (e.g., by using the residual network
as in (5), we take the generating model of (4) and model described above). The same problem occurs in
only allow for unobserved causes of the observed vari- exploratory DAG estimation, in which case a latent vari-
ables. Then, (4) reduces to the well-known factor model able may induce a directed edge in the causal network.
(Brown, 2014). The corresponding model for  now Furthermore, such spurious associations may arise in any
becomes: 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
which can subsequently be inverted to obtain an expres- retrieve the direction of effect, however, when an edge is
sion for the equivalent GGM. Golino & Epskamp (2017) spurious due to the influence of a latent variable, the GGM
provide a detailed derivation of this inverted expression also does not introduce a strong causal hypothesis on what
and show that a factor in the factor model will lead to its would happen under intervention.
indicators to cluster (all nodes connected to each other
with strong edges) in the GGM. This result is in line with
2. Estimating GGMs from different sources of
mathematical equivalences between factor models and
data
network models of binary variables (Epskamp et al., in
press; Kruis & Maris, 2016; Marsman et al., 2017a; Mars-
2.1. Data with independent cases
man, Maris, Bechger, & Glas, 2015). As there is only one
unique inverse to  , there is only one unique GGM for A GGM can be estimated in data sets where cases can be
every factor model. Conversely, however, one GGM may assumed to be independent. Three common examples of
be equivalent to many different factor models (e.g., all such data are cross-sectional data, in which every subject
possible rotations of ). is only measured once on a set of response items, aggre-
Due to these equivalences, network modeling and fac- gated data, in which only one mean score per variable per
tor modeling are closely connected. A natural first step subject is included in the data set, or n = 1 time-series
in performing an exploratory factor analysis would be to data that feature large intervals between measurement
estimate and draw a GGM model and investigate if the occasions. In time-series data featuring shorter intervals,
nodes cluster as would be expected by a factor model. a GGM can be estimated as well; in this case, the network
Cluster-detection algorithms on the GGM could even be could be termed a contemporaneous network. However,
performed to investigate the number of factors to extract as we argue in the next section on temporally ordered
(Golino & Epskamp, 2017). Of note, however, is that data, better methods exist that take temporal information
many GGM estimation methods will always aim to esti- into account in addition to modeling the contemporane-
mate sparse GGM (i.e., K contains exact zeroes), which ous effects in a GGM.
is not expected given a factor model (except when latent
variables are orthogonal). As such, estimating a sparse ... Estimation
network does not provide evidence that a latent vari- In cross-sectional data analysis, only one observation per
able model could not have generated the data (Epskamp, subject is available; thus, we cannot expect to estimate
Kruis, & Marsman, 2017c; Epskamp et al., in press). GGM subject-specific means or GGM networks. It is typically
460 S. EPSKAMP ET AL.

assumed that the subjects all share the same distribution. any programming language and in many statistical pro-
That is, grams by inverting and subsequently standardizing the
sample variance-covariance matrix. In the open-source
y P ∼ N (00,  ) , statistical programming language R (R Core Team, 2017),
automated procedures have been implemented in the
in which y P denotes the random response of subject P on
corpcor package (Schafer et al., 2017) and the qgraph
all items. Similarly, in n = 1 time-series data we can make
(Epskamp et al., 2012) package. The qgraph package
a similar assumption:
also supports thresholding via significance testing or
y T ∼ N (00,  ) , false discovery rates. The GLASSO algorithm is imple-
mented in the glasso (Friedman, Hastie, & Tibshirani,
in which y T denotes the random response of a subject 2014) and huge (Zhao et al., 2015) packages. EBIC-
on all items at time point T . In both cases, the full based tuning parameter selection using the glasso package
likelihood can be readily obtained, and the variance- has been implemented in the qgraph package. The huge
covariance matrix  can reliably be estimated using max- package also allows for selection of the tuning parame-
imum likelihood estimation (MLE), least-square estima- ter using cross validation or EBIC. The parcor package
tion, or Bayesian estimation. (Krämer, Schäfer, & Boulesteix, 2009) implements other
Regularization. The MLE solution of K —the precision LASSO variants to estimate the GGM. The BDgraph pack-
matrix encoding a GGM—can be obtained by standard- age (Mohammadi & Wit, 2015) implements a Bayesian
izing the inverse sample variance-covariance as per Equa- method to estimate the undirected structure. Finally, fit-
tion (2). To obtain a sparse network, model search can be ting an estimated GGM to data can be done in the R pack-
performed by iteratively adding and removing edges and ages ggm (Marchetti, Drton, & Sadeghi, 2015) and lvnet
fitting the corresponding GGM structure (Epskamp et al., (Epskamp et al., 2017d).
2017d). In recent literature, it has become increasingly
popular to use regularization techniques, such as penal-
2.2. Temporally ordered data of a single subject
ized MLE, to jointly estimate model structure and param-
eter values (Costantini et al., 2015; van Borkulo et al., In line with a call for more intraindividual and person-
2014). The least absolute shrinkage and selection operator based research (Molenaar, 2004), an increasingly
(LASSO; Tibshirani, 1996) has been shown to perform popular form of data pertains to n = 1 time series,
well in quickly estimating model structure and parame- in which a single individual is measured repeatedly over
ter estimates of a sparse GGM (Friedman, Hastie, & Tib- a period of time. One such situation is in clinical prac-
shirani, 2008; Meinshausen & Bühlmann, 2006; Yuan & tice (Epskamp et al., 2018; Kroeze et al., 2017), where
Lin, 2007). A particularly popular variant of LASSO is a patient can be measured several times per day over a
the graphical LASSO (GLASSO; Friedman et al., 2008), period of a few weeks. We will limit our discussion to
which directly penalizes elements of the inverse variance- data obtained in a relatively short time-frame so that
covariance matrix (Witten, Friedman, & Simon, 2011; we can reasonably assume the model will remain stable
Yuan & Lin, 2007). The GLASSO algorithm is useful as it over time. Then, we can apply the methodology above to
is typically faster than other GGM estimation algorithms obtain a GGM for the n = 1 data. However, such an anal-
(which conduct multiple separate regressions and then ysis does not take temporal ordering of data into account
combine the results using Equation (3)), and requires only (i.e., relationships between measurement occasions)
an estimate of the variance-covariance matrix rather than and only investigates contemporaneous relationships
raw data (Epskamp & Fried, in press). LASSO utilizes a between variables (e.g., within the same measurement
tuning parameter which can be chosen in a way that opti- occasion). This is important for several reasons. First,
mizes cross-validated prediction accuracy or that mini- valuable information, especially in the context of dynam-
mizes information criteria such as the extended Bayesian ical relationships, might be contained at the temporal
information criterion (EBIC; Chen & Chen, 2008). Esti- level rather than at the contemporaneous level. Second,
mating a GGM with the GLASSO algorithm in combina- not taking temporal ordering into account might bias the
tion with EBIC model selection has been shown to work estimated contemporaneous relationships (see Section
well in retrieving the true network structure (Epskamp, 4 of the supplementary materials). For example, if one
2016; Foygel & Drton, 2010). For an introduction to this variable causes itself and another variable at the next time,
methodology aimed at empirical researchers, we refer the then not taking temporal ordering into account turns
reader to Epskamp & Fried (in press). that variable into a latent cause, which would produce
Software. Several software packages allow for GGM an edge in the GGM. Third, temporal information is
estimation as described above. MLE can be performed in needed when constructing the joint likelihood over time
 MULTIVARIATE BEHAVIORAL RESEARCH 461

(e.g., to obtain the information retained in a system over the effect). Temporal networks may thus highlight poten-
time; Epskamp, 2017a; Quax, Kandhai, & Sloot, 2013). tial causal pathways. While temporal networks are typi-
Finally, temporal information can aide in distinguishing cally cyclic, they can also be interpreted as summarizing
reciprocal and cyclic effects by regarding these as acyclic a DAG unfolding over time.
effects unfolding over time. Contemporaneous networks. In addition to tempo-
Vector Auto-regression. The simplest way to deal with ral effects, VAR analyses also include contemporaneous
temporal ordering of cases is to incorporate the effect effects, which can be modeled as a GGM. We will term
between consecutive measurements (Chatfield, 2016; this modeling framework (a VAR model with contem-
Hamilton, 1994; Shumway & Stoffer, 2010). This is called poraneous effects explicitly modeled and portrayed as a
a Lag-1 model because it includes both measurements at GGM) graphical VAR (GVAR; Wild et al., 2010).13 A use-
the current time point t as well as measurements from the ful equivalent way to denote a GVAR model is by using a
previous time point t − 1. We will focus our discussion on conditional Gaussian distribution:
Lag-1 models, noting that everything below also general-
izes to more complicated models (e.g., Lag-2 models). In  
y T | y T −1 = y t−1 ∼ N By t−1 , .
intraindividual analysis, VAR (Brandt & Williams, 2007;
Rosmalen, Wenting, Roest, de Jonge, & Bos, 2012) has
Which is equivalent to Equation (7). Now, it becomes evi-
gained substantive footing in visualizing temporal infor-
dent that if consecutive cases can be assumed to be inde-
mation through networks. The lag-1 VAR model can be
pendent, and thus B = 0 , the GVAR model is exactly the
denoted as a regression model on the previous measure-
same as the GGM model described above for independent
ment occasion:
cases. Thus, the GVAR model can be seen as a general-
y t = By t−1 + ε t ization of the GGM model to temporally ordered data.
GVAR only differs from regular VAR in that the con-
ε T ∼ N(00, ). (7)
temporaneous structure is modeled and represented as a
GGM, instead of being saturated. This leads to a strikingly
The model matrix B encodes temporal predictive effects
different interpretation of the VAR model; the VAR model
from variables on variables in the next measurement occa-
can be seen as an inclusion of temporal effects on a GGM.
sion, and can be used to obtain a directed network, which
Temporal and contemporaneous information.
we term the temporal network. The variance-covariance
K ( ) = −1 ) to obtain a GGM Figure 2 shows a hypothetical example of the two net-
matrix can be inverted (K
work structures obtained in a GVAR analysis and shows
modeling effect within the same measurement occasion,
how they might plausibly differ. The left panel shows the
after controlling for temporal effects. These can be dis-
temporal network. The self-loop shows that whenever the
played again as a network, which we term the contempo-
subject in question felt energetic (or tired), this person
raneous network.
also felt more (or less) energetic in the next measurement.
Temporal networks. Temporal networks, encoded by
The temporal network also shows us that after exercising,
B , have grown popular in recent psychological literature
this person felt less energetic. The contemporaneous
(e.g., Bringmann et al., 2013, 2015; Bos et al., 2017; Klip-
network in the right panel shows a plausible reverse
pel et al., 2017; Snippe et al., 2017; Wigman et al., 2015). A
relationship: Whenever this person exercised, he or she
temporal network is formed by combining a lagged vari-
felt more energetic in the same measurement occasion.
able yt−1 and current variable yt into a single node, con-
In psychology, there will likely be many causal relation-
nected with directed edges which are weighted according
ships that occur much faster than the lag interval of a
to the regression parameters contained in B .12 Thus, an
typical ESM study; in this case, these pathways will be
edge in the temporal network indicates that a node pre-
captured in the contemporaneous network. For example,
dicts another node (or itself in the common case of self-
if someone is experiencing bodily discomfort, that will
loops) at the next measurement occasion, after control-
immediately negatively affect that person’s ability to enjoy
ling for all other variables at the previous measurement
him or herself (Epskamp et al., 2018). Especially when
occasion. Temporal prediction is central to the concept
the measurement is on blocks of time (e.g., “since the last
of Granger causality in the economic literature (Eichler,
measurement did you feel ...”), such effects are likely to be
2007; Granger, 1969), and it satisfies at least the tempo-
caught in the contemporaneous network.
ral requirement for causation (i.e., the cause must precede

 Note, in graph theory it is common to encode a network using a weights  Wildet al. () do not use the term graphical VAR in the exact same way
matrix in which the row indicates the node of origin and the column indi- we do, and use it more to refer to graphical modeling in a VAR framework,
cates the row of destination. As such, to obtain the directed weights matrix including structural VAR. We use the term here as described because having
to draw a temporal network B needs to be transposed. an explicit term helps in contrasting GVAR from, e.g., structural VAR.
462 S. EPSKAMP ET AL.

Temporal network Contemporaneous network

Exercising Exercising

Energetic Energetic

Figure . A hypothetical example of two network structures obtained from a GVAR analysis. The network on the left indicates the tem-
poral 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.

2.2.1. Estimation variance-covariance matrix (these can be stored using the

BPARAMETERS option).
Estimating saturated (fully connected temporal
When estimating GVAR models, regularization meth-
and contemporaneous networks) GVAR models is
ods can be used similar to the estimation of GGMs on
straightforward. First, one needs to estimate temporal
nontemporally ordered data. Abegaz and Wit (2013) pro-
effects of a regular VAR model by performing multivari-
posed to apply LASSO estimation to jointly estimate the
ate multiple regression of all variables on the previous
temporal and contemporaneous network structures using
measurement occasion,
the multivariate regression with the covariance estimation
y t = By t−1 + ε t , (MRCE) algorithm described by Rothman, Levina, & Zhu
(2010). MRCE involves iteratively optimizing B , using
or by estimating univariate models for every variable, cyclical-coordinate descent, and K ( ) , using the GLASSO
algorithm (Friedman et al., 2008, Friedman et al., 2014).
yti = β iy t−1 + εti , EBIC model selection can be used to obtain the best per-
forming model. This methodology has been implemented
in which β i denotes the ith row of B . Next one can invert in two open source R packages: sparseTSCGM (Abegaz &
the variance-covariance matrix of the residuals to obtain a Wit, 2015), which aims to estimate the model on repeated
GVAR model. Step-wise model selection in latent network multivariate genetic data, and graphicalVAR (Epskamp,
models (Epskamp et al., 2017d) can also be used to esti- 2017c), which was designed to estimate the model on
mate sparse GVAR models. Missing data can be handled the psychological data of a single subject. The graphical-
in default ways of SEM or regression models (e.g., listwise VAR package also allows for unregularized multivariate
deletion and full-information maximum likelihood), or estimation.
by using more sophisticated techniques such as Bayesian An alternative to estimating GVAR models is to esti-
estimation (Schuurman, Grasman, & Hamaker, 2016b) or mate structural VAR (SVAR; Chen et al., 2011) models,
the Kalman filter (Harvey, 1990; Kim & Nelson, 1999). also called unified SEM (Gates, Molenaar, Hillary, Ram,
Novel estimation methods. A promising recent & Rovine, 2010). In SVAR, the contemporaneous effects
method for estimating VAR models is the Bayesian are modeled using a directed network instead of an undi-
dynamical SEM implementation in version 8 of Mplus rected network. The sparsity of the undirected GVAR
(Asparouhov, Hamaker, & Muthén, 2017; Muthén & contemporaneous network corresponds in the same way
Muthén, 2017), which includes handling of missing data, to the sparsity of the directed contemporaneous network
measurement invariance, and latent variables. Mplus in SVAR as how the GGM corresponds to causal models
can be used to estimate saturated GVAR models, and to (edges arise in the GGM due to edges in the causal net-
perform model selection in the temporal network of a work or conditioning on common effects). The temporal
GVAR model. Model selection in the contemporaneous SVAR network is sparser than the temporal GVAR net-
network of a GVAR model is not yet implemented, but work, as contemporaneous mediators can be controlled
credibility intervals around contemporaneous effects can for in SVAR but not in GVAR. A saturated SVAR model
be obtained by manually inverting each sampled residual can be obtained by using regressions on the previous
 MULTIVARIATE BEHAVIORAL RESEARCH 463

time-point as mentioned above, followed by transforming these fixed effects, such as B p − B ∗ , are often called
the contemporaneous variance-covariance matrix (e.g., random effects. Besides the individual network struc-
by using a Cholesky decomposition on its inverse; Lütke- tures, researchers often aim to estimate the structure
pohl, 2005) and subsequently transforming the temporal and parameters of these fixed effects because these tell
effects to take contemporaneous mediators into account. us something about the average intraindividual effect.
This technique of obtaining an SVAR model leads to Researchers also aim to estimate the variance-covariance
multiple solutions (Beltz & Molenaar, 2016). Step-wise structure of the random effects because it tells us some-
model selection can also be used to estimate sparse SVAR, thing about individual differences (Bringmann et al.,
for example, by using model selection in (unified) SEM 2013).
(Gates et al., 2010) or Bayesian dynamical SEM models The random effects can be modeled by assuming a sec-
(Asparouhov et al., 2017; Muthén & Muthén, 2017). ond 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
2.3. Temporally ordered data of multiple hierarchical model for K ( ) needs to take into account
subjects that this matrix must remain positive definite). The inter-
A type of data that is increasingly common due to the pretation of, for example, correlations between differ-
emergence of ESM studies is time series of multiple sub- ent temporal or contemporaneous edges is also difficult.
jects (e.g., Bringmann et al., 2013, 2015; Mõttus et al., Therefore, we only focus here on a subset of the parame-
2017; Schmiedek, Lövdén, & Lindenberger, 2010; Wig- ters where we can easily interpret the second-level model:
man et al., 2015). Such data sets pose a promising gate- the mean structure. As a result, if a multivariate normal is
way to study both intraindividual dynamics and between- assumed for all parameters, then it is also assumed for the
subjects overlap as well as their differences. Here, we marginal distribution of the means—regardless of other
assume that the number of time points might differ per parameters:
person and that measurement occasions are nested in
people. We can model the temporal data of every person μ P ∼ N (00,  ) .
with an individual GVAR model: Again, we can invert the variance-covariance matrix to
 
y [t,p] = μ p + B p y [t−1,p] − μ p + ε [t,p] obtain a GGM,
ε [T,p] ∼ N(00, p ) )
K ( =  −1 ,
−1
p = K (p ) ,
which we will term the between-subjects network, a net-
in which μ p indicates the stationary mean vector of sub- work between stationary means of different subjects.14 As
ject p (which enters the model because we can no longer such, estimating the GVAR model on n > 1 time-series
assume within-subject means are zero without loss of analysis allows for the separation of variance into three
generality), B p encodes the person-specific temporal net- distinct network structures: temporal networks, contem-
work, and K (p ) encodes the person-specific contempora- poraneous networks, and the between-subjects network.
neous GGM.
Multilevel modeling. To gain insight in the general
network structure over subjects, we can investigate the 2.3.1. Estimation
individual networks at a second level. Doing so is termed In this section, we outline several different ways in which
multilevel modeling, explained in more detail in Section individual network structures as well as fixed effect net-
2.1 of the supplementary materials. Let B ∗ and K ∗( ) work structures may be estimated. We first discuss apply-
encode the expected temporal and contemporaneous net- ing the methodology of estimating n = 1 GVARs dis-
work when selecting a person at random. Furthermore, cussed above to both pooled data as well as data of each
we can assume without loss of generality that data are subject separately, followed by a discussion of different
grand-mean centered. We then obtain: multilevel estimation procedures that take clustering of
E (μμP ) = 0
 Ofnote, it is also possible to invert and standardize the full random effects
E (BBP ) = B ∗
 variance-covariance matrix, which would lead to different between-subjects
relationships between the means as well (partial correlations after condi-
( )
E KP = K ∗( ) . tioning on other means and all other between-subject parameters such as
edges). We do not do that here as (a) such a network is hard to interpret,
Here, B ∗ and K ∗( ) now encode the average parame- and (b) most estimation methods we mention do not return the full random
effects variance-covariance matrix (especially the correlations between tem-
ters in the population: the fixed effects. Deviations from poral and contemporaneous edges are hard to obtain).
464 S. EPSKAMP ET AL.

Table . Three methods of estimating GVAR models with n > 1subjects.

Pooled and individual LASSO estimation Bayesian multilevel Two-step frequentist multilevel

Software graphicalVAR (Epskamp, c); MPlus  (Muthén & Muthén, ; mlVAR (Epskamp, Deserno, & Bringmann,
sparseTSCGM (Abegaz & Wit, ). Asparouhov et al., ); mlVAR b).
(wrapper around Mplus).
Estimation () Joint multivariate LASSO estimation MCMC sampling from multivariate () Sequential univariate multilevel
with EBIC model selection (Abegaz & hierarchical model (e.g., Schuurman regression models on previous
Wit, ) of within-subjects centered et al., b). measurement (similar to Bringmann
data to obtain fixed effects temporal et al., ), with within-subject
and contemporaneous networks. () centered lagged variables as
GLASSO algorithm with EBIC model within-subjects level predictors and
selection (Foygel & Drton, ) on sample-means of all other variables as
sample means of subjects to obtain between-subjects predictor. ()
between-subject network. () Step () Sequential multilevel regression
repeated for each individual data set to models using the residuals of ():
obtain subject-specific networks. residuals of one variable are predicted
by residuals of all other variables in the
same measurement occasion.
Pros Fast estimation of fixed effects; scales up Borrowing information in individual Borrowing information in individual
well to large numbers of nodes; model network estimation from other network estimation from other
selection in individual networks; subjects; all model parameters and subjects; scales up well to  nodes
temporal and contemporaneous random-effect (co)variances can be (correlated random effects) or  nodes
networks obtained in the same analysis. estimated; credibility intervals can be (orthogonal random effects); many
obtained for edges and descriptive random effect variances correlations
statistics (e.g., centrality; density); can be estimated; fast estimation of
advanced extensions such as individual networks.
measurement error and latent variable
modeling possible; powerful handling
of missing values.
Cons Fixed effects estimated on pooled data; Relatively slow estimation, especially in Slow estimation in larger data sets; no
Subject specific networks estimated higher dimensional models; no model model selection (fixed effects can be
without borrowing information from selection (thresholding possible via thresholded using significance);
other subjects (no multilevel structure); credibility intervals); complicated to combination of many different models;
between-subjects network estimated in estimate contemporaneous random does not scale up well past  nodes;
a different model; very slow to estimate effects. poor handling of missing values.
subject-specific networks; poor
handling of missing values.

Note: The software listed only concerns user-friendly automated software because all these models could readily be implemented in most programming languages
or Bayesian sampler packages.

the data into account. An overview of these methods is Multilevel estimation. The second and third proce-
also included in Table 2. dures described in Table 2 make use of multilevel model-
Pooled and individual LASSO estimation. First, we ing (Hamaker, 2012). Two main benefits of this approach
can estimate a GVAR model for every subject to obtain are (1) instead of estimating the VAR model in each sub-
subject-specific estimates for the temporal and contem- ject, only the fixed effects and variance-covariance of the
poraneous networks. Similarly, we can estimate fixed- random effects need to be estimated, and (2) afterward,
effects networks by estimating a GVAR model on the estimates of subject-specific parameters can be obtained,
entire within-subjects centered data set, using the sam- which are somewhat pulled together (termed shrinkage).
ple means of every subject on every variable as a plug-in Shrinkage allows the estimation of the model for one sub-
for the within-subject means. Consequently, we can esti- ject to borrow information from other subjects. Multi-
mate the between-subjects network by estimating a GGM level estimation can be performed by specifying the mul-
on the sample means of each subject on all variables. tivariate model using hierarchical Bayesian Monte Carlo
We can readily apply the LASSO regularization meth- sampling methods or by integrating over the distribution
ods described earlier for this purpose: the methodology of the random effects (Gelman & Hill, 2006; Schuurman
outlined by Abegaz and Wit (2013) to estimate tempo- et al., 2016b).
ral and contemporaneous networks and the methodol- Multivariate Bayesian multilevel. Bayesian multivari-
ogy outlined by Foygel & Drton (2010) to estimate the ate estimation has proven to be powerful in estimat-
between-subjects GGM. We term this framework pooled ing multivariate multilevel models, especially given its
and individual LASSO estimation and have implemented flexibility in adding measurement error, latent variables,
it in the R package graphicalVAR (Epskamp, 2017c). The and in handling missing data (Schuurman, Houtveen, &
performance of pooled and individual LASSO estimation Hamaker, 2015). Recently, the dynamic SEM methodol-
is assessed in simulations reported in Section 3 of the sup- ogy implemented in Mplus version 8 (Asparouhov et al.,
plementary materials. 2017; Muthén & Muthén, 2017) has made estimation
 MULTIVARIATE BEHAVIORAL RESEARCH 465

of multivariate multilevel VAR models much faster and 2 of the supplementary materials. In short, we extend
more user-friendly than other Bayesian software routines. the methodology of Bringmann et al. (2013) by within-
Specifying a temporal VAR model with correlated ran- subject centering and by adding subject sample means as
dom effects is straightforward and relatively fast to com- between-subjects predictors (as discussed by, e.g., Curran
pute with a moderate number of variables (e.g., 6). At the & Bauer, 2011; Hoffman & Stawski, 2009; Hamaker &
time of writing, Mplus does not return partial correla- Grasman, 2014). This allows us to estimate between-
tions by default, but these can be obtained by using the subjects networks by collecting regression coefficients as
BPARAMETERS option and manually inverting the in Equation (3) and symmetrizing the resulting matrix.15
sampled variance-covariance matrices. Mplus allows In a second step, we take the residuals of the first anal-
for specifying random effects on the contemporaneous ysis and again perform sequential univariate multilevel
covariances and thus, by extension, allows for estimating regression models to predict each residual from all other
random contemporaneous networks in addition to ran- residuals in the same measurement occasion. Again, these
dom temporal networks. Specifying such a model can be can be collected, as in Equation (3), and symmetrized
done by specifying dummy latent variables for the residual to obtain contemporaneous networks. Networks can be
covariance between each pair of variables (a prior guess thresholded by removing all effects that are not signif-
on the sign of the covariance is needed). Doing so, how- icant. For the between-subjects and contemporaneous
ever, can significantly increase computation length espe- networks, this results in two p-values for every edge—
cially when all random effects are allowed to correlate. either both can be required to be significant (“and”-rule)
To facilitate estimation, we have implemented a function or an edge can be included if one of the two p-values is
generating Mplus code for a multilevel GVAR model and significant (“or”-rule). Using the “and”-rule means erring
subsequently running the model using the MplusAutoma- more on the side of caution (sparser network), whereas
tion package (Hallquist & Wiley, 2017) in version 0.4 of using the “or”-rule means erring more on the side of
the mlVAR package, which can be called using estima- discovery. We have implemented two-step multilevel
tor = ”Mplus” and requires the Mplus program to VAR in the mlVAR package, which can be called using
be installed. estimator = ”lmer” (the default).
Two-step multilevel VAR. A downside of multivari- Choosing the estimation method. The choice of
ate estimation is that the number of random effect covari- which estimator to use is not trivial and depends on
ances to be estimated increases quadratically with the the interests of the researcher. In Table 2, we list some
number of variables. Forcing random effects to be uncor- pros and cons of each of the methodologies. In par-
related helps, but places strict assumptions on the model. ticular, multilevel estimation can be very complicated
Bringmann et al. (2013) proposed to estimate multilevel and is harder in high-dimensional settings. Assuming
VAR models using univariate models instead, using a fre- normally distributed parameters can also be problem-
quentist estimation procedure. In this work, the multilevel atic because doing so imposes that subjects cannot
VAR model is estimated by sequentially estimating uni- differ on the structure of the networks, merely on their
variate multilevel regression models of one variable given parameterization. When a parameter (e.g., a temporal
all lagged variables. Doing so ignores several correlations edge) is zero in some subjects but nonzero in others, then
of random effects because many parameters are not esti- this parameter cannot be normally distributed (the dis-
mated in the same model, simplifying the analysis: only tribution would peak at 0). Therefore, it is currently hard
correlations between incoming edges to the same node to estimate differently structured individual networks
and the intercept of that node are included in an univari- (different edges set to be exactly 0 between subjects) in
ate model. This method scales up well to approximately multilevel estimation. Nonetheless, multilevel estimation
eight variables when estimating correlated random effects particularly shines in that when estimating an individual
and around 20 variables when estimating orthogonal ran- network, researchers can borrow information from other
dom effects (or by using a moving window approach; subjects. We have performed simulation studies to assess
Bringmann et al., 2015). Of note, when specifying orthog- the performance of the two proposed methods in this
onal random effects not all random effects are assumed to paper: pooled and individual LASSO estimation and two-
be uncorrelated, merely the ones used in the same uni- step multilevel VAR. We report the results of these studies
variate model. in Section 3 of the supplementary materials, which shows
The methodology of Bringmann et al. (2013) does that both methods adequately detect the true fixed-effect
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  Standardizing regression parameters from nodewise multilevel models to
of this estimation procedure are explained in Section partial correlation coefficients does not lead to perfectly identical estimates.
466 S. EPSKAMP ET AL.

network structures with increasing sample size. Hav- 3.1. Cross-sectional data analysis
ing more time points per subject helps to estimate the
Within- and between-subjects variation. A type of data
contemporaneous and temporal networks, and having
to which the GGM is currently often applied is data
more subjects helps to estimate the between-subject
belonging to multiple subjects that are all measured only
networks. Two-step multilevel VAR performs well in
once (e.g., Isvoranu et al., 2017; van Borkulo et al., 2015).
estimating intraindividual networks when the number of
Such a data set is often termed cross-sectional data,
observations is low, but does not perform subject-specific
and such an analysis is often termed a between-subjects
model selection: all estimated intra-individual networks
analysis. However, the term between-subjects analysis
are saturated and contain all edges. Pooled and aggre-
might not be warranted, as it is difficult to distinguish
gated LASSO estimation does estimate the structure of
between within-subject variation around an individual’s
intraindividual networks, but performs poorer in intrain-
stable mean and between-subject variation of such sta-
dividual parameter estimation with fewer observations as
ble within-subject means using only cross-sectional data
no information is borrowed from other subjects.
(Hamaker, 2012). It is well known that subjects might
GIMME. Finally, when analyzing n > 1 data, another
respond differently when measured multiple times (Lord,
option is to estimate SVAR models instead. A promising
Novick, & Birnbaum, 1968). As such, the single observa-
estimation procedure to estimate such models over many
tion per subject leads to the time point and the subject
individuals, while dealing with potential heterogeneity,
being random: y [T,P] . We might make the argument that
is “group iterative multiple model estimation” (GIMME;
two distinct sources of variation cause the outcome (Bol-
Gates & Molenaar, 2012), which is implemented in R
ger & Laurenceau, 2013). Repeated measures of a sub-
using the gimme package (Lane, Gates, Molenaar, Hal-
ject (here p) are distributed according to a unique within-
lquist, & Pike, 2016). In GIMME, no multilevel structure
subject model:
is imposed and subject-specific networks are allowed to
differ in structure. Information from other subjects is bor- μ p,  p ).
y [T,p] ∼ N(μ
rowed, however, in that the structure of individual net-
That is, of a particular response, the subject’s score is
works can be based on other subjects (e.g., an edge can be
a composite of the average stationary score μ p and
included because it is present in many other subjects). No
random deviation.16 These average stationary scores also
shrinkage is induced on the parameter estimates that are
differ in the population. Thus, we need to model the aver-
nonzero (as would be the case in multilevel or hierarchical
age stationary scores of a random subject P with a separate
Bayesian modeling). A variant of GIMME that estimates
distribution:
the GVAR or a combination of structural and GVAR mod-
els has not yet been developed, and we note that this may μ P ∼ N(00,  ),
be a promising avenue for further research.
in which we can assume, without loss of generality, an
overall mean of 0 . We can invert the variance-covariance
matrix  to obtain a GGM:
)
K ( =  −1 .
3. Interpreting GGMs estimated from
between-subjects data This GGM corresponds to a between-subjects net-
work. The matrix  p can also be inverted and
This paper describes the estimation of network models on
standardized to a GGM to obtain a within-subject
data from different subjects (both cross-sectional as well
network:
as person-wise average scores). Cross-sectional network
modeling is often criticized for inappropriately taking K p =  −1
p .
cross-sectional results to be reflective of within-person
We will term this network a within-subjects network.
causal processes (e.g., Bos & Wanders, 2016; Bos et al.,
The value of a cross-sectional analysis. It is immedi-
2017), as it can be shown that such results will not equal
ately clear that with only one response per subject we can-
within-person processes except under strong assump-
not hope to estimate subject-specific variance-covariance
tions (Molenaar, 2004). To this end, this section discusses
matrices  p (and as a result individual GGMs). More-
the interpretation of GGMs estimated from data of differ-
over, even if we assume that within-subject effects are
ent subjects. We first argue that cross-sectional data may
be interpreted as a between-subjects analysis, assuming  Section  of the supplementary materials shows that when consecutive cases

that between-subjects variance is dominant, and next dis- (t and t + 1) are assumed dependent, such a zero-order network may result
from a mixture of temporal and contemporaneous effects as described above.
cuss the interpretation of potential causal effects at the The discussion here does not concern estimation of model parameters and
between-subjects level. hence does not require an assumption of independence of cases.
 MULTIVARIATE BEHAVIORAL RESEARCH 467

equal across subjects (denoted with  ∗ below), this still not vary much over time; therefore, a cross-sectional anal-
leaves us without an estimable model because μ is also ysis seems more suitable here. Other examples are ques-
assumed to be normally distributed. The co-variation tionnaires asking participants to rate symptoms over a
between responses thus becomes an unidentified blend period of several weeks or to describe themselves as “I am
of  ∗ and  : A and B may correlate in cross-sectional a person who...” In such cases, the cross-sectional network
data because people who score on average high on A also may be interpreted as a between-subjects network, of
score on average high on B (trait-level variation in  ), which we discuss the interpretation below. Note that some
or because when people deviate from their average on A research indicates that state-level variance (how a person
they also tend to deviate from their average on B (state- feels at the moment) does influence self-reported scores
level variation in  ∗ ). Even when within- and between- given trait-level (how a person feels on average) instruc-
subjects effects are assumed not to correlate, the GGM tions (Brose, Lindenberger, & Schmiedek, 2013), indi-
estimated on such data becomes cating that this interpretation of cross-sectional results
should be taken on a case-by-case basis depending on the
topic studied.
 ∗ +  )−1 ,
K = (

3.2. Within- and between-subjects effects

which is not a simple function of the between-subjects
GGM and the within-subjects GGM. Only when no In contrast to prior work on multilevel VAR modeling
short-term within-subject variation,  ∗ = 0 , or no (e.g., Bringmann et al., 2013, 2015; Pe et al., 2015; Wig-
between-subjects variation,  = 0 , is assumed does the man et al., 2015), in this paper between-subjects effects
cross-sectional GGM correspond exactly to one of the two are conceptualized in addition to the within-subjects
networks. effects in a separate GGM. Furthermore, in contrast to
Cross-sectional data analysis thus cannot disentangle prior work, cross-sectional networks are not interpreted
between-subjects relationships from short-term within- to be reflective of within-subject effects, but rather to
subjects relationships (Hamaker, 2012). For example, potentially reflect a between-subjects structure, assuming
cross-sectional analysis cannot distinguish whether or that observed scores are not dominated by state-like vari-
not fatigue and concentration correlate because when- ance (but see Brose et al., 2013). This raises the question
ever people feel fatigued they also concentrate poorly (a on how such models could be interpreted. In particular, if
within-subjects effect) or because people who are on aver- edges in the GGM are interpreted as generating hypothe-
age fatigued also tend to concentrate poorly on average ses to potential causal pathways, the question is raised
(a between-subjects effect). However, preliminary sim- how such causal effects can occur at the between-subjects
ulation results show that the resulting cross-sectional level. This section therefore discusses the topic of causa-
GGM generally does not contain edges that are not tion at the between-subjects level. Here, we interpret the
present in either the within- or between-subjects net- stationary means as being locally stationary: the average
work (Epskamp, 2018). Depending on the ratio of within- of a subject in a relatively short time span of measure-
to-between person variance, the cross-sectional analysis ment (e.g., a few weeks). As such, we do not interpret
will pick up the within-subject network, the between- the mean vector μ P as a lifetime average. Instead, we
subject network, or a mixture of the two. As such, if assume it could change, potentially due to experimental
one assumes between-subject variance to be dominant, intervention. As a result, we argue that the between-
the cross-sectional results may be interpreted as mainly subjects network can also be indicative of potential causal
reflecting between-subjects relations. pathways—regardless of whether it is estimated from a
Cross-sectional analysis as between-subjects analy- cross-sectional interview concerning variables that are
sis. An important consideration is that a typical cross- not expected to vary much over time or obtained from
sectional questionnaire or interview is vastly different estimating the means from time-series data. To simplify
than a typical ESM questionnaire, and many cross- the argumentation below, we do not discuss separate tem-
sectional studies aim to measure variables that are more poral and contemporaneous networks but only general
stable over time and for which a time-series analysis might within-subjects networks (a GGM of within-subject data
not make sense. Good examples of this are recent network without taking temporal ordering into account).
analyses in the area of schizophrenia (Isvoranu, Bors- Simpson’s paradox. Hamaker (2012) described an
boom, van Os, & Guloksuz, 2016; Isvoranu et al., 2017), example of how within- and between-subject effects can
in which the impact of environmental factors (e.g., child- strongly differ from each other. Suppose we let peo-
hood trauma, urbanization) on psychotic symptoms and ple write several texts, and we measure the number of
general psychopathology was studied. Such variables do spelling errors they make and the number of words per
468 S. EPSKAMP ET AL.

Figure . 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 (). (b) Example based on Hoffman (); Hamaker ().

minute they type (typing speed). We would expect to see to a change in Y . Statistically, the interventionist account
the seemingly paradoxical network structures shown in is compatible with, for example, Pearl’s 2000 semantics
Figure 3, Panel (a). We would expect a positive relation- in terms of a “do-operator.” Here, an intervention on X
ship in the within-subjects network (e.g., typing faster is represented as Do(X = x), and the causal effect on Y
than your average leads to making more errors). Con- is formally expressed as E (Y | Do(X = x)). Pearl distin-
versely, we would expect a negative relationship in the guished this from the classical statistical association, in
between-subject network (e.g., people who type fast, on which no intervention is present, and we get the ordi-
average, generally make fewer spelling errors). This is nary regression E (Y | See(X = x)). This notation is use-
because people who type fast, on average, are likely to be ful here, because it can be used to show how different
more skilled in writing (e.g., a courtroom stenographer) kinds of causal manipulations, each at the intraindividual
and are less prone to make a lot of spelling errors, com- level, can produce a signal in either the between-subjects
pared to someone who types infrequently. Panel (b) of or the within-subjects network.
Figure 3 shows another example in which the structures Cashing out causal effects in terms of interventions
might differ (Hoffman 2015; provided by Hamaker 2017). is useful for understanding the intervention Do(X = x).
These network structures show that when people exert We can think of this in terms of a random shock to the
more physical activity than their average they likely expe- system, which sets X to value x at a particular time point
rience an elevated heart rate, while people who on average and evaluates the effect on another variable Y shortly
are often physically active likely have a lower average heart afterward. If we want to gauge this type of causal relation-
rate. Such a different effect depending on the level of anal- ship, we might look at the within-subjects VAR model.
ysis is well known in the statistical literature as Simpson’s Consider Hamaker’s (2012) example regarding typing
paradox (Simpson, 1951). errors: If a researcher forced a person to type very fast,
Interventionist accounts of causation. The different that researcher would need to evaluate the within-subject
ways of thinking about the effects of manipulations in data, which would show a positive association between
time-series models can be organized in terms recently typing speed and the number of errors. In this example,
developed from interventionist accounts of causation between-subjects data would be misleading because
(Woodward, 2005). According to Woodward, causation is individual differences would probably yield a negative
fleshed out in terms of interventions: X is a cause of Y correlation between speed and accuracy—faster typists
if an intervention (natural or experimental) on X leads are more likely to make less errors.
 MULTIVARIATE BEHAVIORAL RESEARCH 469

Interventions at the mean level. However, we can Five-Factor Model (neuroticism, extraversion, and con-
also think of a manipulation that sets X to value x scientiousness; McCrae & John, 1992) domains were
in a different way, for instance, by inducing a long- administered, as was an additional question that asked
term change in the system that leads it to converge on participants how much they had exercised since the pre-
X = x in expectation. To evaluate the effect of this type ceding measurement occasion. Sample 1 consisted of 26
of intervention, it is important to consider the behavior people providing 1323 observations in total, and Sample 2
of the system as it relates to the changes of the intercept consisted of 62 people providing a total of 2193 obser-
of X. When analyzing time-series data gathered in a rela- vations. Participants in Sample 1 answered questions
tively short time-span, the within-subjects VAR network three times per day, whereas participants in Sample 2
as discussed here cannot represent the relevant effects, answered questions five times per day. In both samples,
because it assumes stationarity. However, such effects the minimum time between measurements was 2 hours.
will be visible in the between-subjects network, which For more information about the samples and the specific
may thus contain important clues to the behavior of the questions asked, we refer readers to Mõttus et al. (2017).
system under potential changes in the intercept of one To obtain an easier and more interpretable example,
variable. In terms of Panel (b) of Figure 3, if we are inter- we first only analyzed questions aimed as measuring the
ested in the effect of changing someone structurally— extraversion trait and the question measuring exercise.
reducing the heart rate of a person on average—our pre- This leads to five variables of interest: questions pertain-
ferred source of hypothesis generation would likely stem ing to feeling outgoing, energetic, adventurous, or happy
from the between-subjects model, as the corresponding and the question measuring participants’ exercise habits.
within-subject model using the methods described in this We analyzed the data using the two-step multilevel VAR
paper only models deviations from the stationary mean. procedure as described in detail in the Section 2 of the
Such hypotheses could then be further investigated by supplementary materials. We used the mlVAR package,
using experimental design or lengthier longitudinal data version 0.4, for the estimation of this model. Because the
analysis. number of variables was small, we estimated the model
Many such examples can be envisioned, especially in using correlated temporal and contemporaneous random
the field of psychopathology. For instance, short-term effects. We ran the model separately for both samples and
deviations from the mean in abusing a substance might computed the fixed effects for the temporal, contempo-
not immediately develop tolerance or lead to one suf- raneous, and between-subjects networks. Correlations of
fering from work or life inferences, but a subject who the edge weights indicated that all three networks showed
abuses a substance on average over a long time period high correspondence between the two samples (temporal
might develop these problems (example based on vari- network: 0.82, contemporaneous network: 0.94, between-
ables used by Rhemtulla et al., 2016). A between-subjects subjects network: 0.70). Owing to the degree of replica-
network could similarly show that loneliness mediates the bility, we combined the two samples and estimated the
effect of losing a spouse on depressive symptoms (Fried model on the combined data.
et al., 2015) or highlight the possible effects of child- Results. Figure 4 shows the estimated fixed effects
hood trauma and urbanization on psychotic symptoms of the temporal, contemporaneous, and between-
(Isvoranu et al., 2016, 2017)—both cases in which within- subjects network. In these figures, only significant edges
subjects networks based on short-term deviations from (α = 0.05) are shown. In the contemporaneous and
the average seem less applicable. This analysis is impor- between-subjects networks, an edge was retained if one
tant because it shows that, even though relevant causal of the two regressions on which the partial correlation is
interventions in psychology will typically operate at the based was significant (the so-called “or” rule; van Borkulo
intra-individual level, evidence for the effect of such inter- et al., 2014). These results are in line with the hypothetical
ventions may arise at either the within- or the between- example shown in Figure 2: People who exercised were
subjects level depending on the nature of the intervention. more energetic while exercising and less energetic after
exercising. In the between-subjects network, no relation-
ship between exercising, feeling energetic, and feeling
4. Empirical examples
adventurous was found. The between-subjects network,
however, showed a strong relationship between feeling
4.1. Reanalysis of Mõttus et al. (2017)
adventurous and exercising: People who, on average,
We reanalyzed the data of Mõttus et al. (2017) to provide exercised more also felt, on average, more adventurous.
an empirical example of the multilevel VAR methods This relationship was not present in the temporal net-
described above. These data consist of two independent work and much weaker in the contemporaneous network.
ESM samples, in which items tapping three of the five Also noteworthy is that people were less outgoing after
470 S. EPSKAMP ET AL.

Temporal Contemporaneous Between−subjects

Adventurous Adventurous Adventurous

Outgoing Outgoing Outgoing

Energetic Energetic Energetic

Exercise Exercise Exercise

Happy Happy Happy

Maximum: 0.2 Maximum: 0.5 Maximum: 0.5

Figure . The estimated fixed effects of the three network structures obtainable in multilevel VAR. The model is based on ESM data of 
people providing a total of  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.

Temporal Contemporaneous

Adventurous Adventurous

Outgoing Outgoing

Energetic Energetic
Exercise Exercise

Happy Happy
Maximum: 0.18 Maximum: 0.01

Figure . 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.

exercising. Figure 5 shows the standard deviation of the network but not in the other networks. Finally, there was
random effects in the temporal and contemporaneous a between-subjects connection between exercising and
networks. The largest individual differences in the tem- feeling self-disciplined: People who, on average, exercised
poral network were found in the auto-regressions, and more also felt, on average, more self-disciplined.
the largest individual differences in the contemporaneous
network were found in relationship between exercising
and feeling energetic. 4.2. Reanalysis of Bringmann et al. (2013)
In addition to using only the extraversion and exer- To showcase additional information that can be obtained
cise items, we also ran the model on all 17 administered using the GGM model, we reanalyzed the data set used
items in the data set. In this analysis, we used orthogo- and made publicly available by Bringmann et al. (2013),
nal random effects to estimate the model because corre- which has been collected by Geschwind et al. (2011).
lated random effects cannot be estimated with such a large This data set contains ESM measures of 129 participants,
number of variables. Figure 6 shows the estimated fixed which was collected in two periods over 6 days each: a
effects of the three network structures; it can be seen that baseline period and a posttreatment period (mindfulness
indicators of the three traits tend to cluster together in all treatment and a control group). Participants answered
three networks. Regarding the node exercise, we found 60 measurements per period. Similar to Figure 1 of
the same relationships between exercise, energetic, and Bringmann et al. (2013), we analyzed only the baseline
adventurous (also found in the previous example) in the data set on the six items selected by Bringmann et al.
larger networks. Furthermore, we noted that exercising (2013). We estimated the networks using three model-
was connected to feeling angry in the between-subjects ing frameworks discussed in Table 2. First, we analyze
 MULTIVARIATE BEHAVIORAL RESEARCH 471

Figure . The estimated fixed effects of the three network structures based on all  variables administered. Only significant edges are
shown. Legend:  = “Worried”;  = “Organized”;  = “Ambitious”;  = “Depressed”;  = “Outgoing”;  = “Self-Conscious”;  = “Self-
Disciplined”;  = “Energetic”;  = “Frustrated”;  = “Focused”;  = “Guilty”;  = “Adventurous”;  = “Happy”;  = “Control”;  =
“Achieved”;  = “Angry”;  = “Exercise.”

data using multilevel Bayesian estimation using Mplus final two analyses, we did not regress the first measure-
version 8 (model generated using the mlVAR package). ment of the day on the last measurement of the previous
We estimated correlated random effects for the temporal day, and removed all pairs of lagged and current variables
effects but only fixed effects for the contemporaneous that contained missing responses. The final sample size
effects (making these random led to slow convergence). was 5927 observations. Edges were retained if they were
The model was estimated using three chains that ran significant at the α = 0.05 level, or if 0 was not included
until convergence. Nights were handled by adding a row in the 95% credibility interval.
of missing values between consecutive days. Second, we Results. Figure 7 shows the resulting network struc-
analyzed the data using two-step multilevel VAR esti- tures, and shows that all three methods are mostly aligned.
mation as implemented in the mlVAR package, using an Unsurprisingly, the temporal networks are very similar
“and”-rule and estimating correlated random temporal to those reported by Bringmann et al. (2013).17 Both the
and contemporaneous effects. Finally, we estimated the
data using pooled and individual LASSO estimation  The networks differ because the estimation of temporal effects differs in that
using the graphicalVAR package, using γ = 0.25. In the measures are within-subjects centered and subject means are included as
Level  predictors.
472 S. EPSKAMP ET AL.

Figure . Reanalysis of the Geschwind et al. () data set used by Bringmann et al. (). (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.

temporal and contemporaneous network are in line with Bayesian multilevel estimation resulted in a sparser tem-
what would be expected under a unidimensional auto- poral network. This difference in sparsity is possibly
correlated latent variable model (many edges selected, because the multivariate Bayesian multilevel approach
low-rank structure, edges of expected sign) with the more accurately represents the uncertainties in parameter
exception of the positive temporal edge from “fearful” estimation, while two-step multilevel VAR estimates the
to “pleasant” in the two-step multilevel network (which model piecewise and pooled LASSO estimation does not
was not selected by the other methods). Of note is that take the multilevel structure into account. Remarkable is
 MULTIVARIATE BEHAVIORAL RESEARCH 473

the positive edge between “sad” and “relaxed” in the two- n = 1 time series), the graphical VAR (GVAR; Wild et al.,
step multilevel between-subjects network, which is based 2010) model generalizes the GGM to incorporate tempo-
on two significant positive Level 2 regression coefficients ral effects. We showed how two network structures can
(β = 0.202, p = 0.046 and β = 0.151, p = 0.036) where be obtained: a temporal network, which is a directed net-
the estimated between-subjects correlation is strongly work of regression coefficients between lagged and cur-
negative (−0.53). This edge is especially remarkable since rent variables, and a contemporaneous network, which is
both nodes are strongly connected to other nodes in the a GGM describing the relationships that remain after con-
network. The Bayesian multilevel between-subjects net- trolling for temporal effects. In temporally ordered data
work showed a similar positive edge between “cheerful” of multiple subjects (e.g., n > 1 time series), the natural
and “worry.” This is noteworthy because under a unidi- combination of cross-sectional and time-series data came
mensional factor model, we would not expect partial cor- by adding a third network structure: the between-subjects
relation coefficients to switch sign from marginal correla- network, which is a GGM that describes relationships
tion coefficients (Holland & Rosenbaum, 1986; van Bork, between the stationary means of subjects. We proposed
Grasman, & Waldorp, 2016). A possible way the partial two methods to estimate the three network structures: (1)
correlation coefficient switches sign is if it has been con- two-step multilevel estimation, which we implemented in
ditioned on one or more common effects between the the open source R package mlVAR, and (2) pooled and
two variables of interest (in this case, potentially “worry,” individual VAR model estimations using LASSO regu-
“pleasant,” or “fearful”). Of course, these effects must be larization, which we implemented in the open source R
interpreted with great care, especially given the high p- package graphicalVAR.
values; we did not control for multiple comparisons, and
the same edges are not retained in the other methods. Still,
5.1. Limitations and challenges
it is noteworthy that if this edge is weak or nonexistent, the
between-subjects structure is still not in line with a unidi- Multilevel estimation. The presented methods are not
mensional factor model. In such a factor model, “sad” and without problems and have several limitations. With
“relaxed”(which feature the most connections) would be regard to multilevel estimation, first, multivariate estima-
expected to have a strong negative edge between them (a tion of the multilevel VAR model is not yet feasible for
depression factor would lead to “sad” having a strong pos- larger data sets. As such, the proposed two-step multi-
itive factor loading and “relaxed” having a strong negative level VAR combines univariate models. Doing so, how-
factor loading). ever, means that not all parameters are in the same model.
In addition, univariate models do not readily provide esti-
mates of the contemporaneous networks, which must be
5. Discussion
estimated in a second step. Second, even when multivari-
We discussed the Gaussian graphical model (GGM; Lau- ate estimation is possible, it is still challenging to estimate
ritzen, 1996), an undirected network model of partial a multilevel model on contemporaneous networks due to
correlation coefficients, and discussed its utility in the the requirement of positive definite matrices. Third, when
analysis of psychological data sets. The GGM presents a more than approximately eight variables are measured,
promising exploratory data analysis tool that allows for estimating the multilevel models with correlated random
different levels of interpretation: (1) Edges in the GGM effects is no longer feasible in open source software. In
can be interpreted without reliance on a causal interpre- this case, orthogonal random effects can be used, which
tation and merely used to show which variables predict induce a level of parsimony that may not be substan-
each other. (2) Causal effects between variables result in tively plausible. Finally, even when orthogonal estima-
an edge, whereas the lack of a causal effect results in no tion is used, multilevel analysis runs very slowly in mod-
edge, except in the presence of latent variables or a com- els with more than 20 variables. As such, multilevel VAR
mon effect. The GGM can, therefore, be seen as hypoth- analysis of high-dimensional data sets is not yet feasible.
esis generating structures that highlight potential causal To this end, we discussed pooling within-subject centered
pathways. (3) Undirected models can be used and inter- data and estimating fixed-effects models using LASSO
preted as causal data-generating process and have been regularization (Abegaz & Wit, 2013). This performed on
used as such in several fields of research. par with multilevel estimation in higher sample sizes and
The GGM can readily be estimated on any data set allows researchers to scale up the analysis. However, indi-
that contains multiple observations of the same vari- vidual network estimation using separate VAR models
ables (e.g., multiple people in cross-sectional data or does not borrow information from other subjects and
multiple responses in time-series data). LASSO regular- performs poorly in low sample sizes. Promising devel-
ization methods perform especially well in estimating opments are new LASSO methods in which shrinkage
such a GGM structure. In temporally ordered data (e.g., from subject-specific parameters to their mean is attained
474 S. EPSKAMP ET AL.

through penalization rather than hierarchical modeling variables as well as binary and Gaussian variables
(Hastie et al., 2015). Future research should investigate the (Haslbeck & Waldorp, 2016b). Such models have yet to
utility of such models in estimating individual network be extended to time-series analysis, especially in sep-
structures that might differ in structure but borrow infor- arating temporal and contemporaneous effects as the
mation from other subjects in its estimation. GVAR model does. When data are continuous but not
VAR modeling assumptions. These limitations on normal, multiple reasons can (again) contribute to this.
the estimation methods come with more limitations in When the underlying process is normal but the mea-
the statistical models themselves. VAR modeling, espe- sured variables are on a transformed scale, transforming
cially, is not without problems and faces severe challenges data back to normal should offer a solution (Liu, Laf-
(Hamaker, Ceulemans, Grasman, & Tuerlinckx, 2015; ferty, & Wasserman, 2009), but when the process itself
Hamaker & Wichers, 2017). We made several assump- is nonnormal, such as skewed residuals, the entire mod-
tions that can be problematic. For instance, in character- eling framework does not correctly capture the likeli-
izing the likelihood of time-series data, we need to assume hood. Finally, multivariate normality assumes all rela-
that the conditional distribution of variables at time t tionships between variables are linear. When this is not
given time t − 1 are the same for all t. That raises two the case, the GGM and VAR model (which fit linear
distinct assumptions: (1) The difference in time between effects) will not properly describe the data. We encour-
measurements are roughly equal, and (2) the parameters age future researchers to focus on the problem of nor-
do not change over time. Equidistance in time is espe- mality and to develop new methods of overcoming these
cially important for the interpretation of temporal net- challenges.
works. Promising work is being done in this area where Interpretation. Finally, it should be noted that when
VAR networks can be estimated on nonequidistant data taking a causal interpretation of edges, all methods dis-
sets (Driver, Oud, & Voelkle, 2017; Oravecz, Tuerlinckx, cussed in this paper are exploratory in nature and can only
& Vandekerckhove, 2009; Oud & Jansen, 2000). The generate hypotheses—they do not confirm causal rela-
assumption of stationarity is needed to estimate structures tions. The analyses showcased in this paper can also be
when data are limited but might not be tenable especially used without relying on a causal interpretation and allow
in longer time series (Rovine & Walls, 2006). Promising researchers to obtain insights into the predictive rela-
time-varying estimation procedures are being developed tionships present in the data—regardless of theory with
(Bringmann et al., 2016; Haslbeck & Waldorp, 2016a), but respect to the data-generating model. Under the assump-
are not yet extended to the GVAR framework. Further- tions of multivariate normality, stationarity, and the Lag-
more, the interpretation of temporal coefficients when 1 factorization, the networks show how variables predict
represented as a network is not without discussion, and each other over time (temporal network), within time
several different methods for standardization exist (Bul- (contemporaneous network), and on average (between-
teel, Tuerlinckx, Brose, & Ceulemans, 2016; Schuurman, subjects network). Furthermore, during the thresholding
Ferrer, de Boer-Sonnenschein, & Hamaker, 2016a).18 of edges in the multilevel analyses, we did not apply a
Normality. Another particularly important assump- correction for multiple testing by default. We deliberately
tion made in this paper is that of multivariate normality. chose this because our aim was to present exploratory
Indeed, Equation (1) makes this assumption and all other hypothesis-generating structures, and not correcting for
equations follow from this. The assumption of normality multiple testing yields greater sensitivity.
is not without problems (Terluin, de Boer, & de Vet, 2016).
However, it is not always straightforward to deal with
6. Conclusion
these issues, because violations of normality may arise for
many different reasons. When data are not normally dis- This paper provides a methodological overview of how
tributed, then they cannot be represented properly using the GGM can be used in various different kinds of
only the means vector and variance-covariance matrix. As psychological data. The GGM can be used to map out
a result, the GGM does not properly characterize the joint unique variance in cross-sectional data or at the con-
likelihood function. When data are measured on a differ- temporaneous and between-subjects levels of time-series
ent scale (Stevens, 1946), a different graphical model can analysis. We contrasted this method to exploratory esti-
be used, such as the Ising model for binary data (Epskamp mation of causal models. While losing information on the
et al., in press; van Borkulo et al., 2014) or a mixed direction of effect, estimating GGMs offers an attractive
graphical model for categorical and Poisson-distributed alternative in that these models are uniquely identified,
well parameterized, closely related to causal models
 Westandardized every data set before analyzing and used the standardiza- and also offer exploratory insight on predictive effects
tion of Wild et al. () for temporal networks in n = 1 and pooled temporal
networks. GGMs are readily standardized by using partial correlation coeffi-
between observed variables. When the aim is to discover
cients (Equation ()), which have been used in all GGMs shown in this paper. psychological dynamics, the GGM can be used as a
 MULTIVARIATE BEHAVIORAL RESEARCH 475

hypothesis generating technique inspiring future research Asparouhov, T., Hamaker, E., & Muthén, B. (2017). Dynamic
or therapy directions (Epskamp et al., 2018; Kroeze et al., structural equation models. Structural Equation Mod-
2017). For example, an effect found in a cross-sectional eling: A Multidisciplinary Journal, 25(3), 359–388. doi:
10.1080/10705511.2017.1406803
analysis could inspire a time-series study, a contempora- Beltz, A. M., & Molenaar, P. C. (2016). Dealing with multi-
neous effect could inspire a shorter time-lag time-series ple solutions in structural vector autoregressive models.
study and a between-subjects effect could inspire lengthy Multivariate Behavioral Research, 51(2–3), 357–373. doi:
longitudinal studies. All network structures may inspire 10.1080/00273171.2016.1151333.
experimental design, or to gather a mixture of observa- Bolger, N., & Laurenceau, J. (2013). Intensive longitudinal meth-
ods. New York, NY, USA: Guilford.
tional and experimental data (Magliacane et al., 2017).
Borsboom, D., & Cramer, A. O. J. (2013). Network analysis:
The GGM thus provides a powerful addition to the An integrative approach to the structure of psychopathol-
exploratory toolbox in behavioral research. ogy. Annual Review of Clinical Psychology, 9, 91–121. doi:
10.1146/annurev-clinpsy-050212-185608.
Borsboom, D., Cramer, A. O. J., Schmittmann, V. D.,
Article Information Epskamp, S., & Waldorp, L. J. (2011). The small
world of psychopathology. PloS One, 6(11), e27407.
Conflict of Interest Disclosures: Each author signed a doi: 10.1371/[Link].0027407.
form for disclosure of potential conflicts of interest. No Bos, E. H., & Wanders, R. B. (2016). Group-level symptom net-
authors reported any financial or other conflicts of inter- works in depression. JAMA Psychiatry, 73(4), 411–411. doi:
est in relation to the work described. 10.1001/jamapsychiatry.2015.3103.
Bos, F. M., Snippe, E., de Vos, S., Hartmann, J. A., Simons,
Ethical Principles: The authors affirm having followed C. J. P., van der Krieke, L., .. Wichers, M. (2017). Can we
professional ethical guidelines in preparing this work. jump from cross-sectional to dynamic interpretations of
networks? Implications for the network perspective in psy-
These guidelines include obtaining informed consent
chiatry. Psychotherapy and Psychosomatics, 86, 185–177. doi:
from human participants, maintaining ethical treatment 10.1159/000453583.
and respect for the rights of human or animal partici- Brandt, P. T., & Williams, J. T. (2007). Multiple time series mod-
pants, and ensuring the privacy of participants and their els (vol. 148). Thousand Oaks, CA, USA: Sage Publications,
data, such as ensuring that individual participants cannot Inc.
be identified in reported results or from publicly available Bringmann, L. F., Hamaker, E. L., Vigo, D. E., Aubert, A.,
Borsboom, D., & Tuerlinckx, F. (2016). Changing dynam-
original or archival data. ics: Time-varying autoregressive models using generalized
Funding: This work was supported by Grant 406-11-066 additive modeling. Psychological Methods, 22(3), 409–425.
doi: 10.1037/met0000085.
from the NWO.
Bringmann, L. F., Lemmens, L. H., Huibers, M. J., Bors-
Role of the Funders/Sponsors: None of the funders or boom, D., & Tuerlinckx, F. (2015). Revealing the
sponsors of this research had any role in the design dynamic network structure of the Beck Depression
Inventory-II. Psychological Medicine, 45(04), 747–757. doi:
and conduct of the study; collection, management, anal- 10.1017/S0033291714001809.
ysis, and interpretation of data; preparation, review, or Bringmann, L. F., Vissers, N., Wichers, M., Geschwind,
approval of the manuscript; or decision to submit the N., Kuppens, P., Peeters, F., .. Tuerlinckx, F. (2013). A
manuscript for publication. network approach to psychopathology: New insights
into clinical longitudinal data. PLOS ONE, 8(4),
Acknowledgments: The authors would like to thank e60188.
Laura Bringmann, Noémi Schuurman, Oisín Ryan, and Brose, A., Lindenberger, U., & Schmiedek, F. (2013). Affective
Ellen Hamaker for helpful tips and invigorating discus- states contribute to trait reports of affective well-being. Emo-
sions, and Katharina Jorgensen for valuable comments on tion, 13(5), 940.
Brown, T. A. (2014). Confirmatory factor analysis for applied
earlier versions of this paper. An earlier version of this research. London, UK: Guilford Publications.
paper has been adapted as a chapter in the dissertation Bulteel, K., Tuerlinckx, F., Brose, A., & Ceulemans, E. (2016).
of the main author (Epskamp, 2017b). Using raw var regression coefficients to build networks can
be misleading. Multivariate Behavioral Research, 51(2–3),
330–344. doi: 10.1080/00273171.2016.1150151.
Chatfield, C. (2016). The analysis of time series: An introduction.
References Boca Raton, FL, USA: CRC press.
Chen, G., Glen, D. R., Saad, Z. S., Hamilton, J. P., Thoma-
Abegaz, F., & Wit, E. (2013). Sparse time series chain graphi-
son, M. E., Gotlib, I. H., .. Cox, R. W. (2011). Vec-
cal models for reconstructing genetic networks. Biostatistics,
tor autoregression, structural equation modeling, and
14(3), 586–599. doi: 10.1093/biostatistics/kxt005.
their synthesis in neuroimaging data analysis. Comput-
Abegaz, F., & Wit, E. (2015). SparseTSCGM: Sparse time series
ers in Biology and Medicine, 41(12), 1142–1155. doi:
chain graphical models (R package version 2.2). Retrieved
10.1016/[Link].2011.09.004.
from [Link]
476 S. EPSKAMP ET AL.

Chen, J., & Chen, Z. (2008). Extended Bayesian informa- Epskamp, S., Cramer, A., Waldorp, L., Schmittmann, V. D.,
tion criteria for model selection with large model spaces. & Borsboom, D. (2012). qgraph: Network visualizations
Biometrika, 95(3), 759–771. doi: 10.1093/biomet/asn034. of relationships in psychometric data. Journal of Statistical
Costantini, G., Epskamp, S., Borsboom, D., Perugini, M., Mõt- Software, 48(1), 1–18. doi: 10.18637/jss.v048.i04.
tus, R., Waldorp, L. J., .. Cramer, A. O. J. (2015). State of the Epskamp, S., Deserno, M. K., & Bringmann, L. F. (2017b).
art personality research: A tutorial on network analysis of mlVAR: Multi-level vector autoregression (R pack-
personality data in R. Journal of Research in Personality, 54, age version 0.4). Retrieved from [Link]
13–29. doi: 10.1016/[Link].2014.07.003. [Link]/package=mlVAR.
Cramer, A. O., van Borkulo, C. D., Giltay, E. J., van der Maas, Epskamp, S., & Fried, E. I. (in press). A tutorial on estimating
H. L., Kendler, K. S., Scheffer, M., .. Borsboom, D. (2016). regularized psychological networks. Psychological Methods.
Major depression as a complex dynamic system. PloS One, DOI: 10.1037/met0000167
11(12), e0167490. Epskamp, S., Kruis, J., & Marsman, M. (2017c). Estimating psy-
Cramer, A. O. J., Sluis, S., Noordhof, A., Wichers, M., chopathological networks: Be careful what you wish for.
Geschwind, N., Aggen, S. H., .. Borsboom, D. (2012). PloS One, 12(6), e0179891.
Dimensions of normal personality as networks in search Epskamp, S., Maris, G., Waldorp, L., & Borsboom, D. (in press).
of equilibrium: You can’t like parties if you don’t like peo- Network psychometrics. In P. Irwing, D. Hughes, & T. Booth
ple. European Journal of Personality, 26(4), 414–431. doi: (Eds.), Handbook of psychometrics. New York, NY, USA:
10.1002/per.1866. Wiley. Retrieved from [Link]/abs/1609.02818.
Cramer, A. O. J., Waldorp, L., van der Maas, H., & Bors- Epskamp, S., Rhemtulla, M., & Borsboom, D. (2017d). General-
boom, D. (2010). Comorbidity: A Network Perspective. ized network psychometrics: Combining network and latent
Behavioral and Brain Sciences, 33(2–3), 137–150. doi: variable models. Psychometrika, 82(4), 904–927.
10.1017/S0140525X09991567. Epskamp, S., van Borkulo, C. D., Isvoranu, A. M., Ser-
Curran, P. J., & Bauer, D. J. (2011). The disaggregation of within- vaas, M. N., van der Veen, D. C., Riese, H., ... Cramer,
person and between-person effects in longitudinal models A. O. J. (2018). Personalized network modeling in psy-
of change. Annual Review of Psychology, 62, 583–619. doi: chopathology: The importance of contemporaneous and
10.1146/[Link].093008.100356. temporal connections. Clinical Psychological Science. doi:
Dalege, J., Borsboom, D., van Harreveld, F., van den Berg, H., 10.1177/2167702617744325.
Conner, M., & van der Maas, H. L. (2016). Toward a formal- Foygel, R., & Drton, M. (2010). Extended Bayesian information
ized account of attitudes: The causal attitude network (can) criteria for Gaussian graphical models. Advances in Neural
model. Psychological Review, 123(1), 2–22. Information Processing Systems, 23, 2020–2028. Retrieved
Driver, C. C., Oud, J. H. L., & Voelkle, M. C. (2017). Con- from [Link]
tinuous time structural equation modelling with r pack- Fried, E. I., Bockting, C., Arjadi, R., Borsboom, D., Amshoff,
age ctsem. Journal of Statistical Software, 77, 1–35. doi: M., Cramer, O. J., .. Stroebe, M. (2015). From loss to loneli-
10.18637/jss.v077.i05. ness: The relationship between bereavement and depressive
Drton, M., & Maathuis, M. H. (2017). Structure learning in symptoms. Journal of Abnormal Psychology, 124(2), 256–
graphical modeling. Annual Review of Statistics and Its 265. doi: 10.1037/abn0000028.
Application, 4, 365–393. doi: 10.1146/annurev-statistics- Fried, E. I., Epskamp, S., Nesse, R. M., Tuerlinckx, F., & Bors-
060116-053803. boom, D. (2016). What are ‘good’ depression symptoms?
Eichler, M. (2007). Granger causality and path diagrams for Comparing the centrality of DSM and non-DSM symptoms
multivariate time series. Journal of Econometrics, 137(2), of depression in a network analysis. Journal of Affective Dis-
334–353. doi: 10.1016/[Link].2005.06.032. orders, 189, 314–320. doi: 10.1016/[Link].2015.09.005.
Epskamp, S. (2016). Regularized Gaussian psycholog- Friedman, J. H., Hastie, T., & Tibshirani, R. (2008). Sparse
ical networks: Brief report on the performance inverse covariance estimation with the graphical lasso. Bio-
of extended BIC model selection. Retrieved from statistics, 9(3), 432–441. doi: 10.1093/biostatistics/kxm045.
[Link] Friedman, J. H., Hastie, T., & Tibshirani, R. (2014). glasso:
Epskamp, S. (2017a). Discussion: The road ahead. In Net- Graphical lasso-estimation of Gaussian graphical models
work psychometrics, 237–248. Retrieved from http:// (R package version 1.8). Retrieved from [Link]
[Link]/dissertation/[Link]. [Link]/package=glasso.
Epskamp, S. (2017b). Discovering psychological dynam- Gates, K. M., & Molenaar, P. C. (2012). Group search algo-
ics. In Network psychometrics, 85–114. Retrieved from rithm recovers effective connectivity maps for individuals
[Link] in homogeneous and heterogeneous samples. NeuroImage,
Epskamp, S. (2017c). graphicalVAR: Graphical VAR for experi- 63(1), 310–319. doi: 10.1016/[Link].2012.06.026.
ence sampling data (R package version 0.1.6). Retrieved from Gates, K. M., Molenaar, P. C., Hillary, F. G., Ram, N., &
[Link] Rovine, M. J. (2010). Automatic search for fMRI connec-
Epskamp, S. (2018). Preliminary simulations on the interpreta- tivity mapping: An alternative to Granger causality test-
tion of cross-sectional Gaussian graphical models. PsyArXiv ing using formal equivalences among SEM path modeling,
preprint, 54xrs. doi: 10.17605/[Link]/54XRS. var, and unified SEM. Neuroimage, 50(3), 1118–1125. doi:
Epskamp, S., Borsboom, D., & Fried, E. I. (2017a). Estimat- 10.1016/[Link].2009.12.117.
ing psychological networks and their accuracy: A tuto- Gelman, A., & Hill, J. (2006). Data analysis using regression and
rial paper. Behavior Research Methods, 50(1), 195–212. doi: multilevel/hierarchical models. New York, NY, USA: Cam-
10.3758/s13428-017-0862-1. bridge University Press.
 MULTIVARIATE BEHAVIORAL RESEARCH 477

Geschwind, N., Peeters, F., Drukker, M., van Os, J., & Wich- Hoffman, L. (2015). Longitudinal analysis: Modeling within-
ers, M. (2011). Mindfulness training increases momentary person fluctuation and change. New York, NY, USA: Rout-
positive emotions and reward experience in adults vulner- ledge.
able to depression: a randomized controlled trial. Journal Hoffman, L., & Stawski, R. S. (2009). Persons as contexts: Eval-
of Consulting and Clinical Psychology, 79(5), 618–628. doi: uating between-person and within-person effects in longi-
10.1037/a0024595. tudinal analysis. Research in Human Development, 6(2–3),
Golino, H. F., & Epskamp, S. (2017). Exploratory graph analy- 97–120. doi: 10.1080/15427600902911189.
sis: A new approach for estimating the number of dimen- Holland, P. W., & Rosenbaum, P. R. (1986). Conditional asso-
sions in psychological research. PlosOne, 2(6), e0174035. ciation and unidimensionality in monotone latent variable
doi: 10.1371/[Link].0174035. models. The Annals of Statistics, 14(4), 1523–1543. doi:
Granger, C. W. J. (1969). Investigating causal relations by econo- 10.1214/aos/1176350174.
metric models and cross-spectral methods. Econometrica: Ising, E. (1925). Beitrag zur theorie des ferromagnetismus.
Journal of the Econometric Society, 37(3), 424–438. Zeitschrift für Physik A Hadrons and Nuclei, 31(1), 253–258.
Guttman, L. (1938). A note on the derivation of formu- doi: 10.1007/BF02980577.
lae for multiple and partial correlation. The Annals of Isvoranu, A. M., Borsboom, D., van Os, J., & Guloksuz, S. (2016).
Mathematical Statistics, 9(4), 305–308. Retrieved from A network approach to environmental impact in psychotic
[Link] disorders: Brief theoretical framework. Schizophrenia Bul-
Hallquist, M., & Wiley, J. (2017). MplusAutomation: Automat- letin, 42(4), 870–873. doi: 10.1093/schbul/sbw049.
ing Mplus model estimation and interpretation [Computer Isvoranu, A. M., van Borkulo, C. D., Boyette, L., Wigman, J. T.
software manual] (R package version 0.7). Retrieved from W., Vinkers, C. H., Borsboom, D., & GROUP Investigators,
[Link] (2017). A network approach to psychosis: Pathways between
Hamaker, E., Ceulemans, E., Grasman, R., & Tuerlinckx, F. childhood trauma and psychotic symptoms. Schizophrenia
(2015). Modeling affect dynamics: State of the art and Bulletin, 43(1), 187–196. doi: 10.1093/schbul/sbw055.
future challenges. Emotion Review, 7(4), 316–322. doi: Kalisch, M., & Bühlmann, P. (2007). Estimating
10.1177/1754073915590619. high-dimensional directed acyclic graphs with
Hamaker, E. L. (2012). Why researchers should think the pc-algorithm. Journal of Machine Learn-
“within-person”: A paradigmatic rationale. In M. ing Research, 8(Mar), 613–636. Retrieved from
R. Mehl, & T. S. Conner (Eds.), Handbook of [Link]
research methods for studying daily life (pp. 43– Kalisch, M., Mächler, M., Colombo, D., Maathuis, M. H., &
61). New York, NY: Guilford Press. Retrieved from Bühlmann, P. (2012). Causal inference using graphical
[Link] models with the R package pcalg. Journal of Statistical Soft-
Hamaker, E. L. (2017). A brief history of dynamic modeling in ware, 47(11), 1–26. doi: 10.18637/jss.v047.i11.
psychology. 82th Annual Meeting of the Psychometric Soci- Kalisch, M., Maechler, M., & Colombo, D. (2017). pcalg: Esti-
ety (IMPS), Zurich: Switzerland. mation of CPDAG/PAG and causal inference using the
Hamaker, E. L., & Grasman, R. P. (2014). To center or not IDA algorithm. (R package version 2.5-0). Retrieved from
to center? Investigating inertia with a multilevel autore- [Link]
gressive model. Frontiers in Psychology, 5, 1492. doi: Kaplan, D. (2000). Structural equation modeling: Foundations
10.3389/fpsyg.2014.01492. and extensions. CA, USA: Sage, Thousand Oaks.
Hamaker, E. L., & Wichers, M. (2017). No time like the present. Kim, C.-J., & Nelson, C. R. (1999). State-space models with
Current Directions in Psychological Science, 26(1), 10–15. regime switching: Classical and Gibbs-sampling approaches
doi: 10.1177/0963721416666518. with applications. Cambridge, MA, USA: The MIT Press.
Hamilton, J. D. (1994). Time series analysis (vol. 2). Princeton, Klippel, A., Viechtbauer, W., Reininghaus, U., Wigman, J., van
NJ, USA: Princeton university press. Borkulo, C., MERGE, .. Wichers, M. (2017). The cascade of
Harvey, A. C. (1990). Forecasting, structural time series models stress: A network approach to explore differential dynamics
and the Kalman filter. Cambridge, UK: Cambridge univer- in populations varying in risk for psychosis. Schizophrenia
sity press. Bulletin, 44(2), 328–337. doi: 10.1093/schbul/sbx037.
Haslbeck, J. M. B., & Waldorp, L. J. (2016a). mgm: Structure esti- Koller, D., & Friedman, N. (2009). Probabilistic graphical mod-
mation for time-varying mixed graphical models in high- els: Principles and techniques. Cambridge, MA, USA: MIT
dimensional data. arXiv preprint, page arXiv:1510.06871. Press.
Haslbeck, J. M. B., & Waldorp, L. J. (2016b). Structure estima- Kossakowski, J. J., Epskamp, S., Kieffer, J. M., van Borkulo, C.
tion for mixed graphical models in high dimensional data. D., Rhemtulla, M., & Borsboom, D. (2015). The applica-
arXiv preprint, page arXiv:1510.05677. tion of a network approach to health-related quality of life
Hastie, T., Tibshirani, R., & Wainwright, M. (2015). Statistical (HRQoL): Introducing a new method for assessing hrqol in
learning with sparsity: The Lasso and generalizations. Boca healthy adults and cancer patient. Quality of Life Research,
Raton, FL, USA: CRC Press. 25, 781–792. doi: 10.1007/s11136-015-1127-z.
Hayduk, L. A. (2009). Finite feedback cycling in structural equa- Krämer, N., Schäfer, J., & Boulesteix, A.-L. (2009). Regularized
tion models. Structural Equation Modeling, 16(4), 658–675. estimation of large-scale gene association networks using
doi: 10.1080/10705510903206030. graphical Gaussian models. BMC Bioinformatics, 10(1), 1–
Heiser, W. J. (2017). Early psychometric contributions to Gaus- 24. doi: 10.1186/1471-2105-10-384.
sian graphical modeling: A tribute to Louis Guttman (1916- Kroeze, R., Van Veen, D., Servaas, M., Bastiaansen, J., Oude
1987). Zurich, Switzerland: 82th Annual Meeting of the Psy- Voshaar, R., Borsboom, D., .. Riese, H. (2017). Personal-
chometric Society (IMPS). ized feedback on symptom dynamics of psychopathology:
478 S. EPSKAMP ET AL.

A proof-of-principle study. Journal for Person-Oriented Molenaar, P. C. (2004). A manifesto on psychology as idio-
Research, 3(1), 1–10. doi: 10.17505/jpor.2017.01. graphic science: Bringing the person back into scientific
Kruis, J., & Maris, G. (2016). Three representations of the ising psychology, this time forever. Measurement, 2(4), 201–218.
model. Scientific Reports, 6, 34175. doi: 10.1207/s15366359mea0204_1.
Lane, S., Gates, K., Molenaar, P., Hallquist, M., & Pike, H. Mõttus, R., Epskamp, S., & Francis, A. (2017). Within-and
(2016). Gimme: Group iterative multiple model estimation. between individual variability of personality characteristics
(R package version 0.1-7). Retrieved from [Link] and physical exercise. Journal of Research in Personality, 69,
[Link]/package=gimme. 139–148. doi: 10.1016/[Link].2016.06.017.
Lauritzen, S. L. (1996). Graphical models. Oxford, UK: Claren- Murphy, K. P. (2012). Machine learning: A probabilistic perspec-
don Press. tive. Cambridge, MA, USA: MIT Press.
Liu, H., Lafferty, J. D., & Wasserman, L. (2009). The nonpara- Muthén, L. K., & Muthén, B. O. (2017). Mplus user’s
normal: Semiparametric estimation of high dimen- guide: Statistical analysis with latent variables. Statis-
sional undirected graphs. The Journal of Machine tical analysis with latent variables. Version. Retrieved
Learning Research, 10, 2295–2328. Retrieved from from [Link]
[Link] /MplusUserGuideVer_8.pdf.
Lord, F. M., Novick, M. R., & Birnbaum, A. (1968). Statistical Myin-Germeys, I., Oorschot, M., Collip, D., Lataster, J.,
theories of mental test scores. Oxford, UK: Addison-Wesley. Delespaul, P., & van Os, J. (2009). Experience sampling
Lütkepohl, H. (2005). New introduction to multiple time series research in psychopathology: Opening the black box of
analysis. Berlin, Germany: Springer-Verlag. daily life. Psychological Medicine, 39(09), 1533–1547. doi:
MacCallum, R. C., Wegener, D. T., Uchino, B. N., & Fabrigar, 10.1017/S0033291708004947.
L. R. (1993). The problem of equivalent models in applica- Oravecz, Z., Tuerlinckx, F., & Vandekerckhove, J. (2009). A
tions of covariance structure analysis. Psychological Bulletin, hierarchical Ornstein–Uhlenbeck model for continuous
114(1), 185–199. doi: 10.1037/0033-2909.114.1.185. repeated measurement data. Psychometrika, 74(3), 395–418.
Magliacane, S., van Ommen, T., Claassen, T., Bongers, S., Ver- doi: 10.1007/s11336-008-9106-8.
steeg, P., & Mooij, J. M. (2017). Causal transfer learning. Oud, J. H., & Jansen, R. A. (2000). Continuous time state space
arXiv preprint arXiv:1707.06422. modeling of panel data by means of SEM. Psychometrika,
Marchetti, G. M., Drton, M., & Sadeghi, K. (2015). ggm: 65(2), 199–215. doi: 10.1007/BF02294374.
Functions for graphical Markov models (R pack- Pan, J., Ip, E. H., & Dubé, L. (2017). An alternative to post
age version 2.3). Retrieved from [Link] hoc model modification in confirmatory factor analysis:
[Link]/package=ggm. The Bayesian lasso. Psychological Methods, 22(4), 687. doi:
Marsman, M., Borsboom, D., Kruis, J., Epskamp, S., van Bork, 10.1037/met0000112.
R., Waldorp, L., ... Maris, G. (2017a). An introduction to net- Pe, M. L., Kircanski, K., Thompson, R. J., Bringmann, L.
work psychometrics: Relating ising network models to item F., Tuerlinckx, F., Mestdagh, M., .. Jonides, J. (2015).
response theory models. Multivariate Behavioral Research, Emotion-network density in major depressive disor-
53(1), 15–35. doi: 10.1080/00273171.2017.1379379. der. Clinical Psychological Science, 3(2), 292–300. doi:
Marsman, M., Maris, G., Bechger, T., & Glas, C. (2015). 10.1177/2167702614540645.
Bayesian inference for low-rank ising networks. Scientific Pearl, J. (2000). Causality: Models, reasoning, and inference. New
Reports, 5(9050), 1–7. doi: 10.1038/srep09050. York, NY: Cambridge University Press.
Marsman, M., Waldorp, L., & Maris, G. (2017b). A note on Pourahmadi, M. (2011). Covariance estimation: The glm and
large-scale logistic prediction: Using an approximate graph- regularization perspectives. Statistical Science, 26(3), 369–
ical model to deal with collinearity and missing data. Behav- 387. doi: 10.1214/11-STS358.
iormetrika, 44(2), 513–534. doi: 10.1007/s41237-017-0024- Quax, R., Kandhai, D., & Sloot, P. M. (2013). Information dissi-
x. pation as an early-warning signal for the Lehman brothers
McCrae, R. R., & John, O. P. (1992). An introduction to the collapse in financial time series. Scientific Reports, 3, 1898.
five-factor model and its applications. Journal of Personality, doi: 10.1038/srep01898.
60(2), 175–215. doi: 10.1111/j.1467-6494.1992.tb00970.x. R Core Team, (2017). R: A Language and Environment for Sta-
McNally, R. J., Robinaugh, D. J., Wu, G. W., Wang, L., Deserno, tistical Computing. Vienna, Austria: R Foundation for Sta-
M. K., & Borsboom, D. (2015). Mental disorders as causal tistical Computing. Retrieved from [Link]
systems a network approach to posttraumatic stress dis- org/.
order. Clinical Psychological Science, 3(6), 836–849. doi: Rhemtulla, M., Fried, E. I., Aggen, S. H., Tuerlinckx, F.,
10.1177/2167702614553230. Kendler, K. S., & Borsboom, D. (2016). Network
Meehl, P. E. (1990). Why summaries of research on psychologi- analysis of substance abuse and dependence symp-
cal theories are often uninterpretable. Psychological Reports, toms. Drug and Alcohol Dependence, 161, 230–237. doi:
66(1), 195–244. doi: 10.2466/pr0.1990.66.1.195. 10.1016/[Link].2016.02.005.
Meinshausen, N., & Bühlmann, P. (2006). High- Rigdon, E. E. (1995). A necessary and sufficient identifi-
dimensional graphs and variable selection with the cation rule for structural models estimated in practice.
lasso. The Annals of Statistics, 34(3), 1436–1462. doi: Multivariate Behavioral Research, 30(3), 359–383. doi:
10.1214/009053606000000281. 10.1207/s15327906mbr3003_4.
Mohammadi, A., & Wit, E. C. (2015). BDgraph: An R pack- Rosmalen, J. G., Wenting, A. M., Roest, A. M., de Jonge, P., &
age for Bayesian structure learning in graphical mod- Bos, E. H. (2012). Revealing causal heterogeneity using time
els. arXiv preprint, page arXiv:1501.05108. Retrieved from series analysis of ambulatory assessments: application to the
[Link] association between depression and physical activity after
 MULTIVARIATE BEHAVIORAL RESEARCH 479

myocardial infarction. Psychosomatic Medicine, 74(4), 377– Stevens, S. S. (1946). On the theory of scales of mea-
386. doi: 10.1097/PSY.0b013e3182545d47. surement. Science, New Series, 103(2684), 677–680. doi:
Rothman, A. J., Levina, E., & Zhu, J. (2010). Sparse multi- 10.1126/science.103.2684.677.
variate regression with covariance estimation. Journal of Terluin, B., de Boer, M. R., & de Vet, H. C. (2016). Differences
Computational and Graphical Statistics, 19(4), 947–962. doi: in connection strength between mental symptoms might be
10.1198/jcgs.2010.09188. explained by differences in variance: Reanalysis of network
Rovine, M. J., & Walls, T. A. (2006). Multilevel autoregres- data did not confirm staging. PloS One, 11(11), e0155205.
sive modeling of interindividual differences in the sta- doi: 10.1371/[Link].0155205.
bility of a process. In Theodore A. Walls & Joseph L. Tibshirani, R. (1996). Regression shrinkage and selection
Schafer (Eds.), Models for Intensive Longitudinal Data via the lasso. Journal of the Royal Statistical Society.
(pp. 124–147). Oxford, UK: Oxford University Press. doi: Series B (Methodological), 58, 267–288. Retrieved from
10.3389/fpsyg.2017.00262. [Link]
Schafer, J., Opgen-Rhein, R., Zuber, V., Ahdesmaki, M., van Bork, R., Grasmans, R. P. P. P., & Waldorp, L. J. (2016).
Silva, A. P. D., & Strimmer, K. (2017). corpcor: Effi- Unidimensional factor models imply weaker partial corre-
cient estimation of covariance and (partial) correlation (R lations than zero-order correlations. arXiv preprint, in press.
package version 1.6.9). Retrieved from [Link] page arXiv:1610.03375.
[Link]/package=corpcor. van Borkulo, C. D., Borsboom, D., Epskamp, S., Blanken,
Schmiedek, F., Lövdén, M., & Lindenberger, U. (2010). Hundred T. F., Boschloo, L., Schoevers, R. A., .. Waldorp, L.
days of cognitive training enhance broad cognitive abilities J. (2014). A new method for constructing networks
in adulthood: Findings from the cogito study. Frontiers in from binary data. Scientific Reports, 4(5918), 1–10. doi:
Aging Neuroscience, 2, 27. doi: 10.3389/fnagi.2010.00027. 10.1038/srep05918.
Schmittmann, V. D., Cramer, A. O. J., Waldorp, L. J., Epskamp, van Borkulo, C. D., Boschloo, L., Borsboom, D., Penninx, B.
S., Kievit, R. A., & Borsboom, D. (2013). Deconstruct- W. J. H., Waldorp, L. J., & Schoevers, R. A. (2015). Asso-
ing the construct: A network perspective on psychologi- ciation of symptom network structure with the course
cal phenomena. New Ideas in Psychology, 31(1), 43–53. doi: of depression. JAMA Psychiatry, 72(12), 1219–1226. doi:
10.1016/[Link].2011.02.007. 10.1001/jamapsychiatry.2015.2079.
Schuurman, N. K. (2016). Measurement error and person- van der Maas, H. L., Dolan, C. V., Grasman, R. P., Wicherts, J. M.,
specific reliabilities in multilevel autoregressive model- Huizenga, H. M., & Raijmakers, M. E. (2006). A dynamical
ing. In Multilevel autoregressive modeling in psychology: model of general intelligence: The positive manifold of intel-
Snags and solutions, 139–178. Retrieved from http:// ligence by mutualism. Psychological Review, 113(4), 842–
[Link]/NKSchuurman_dissertation.pdf. 861. doi: 10.1037/0033-295X.113.4.842.
Schuurman, N. K., Ferrer, E., de Boer-Sonnenschein, M., & Wigman, J., van Os, J., Borsboom, D., Wardenaar, K., Epskamp,
Hamaker, E. L. (2016a). How to compare cross-lagged asso- S., Klippel, A., .. Wichers, M. (2015). Exploring the
ciations in a multilevel autoregressive model. Psychological underlying structure of mental disorders: Cross-
Methods, 21(2), 206. doi: 10.1037/met0000062. diagnostic differences and similarities from a network
Schuurman, N. K., Grasman, R. P. P. P., & Hamaker, E. perspective using both a top-down and a bottom-up
L. (2016b). A comparison of inverse-Wishart prior approach. Psychological Medicine, 45(11), 2375–2387. doi:
specifications for covariance matrices in multilevel autore- 10.1017/S0033291715000331.
gressive models. Multivariate Behavioral Research, 51(2–3), Wild, B., Eichler, M., Friederich, H.-C., Hartmann, M.,
185–206. doi: 10.1080/00273171.2015.1065398. Zipfel, S., & Herzog, W. (2010). A graphical vector
Schuurman, N. K., Houtveen, J. H., & Hamaker, E. L. (2015). autoregressive modeling approach to the analysis of elec-
Incorporating measurement error in n= 1 psychological tronic diary data. BMC Medical Research Methodology,
autoregressive modeling. Frontiers in Psychology, 6. doi: 10(1), 28. doi: 10.1186/1471-2288-10-28.
10.3389/fpsyg.2015.01038. Witten, D. M., Friedman, J. H., & Simon, N. (2011). New insights
Scutari, M. (2010). Learning Bayesian networks with the and faster computations for the graphical lasso. Journal of
bnlearn R package. Journal of Statistical Software, 35(3), 1– Computational and Graphical Statistics, 20(4), 892–900. doi:
22. doi: 10.18637/jss.v035.i03. 10.1198/jcgs.2011.11051a.
Shumway, R. H., & Stoffer, D. S. (2010). Time series analysis Woodward, J. (2005). Making things happen: A theory of
and its applications: With R examples. New York, NY, USA: causal explanation. Oxford, UK: Oxford University
Springer Science & Business Media. doi: 10.1007/978-3- Press. Retrieved from [Link]
319-52452-8. 611.
Simpson, E. H. (1951). The interpretation of interaction in Wright, S. (1921). Correlation and causation. Journal of Agricul-
contingency tables. Journal of the Royal Statistical Society. tural Research, 20(7), 557–585.
Series B (Methodological), 13(2), 238–241. Retrieved from Yuan, M., & Lin, Y. (2007). Model selection and estimation
[Link] in the Gaussian graphical model. Biometrika, 94(1), 19–35.
Snippe, E., Viechtbauer, W., Geschwind, N., Klippel, A., De doi: 10.1093/biomet/asm018.
Jonge, P., & Wichers, M. (2017). The impact of treatments Zhao, T., Li, X., Liu, H., Roeder, K., Lafferty, J., & Wasserman, L.
for depression on the dynamic network structure of mental (2015). huge: High-dimensional undirected graph estimation
states: Two randomized controlled trials. Scientific Reports, (R package version 1.2.7). Retrieved from [Link]
7, 46523. doi: 10.1038/srep46523. [Link]/package=huge
480 S. EPSKAMP ET AL.

Appendix: Glossary of terms

Term Explanation

Undirected network A network model in which nodes are connected by edges (also termed links) without arrowheads.
Directed network A network model in which nodes are connected by edges with arrowheads, assumed to display causal effects or
temporal prediction.
Gaussian graphical model An undirected network model in which observed variables are represented with nodes. Nodes are connected with
an edge if two variables are not independent after conditioning on all other observed variables. Edges are
parameterized by using partial correlation coefficients.
Causal model A causal model of observed and unobserved variables that is assumed to generate the data.
Directed acyclic graph A directed network in which one node does not eventually point to itself.
Within-subjects network A network model explaining within-subject (co)variation from the stationary mean.
Between-subjects network A network model explaining (co)variation between stationary means of different persons.
Cross-sectional network A network model estimated on cross-sectional data. Can be shown to be a blend of the within-subjects and
between-subjects networks. Can be interpreted as representative of within-subjects or between-subjects network
based on the way in which data are gathered.
Vector auto-regression (VAR) Multivariate regression of a set of variables on previous realizations of that set of variables.
Temporal network A within-subject network model of effects between different measurement occasions, showing temporal
prediction or potential causal pathways.
Contemporaneous network A within-subject undirected network model of effects between variables in the same measurement occasion, after
taking temporal effects into account.